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 FinalState(const T& rhs); template 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 + +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(200, 0.0, 20.0); + _histCentralityControl = make_shared(200, 0.0, 20.0); + break; + + case Multiplicity: + _methodName = "Multiplicity"; + _histCentralityCalibration = make_shared(10000, 0.0, 10000.0); + _histCentralityControl = make_shared(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(p); + const FinalState& fs = getProjection("FS"); + const FinalState& otherfs = other.getProjection("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(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(*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 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(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(*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(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::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