Page MenuHomeHEPForge
Files F10881204

No OneTemporary
Actions

View Options

diff --git a/include/Rivet/Makefile.am b/include/Rivet/Makefile.am
--- a/include/Rivet/Makefile.am
+++ b/include/Rivet/Makefile.am
@@ -1,156 +1,157 @@
## Internal headers - not to be installed
nobase_dist_noinst_HEADERS =
## Public headers - to be installed
nobase_pkginclude_HEADERS =
## Rivet interface
nobase_pkginclude_HEADERS += \
Rivet.hh \
Run.hh \
Event.hh \
ParticleBase.hh \
Particle.fhh Particle.hh \
Jet.fhh Jet.hh \
Projection.fhh Projection.hh \
ProjectionApplier.hh \
ProjectionHandler.hh \
Analysis.hh \
AnalysisHandler.hh \
AnalysisInfo.hh \
AnalysisBuilder.hh \
AnalysisLoader.hh
## Build config stuff
nobase_pkginclude_HEADERS += \
Config/RivetCommon.hh \
Config/RivetConfig.hh \
Config/BuildOptions.hh
## Projections
nobase_pkginclude_HEADERS += \
Projections/AxesDefinition.hh \
Projections/Beam.hh \
Projections/BeamThrust.hh \
Projections/CentralEtHCM.hh \
+ Projections/Centrality.hh \
Projections/ChargedFinalState.hh \
Projections/ChargedLeptons.hh \
Projections/ConstLossyFinalState.hh \
Projections/DISFinalState.hh \
Projections/DISKinematics.hh \
Projections/DISLepton.hh \
Projections/DressedLeptons.hh \
Projections/FastJets.hh \
Projections/FinalPartons.hh \
Projections/FinalState.hh \
Projections/FoxWolframMoments.hh \
Projections/FParameter.hh \
Projections/HadronicFinalState.hh \
Projections/HeavyHadrons.hh \
Projections/Hemispheres.hh \
Projections/IdentifiedFinalState.hh \
Projections/InitialQuarks.hh \
Projections/InvMassFinalState.hh \
Projections/JetAlg.hh \
Projections/JetShape.hh \
Projections/LeadingParticlesFinalState.hh \
Projections/LossyFinalState.hh \
Projections/MergedFinalState.hh \
Projections/MissingMomentum.hh \
Projections/NeutralFinalState.hh \
Projections/NonHadronicFinalState.hh \
Projections/NonPromptFinalState.hh \
Projections/ParisiTensor.hh \
Projections/ParticleFinder.hh \
Projections/PartonicTops.hh \
Projections/PrimaryHadrons.hh \
Projections/PromptFinalState.hh \
Projections/SmearedParticles.hh \
Projections/SmearedJets.hh \
Projections/SmearedMET.hh \
Projections/Sphericity.hh \
Projections/Spherocity.hh \
Projections/TauFinder.hh \
Projections/Thrust.hh \
Projections/TriggerCDFRun0Run1.hh \
Projections/TriggerCDFRun2.hh \
Projections/TriggerUA5.hh \
Projections/UnstableFinalState.hh \
Projections/VetoedFinalState.hh \
Projections/VisibleFinalState.hh \
Projections/WFinder.hh \
Projections/ZFinder.hh
## Meta-projection convenience headers
nobase_pkginclude_HEADERS += \
Projections/FinalStates.hh \
Projections/Smearing.hh
## Analysis base class headers
nobase_pkginclude_HEADERS += \
Analyses/MC_ParticleAnalysis.hh \
Analyses/MC_JetAnalysis.hh \
Analyses/MC_JetSplittings.hh
## Tools
nobase_pkginclude_HEADERS += \
Tools/BeamConstraint.hh \
Tools/BinnedHistogram.hh \
Tools/CentralityBinner.hh \
Tools/Cmp.fhh \
Tools/Cmp.hh \
Tools/Cutflow.hh \
Tools/Cuts.fhh \
Tools/Cuts.hh \
Tools/Exceptions.hh \
Tools/JetUtils.hh \
Tools/Logging.hh \
Tools/ParticleBaseUtils.hh \
Tools/ParticleIdUtils.hh \
Tools/ParticleUtils.hh \
Tools/ParticleName.hh \
Tools/PrettyPrint.hh \
Tools/RivetPaths.hh \
Tools/RivetSTL.hh \
Tools/RivetFastJet.hh \
Tools/RivetHepMC.hh \
Tools/RivetYODA.hh \
Tools/RivetMT2.hh \
Tools/SmearingFunctions.hh \
Tools/TypeTraits.hh \
Tools/Utils.hh
nobase_dist_noinst_HEADERS += \
Tools/osdir.hh
## Maths
nobase_pkginclude_HEADERS += \
Math/Matrices.hh \
Math/Vector3.hh \
Math/VectorN.hh \
Math/MatrixN.hh \
Math/MatrixDiag.hh \
Math/MathHeader.hh \
Math/Vectors.hh \
Math/LorentzTrans.hh \
Math/Matrix3.hh \
Math/MathUtils.hh \
Math/Vector4.hh \
Math/Math.hh \
Math/Units.hh \
Math/Constants.hh \
Math/eigen/util.h \
Math/eigen/regressioninternal.h \
Math/eigen/regression.h \
Math/eigen/vector.h \
Math/eigen/ludecompositionbase.h \
Math/eigen/ludecomposition.h \
Math/eigen/linearsolver.h \
Math/eigen/linearsolverbase.h \
Math/eigen/matrix.h \
Math/eigen/vectorbase.h \
Math/eigen/projective.h \
Math/eigen/matrixbase.h
diff --git a/include/Rivet/Projections/Centrality.hh b/include/Rivet/Projections/Centrality.hh
new file mode 100644
--- /dev/null
+++ b/include/Rivet/Projections/Centrality.hh
@@ -0,0 +1,96 @@
+// -*- C++ -*-
+#ifndef RIVET_Centrality_HH
+#define RIVET_Centrality_HH
+
+#include "Rivet/Projection.hh"
+#include "YODA/WriterYODA.h"
+
+namespace Rivet {
+
+
+ class Centrality : public Projection {
+ public:
+
+ enum Method {
+ ImpactParameter,
+ Multiplicity
+ };
+
+ Centrality (const FinalState &fs, const Method method, const size_t numEventsRequired = 10000);
+
+ /// Clone on the heap
+ DEFAULT_RIVET_PROJ_CLONE(Centrality);
+
+ /// Is the centrality valid?
+ bool isValid() const { return (_centrality >= 0.); }
+
+ /// Add calibration file
+ void addCalibrationFile(const string& filename);
+
+ /// Get centrality of the event
+ double getCentrality() const { return _centrality; }
+
+ /// Finalize the centrality calculation; save histograms
+ void finalize() const;
+
+ /// Add cut to the projection
+ void addCut(const Cut& cut);
+
+ protected:
+
+ /// Apply the projection to the event
+ virtual void project(const Event& e);
+
+ /// Compare projections
+ virtual int compare(const Projection& p) const;
+
+ /// Calculate the value of the observable using the selected method
+ float calculateObservable(const Event& e) const;
+
+ /// Check if observable is correct
+ void checkObservable(const float observable) const;
+
+ /// Calculate centrality percentile
+ void calculateCentralityPercentile(const float observable);
+
+ /// Save 1D histogram
+ void saveHistogram(const string& filename, const Histo1DPtr hist) const;
+
+ /// Trim string by removing some special characters
+ void trimString(string& str) const;
+
+ /// Add cut string to the projection
+ void addCutString(string cutString);
+
+ /// When there is no implementation for requested method
+ void method_not_found() const;
+
+ private:
+
+ /// Histogram for centrality calibration
+ Histo1DPtr _histCentralityCalibration;
+ /// Histogram for centrality control. It may be used to compare distribution
+ /// in the current run to the one provided in calibration histogram.
+ Histo1DPtr _histCentralityControl;
+
+ /// String with the cuts
+ string _cutString;
+
+ /// Method of the centrality calibration
+ Method _method;
+
+ /// Number of events required for a selected method of centrality calibration
+ size_t _numEventsRequired;
+
+ /// Centrality value
+ float _centrality;
+
+ /// Name of the centrality calibration method
+ string _methodName;
+
+ };
+
+}
+
+
+#endif
diff --git a/include/Rivet/Projections/FinalState.hh b/include/Rivet/Projections/FinalState.hh
--- a/include/Rivet/Projections/FinalState.hh
+++ b/include/Rivet/Projections/FinalState.hh
@@ -1,55 +1,57 @@
// -*- C++ -*-
#ifndef RIVET_FinalState_HH
#define RIVET_FinalState_HH
#include "Rivet/Projections/ParticleFinder.hh"
namespace Rivet {
/// @brief Project out all final-state particles in an event.
/// Probably the most important projection in Rivet!
class FinalState : public ParticleFinder {
private:
// hide lossy copy constructors for all classes derived from FinalState
template<typename T> FinalState(const T& rhs);
template<typename T> FinalState const& operator=(T const& rhs);
public:
/// @name Standard constructors etc.
//@{
/// Construction using Cuts object
FinalState(const Cut& c=Cuts::open());
// /// Construction using Cuts object and another FinalState
// FinalState(const Cut& c=Cuts::open(), const FinalState& fsp=FinalState());
/// Old constructor with numeric cut arguments, retained for compatibility
/// @deprecated Use the versions with Cut arguments
FinalState(double mineta, double maxeta, double minpt=0.0*GeV);
/// Clone on the heap.
DEFAULT_RIVET_PROJ_CLONE(FinalState);
//@}
/// Apply the projection to the event.
virtual void project(const Event& e);
/// Compare projections.
virtual int compare(const Projection& p) const;
/// Decide if a particle is to be accepted or not.
/// @todo Rename to _accept or acceptFinal?
virtual bool accept(const Particle& p) const;
+ /// Get the cut object
+ virtual Cut getCuts() const { return _cuts; }
};
}
#endif
diff --git a/src/Projections/Centrality.cc b/src/Projections/Centrality.cc
new file mode 100644
--- /dev/null
+++ b/src/Projections/Centrality.cc
@@ -0,0 +1,258 @@
+// -*- C++ -*-
+#include "Rivet/Particle.hh"
+#include "Rivet/Event.hh"
+#include "Rivet/Tools/RivetYODA.hh"
+#include "Rivet/Analysis.hh"
+#include "YODA/WriterYODA.h"
+#include "YODA/ReaderYODA.h"
+#include "Rivet/Projections/FinalState.hh"
+#include "Rivet/Projections/Centrality.hh"
+
+#include <algorithm>
+
+namespace Rivet {
+
+ Centrality::Centrality(const FinalState &fs, const Method method, const size_t numEventsRequired) :
+ _cutString("true"),
+ _method(method),
+ _numEventsRequired(numEventsRequired),
+ _centrality(-1.)
+ {
+ addProjection(fs, "FS");
+
+ switch (_method) {
+
+ case ImpactParameter:
+ _methodName = "ImpactParameter";
+ _histCentralityCalibration = make_shared<Histo1D>(200, 0.0, 20.0);
+ _histCentralityControl = make_shared<Histo1D>(200, 0.0, 20.0);
+ break;
+
+ case Multiplicity:
+ _methodName = "Multiplicity";
+ _histCentralityCalibration = make_shared<Histo1D>(10000, 0.0, 10000.0);
+ _histCentralityControl = make_shared<Histo1D>(10000, 0.0, 10000.0);
+ break;
+
+ default:
+ method_not_found();
+
+ }
+
+ addCut(fs.getCuts());
+ MSG_INFO("Using centrality projection for cuts " + _cutString + ", calibration method " + _methodName + ", and " << _numEventsRequired << " required events");
+ }
+
+
+ int Centrality::compare(const Projection& p) const {
+ // @note Is it okay to cast projetion to centrality? What if it fails?
+ const Centrality& other = dynamic_cast<const Centrality&>(p);
+ const FinalState& fs = getProjection<FinalState>("FS");
+ const FinalState& otherfs = other.getProjection<FinalState>("FS");
+ return \
+ (cmp(fs.getCuts(), otherfs.getCuts()) ||
+ cmp(_method, other._method) ||
+ cmp(_numEventsRequired, other._numEventsRequired) ||
+ mkNamedPCmp(p, "FS"));
+ }
+
+
+ void Centrality::project(const Event& e) {
+
+ // Set the default value of centrality
+ _centrality = -1.;
+
+ // Initialize value of the event observable to be filled to the histograms
+ const float observable = calculateObservable(e);
+
+ // Check if there are enough events in the calibration histogram
+ if (_histCentralityCalibration->numEntries() >= _numEventsRequired) {
+ // Calculate centrality as percentile using observable calculated with the selected method
+ calculateCentralityPercentile(observable);
+ // Fill the control histogram with the impact parameter value and event weight
+ MSG_DEBUG("Adding control point nr " << _histCentralityControl->numEntries() << ": (" << observable << ", " << e.weight() << ")");
+ _histCentralityControl->fill(observable, e.weight());
+ }
+ // Otherwise, if there are not enough events in the calibration histogram
+ else {
+ // Fill the calibration histogram with the impact parameter value and event weight
+ MSG_DEBUG("Adding calibration point nr " << _histCentralityCalibration->numEntries() << ": (" << observable << ", " << e.weight() << ")");
+ _histCentralityCalibration->fill(observable, e.weight());
+ }
+
+ }
+
+
+ float Centrality::calculateObservable(const Event& e) const {
+
+ // Initialize observable
+ float observable = -1.;
+
+ // Calculate its value according to the selected method
+ switch (_method) {
+
+ case ImpactParameter:
+ // Get impact parameter of the event
+ observable = e.genEvent()->heavy_ion() ? e.genEvent()->heavy_ion()->impact_parameter() : -1.;
+ break;
+
+ case Multiplicity:
+ {
+ // Get multiplicity of the event
+ const FinalState& fs = applyProjection<FinalState>(e, "FS");
+ observable = e.genEvent()->heavy_ion() ? fs.particles().size() : -1.;
+ break;
+ }
+
+ default:
+ method_not_found();
+ }
+
+ // Check if observable is correct. If not, skip this event
+ checkObservable(observable);
+
+ return observable;
+ }
+
+
+ void Centrality::checkObservable(const float observable) const {
+
+ try {
+ // Check if the observable is greater than 0
+ if (observable < 0.)
+ throw UserError("Calculated observable is lower than 0!");
+ // Check if the observable is inside histogram ranges
+ if ((observable < _histCentralityCalibration->xMin()) ||
+ (observable > _histCentralityCalibration->xMax()) ||
+ (_histCentralityCalibration->binIndexAt(observable) < 0))
+ throw RangeError("Calculated observable is out of bounds!");
+ } catch (...) {
+ MSG_WARNING("Skipping this event...");
+ }
+
+ }
+
+
+ void Centrality::calculateCentralityPercentile(const float observable) {
+
+ // Default integration ranges
+ size_t bin1 = 0;
+ size_t bin2 = _histCentralityCalibration->numBins() - 1;
+
+ // Select proper integration ranges according to the selected method
+ switch (_method) {
+
+ case ImpactParameter:
+ // Calculate centrality using impact parameter distribution
+ bin2 = _histCentralityCalibration->binIndexAt(observable);
+ break;
+
+ case Multiplicity:
+ // Calculate centrality using multiplicity distribution
+ bin1 = _histCentralityCalibration->binIndexAt(observable);
+ break;
+
+ default:
+ method_not_found();
+ }
+
+ // Calculate centrality as percentile
+ if ((0 <= bin1) && (bin1 <= bin2) && (bin2 < _histCentralityCalibration->numBins())) {
+ _centrality = _histCentralityCalibration->integralRange(bin1, bin2) * 100.0 / _histCentralityCalibration->integral();
+ }
+
+ }
+
+
+ void Centrality::saveHistogram(const string& filename, const Histo1DPtr hist) const {
+ Histo1DPtr histCopy = make_shared<Histo1D>(*hist);
+ histCopy->scaleW(1. / hist->numEntries());
+ try {
+ MSG_INFO("Saving the centrality histogram to " << filename << "...");
+ YODA::WriterYODA::write(filename, *histCopy);
+ } catch (...) {
+ throw UserError("Unexpected error in writing file to: " + filename);
+ }
+ histCopy.reset();
+ }
+
+
+ void Centrality::finalize() const {
+
+ // Save histograms
+ saveHistogram(name() + "_calibration.yoda", _histCentralityCalibration);
+ if (_histCentralityControl->numEntries() > 0)
+ saveHistogram(name() + "_control.yoda", _histCentralityControl);
+
+ // @note Is there anything else that should be done here?
+
+ }
+
+
+ void Centrality::addCalibrationFile(const string& filename) {
+
+ // Open the calibration file and read all objects
+ std::vector<AnalysisObject*> aos;
+ YODA::Reader& reader = YODA::mkReader(filename);
+ std::ifstream istring(filename);
+ try {
+ MSG_INFO("Reading the calibration file: " << filename << "...");
+ reader.read(istring, aos);
+ } catch (const YODA::ReadError& err) {
+ throw UserError("Unexpected error in reading the file " + filename + ": " + err.what());
+ }
+
+ // Create a histogram path template to look for
+ const string histPath = "/Centrality/" + _methodName + "/" + _cutString + "/calibration";
+ MSG_INFO("Searching for: " << histPath);
+
+ // Extract the proper calibration histogram from the file
+ for (const auto ao : aos) {
+ AnalysisObjectPtr aoPtr(ao);
+ Histo1DPtr hist = dynamic_pointer_cast<Histo1D>(aoPtr);
+ if (hist != nullptr) {
+ MSG_DEBUG(hist->path() << " ? " << histPath);
+ if (hist->path().compare(histPath) == 0) {
+ // Check if calibration histogram already exists
+ if (_histCentralityCalibration && (_histCentralityCalibration->numEntries() > 0)) {
+ MSG_WARNING("Calibration histogram already exists! The histogram found in " << filename << " file will be omitted!");
+ return;
+ }
+ MSG_INFO("Histogram found for centrality projection. Setting the new calibration histogram...");
+ _histCentralityCalibration = make_shared<Histo1D>(*hist);
+ }
+ }
+ else {
+ MSG_ERROR("Unsupported type of histogram");
+ }
+ }
+
+ }
+
+ void Centrality::trimString(string& str) const {
+ // @note Which characters should be removed?
+ const char characters[] = "\\ ";
+ for (const auto character : characters) {
+ str.erase (std::remove (str.begin(), str.end(), character), str.end());
+ }
+ }
+
+ void Centrality::addCut(const Cut& cut) {
+ addCutString(cut->description());
+ }
+
+ void Centrality::addCutString(string cutString) {
+ trimString(cutString);
+ _cutString = cutString;
+ _histCentralityCalibration->setPath("/Centrality/" + _methodName + "/" + _cutString + "/calibration");
+ _histCentralityCalibration->setTitle("/Centrality/" + _methodName + "/" + _cutString + "/calibration");
+ _histCentralityControl->setPath("/Centrality/" + _methodName + "/" + _cutString + "/control");
+ _histCentralityControl->setTitle("/Centrality/" + _methodName + "/" + _cutString + "/control");
+ setName("Centrality_" + _methodName + "_" + _cutString);
+ }
+
+ void Centrality::method_not_found() const {
+ throw Exception("Unimplemented method of centrality calibration.");
+ }
+
+}
diff --git a/src/Projections/ChargedFinalState.cc b/src/Projections/ChargedFinalState.cc
--- a/src/Projections/ChargedFinalState.cc
+++ b/src/Projections/ChargedFinalState.cc
@@ -1,48 +1,50 @@
// -*- C++ -*-
#include "Rivet/Projections/ChargedFinalState.hh"
namespace Rivet {
ChargedFinalState::ChargedFinalState(const FinalState& fsp) {
setName("ChargedFinalState");
addProjection(fsp, "FS");
}
- ChargedFinalState::ChargedFinalState(const Cut& c) {
+ ChargedFinalState::ChargedFinalState(const Cut& c) :
+ FinalState(c)
+ {
setName("ChargedFinalState");
addProjection(FinalState(c), "FS");
}
ChargedFinalState::ChargedFinalState(double mineta, double maxeta, double minpt) {
setName("ChargedFinalState");
addProjection(FinalState(mineta, maxeta, minpt), "FS");
}
int ChargedFinalState::compare(const Projection& p) const {
return mkNamedPCmp(p, "FS");
}
}
namespace {
inline bool chargedParticleFilter(const Rivet::Particle& p) {
return Rivet::PID::threeCharge(p.pdgId()) == 0;
}
}
namespace Rivet {
void ChargedFinalState::project(const Event& e) {
const FinalState& fs = applyProjection<FinalState>(e, "FS");
_theParticles.clear();
std::remove_copy_if(fs.particles().begin(), fs.particles().end(),
std::back_inserter(_theParticles), chargedParticleFilter);
MSG_DEBUG("Number of charged final-state particles = " << _theParticles.size());
if (getLog().isActive(Log::TRACE)) {
for (vector<Particle>::iterator p = _theParticles.begin(); p != _theParticles.end(); ++p) {
MSG_TRACE("Selected: " << p->pdgId() << ", charge = " << PID::threeCharge(p->pdgId())/3.0);
}
}
}
}
diff --git a/src/Projections/Makefile.am b/src/Projections/Makefile.am
--- a/src/Projections/Makefile.am
+++ b/src/Projections/Makefile.am
@@ -1,45 +1,46 @@
noinst_LTLIBRARIES = libRivetProjections.la
libRivetProjections_la_SOURCES = \
Beam.cc \
BeamThrust.cc \
ChargedFinalState.cc \
ChargedLeptons.cc \
CentralEtHCM.cc \
+ Centrality.cc \
DISFinalState.cc \
DISKinematics.cc \
DISLepton.cc \
DressedLeptons.cc \
FastJets.cc \
FinalPartons.cc \
FinalState.cc \
FParameter.cc \
HadronicFinalState.cc \
HeavyHadrons.cc \
Hemispheres.cc \
IdentifiedFinalState.cc \
InitialQuarks.cc \
InvMassFinalState.cc \
JetAlg.cc \
JetShape.cc \
LeadingParticlesFinalState.cc \
MergedFinalState.cc \
MissingMomentum.cc \
NeutralFinalState.cc \
NonHadronicFinalState.cc \
NonPromptFinalState.cc \
ParisiTensor.cc \
ParticleFinder.cc \
PrimaryHadrons.cc \
PromptFinalState.cc \
Sphericity.cc \
Spherocity.cc \
TauFinder.cc \
Thrust.cc \
TriggerCDFRun0Run1.cc \
TriggerCDFRun2.cc \
TriggerUA5.cc \
UnstableFinalState.cc \
VetoedFinalState.cc \
VisibleFinalState.cc \
WFinder.cc \
ZFinder.cc

File Metadata

Mime Type
text/x-diff
Expires
Sat, May 3, 6:00 AM (1 d, 3 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4982905
Default Alt Text
(20 KB)

Event Timeline