diff --git a/analyses/pluginALICE/ALICE_2012_I1127497.cc b/analyses/pluginALICE/ALICE_2012_I1127497.cc --- a/analyses/pluginALICE/ALICE_2012_I1127497.cc +++ b/analyses/pluginALICE/ALICE_2012_I1127497.cc @@ -1,180 +1,185 @@ // -*- C++ -*- #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Projections/SingleValueProjection.hh" #include "Rivet/Tools/AliceCommon.hh" #include "Rivet/Projections/AliceCommon.hh" +#include "Rivet/Projections/HepMCHeavyIon.hh" #include #define _USE_MATH_DEFINES #include namespace Rivet { /// Analysis class ALICE_2012_I1127497 : public Analysis { public: DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2012_I1127497); /// Book histograms and initialise projections before the run void init() { // Declare centrality projection if (getOption("cent") == "REF") { _mode = 1; declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M", "V0M"); } else if (getOption("cent") == "IMP") { _mode = 2; declareCentrality(ALICE::V0MMultiplicity(), "ALICE_2015_PBPBCentrality", "V0M", "V0M_IMP"); } // Charged final states with |eta| < 0.5 and pT > 150 MeV const Cut& cut = Cuts::abseta < 0.5 && Cuts::pT > 150*MeV; const ChargedFinalState cfs(cut); addProjection(cfs,"CFS"); + // Access the HepMC heavy ion info + declare(HepMCHeavyIon(), "HepMC"); + // Loop over all histograms for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { // Initialize PbPb objects _histNch[PBPB][ihist] = bookHisto1D(ihist+1, 1, 1); std::string nameCounterPbPb = "Counter_PbPb_" + std::to_string(ihist); _counterSOW[PBPB][ihist] = bookCounter(nameCounterPbPb, "Sum of weights counter for PbPb"); std::string nameCounterNcoll = "Counter_Ncoll_" + std::to_string(ihist); _counterNcoll[ihist] = bookCounter(nameCounterNcoll, "Ncoll counter for PbPb"); // Initialize pp objects. In principle, only one pp histogram would be needed since // centrality does not make any difference here. However, in some cases in this analysis // the binning differ from each other, so this is easy-to-implement way to account for that. std::string namePP = _histNch[PBPB][ihist]->name() + "-pp"; _histNch[PP][ihist] = bookHisto1D(namePP, refData(ihist+1, 1, 1)); // binning taken from ref data std::string nameCounterpp = "Counter_pp_" + std::to_string(ihist); _counterSOW[PP][ihist] = bookCounter(nameCounterpp, "Sum of weights counter for pp"); } // Centrality regions keeping boundaries for a certain region. // Note, that some regions overlap with other regions. _centrRegions.clear(); _centrRegions = {{0., 5.}, {5., 10.}, {10., 20.}, {20., 30.}, {30., 40.}, {40., 50.}, {50., 60.}, {60., 70.}, {70., 80.}, {0., 10.}, {0., 20.}, {20., 40.}, {40., 60.}, {40., 80.}, {60., 80.}}; } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Final state particles with at least pT = 150 MeV in eta range of |eta| < 0.5 const ChargedFinalState& charged = applyProjection(event, "CFS"); Particles chargedParticles = charged.particlesByPt(); // Check type of event. This may be not the perfect way to check for the type of event as // there might be some weird conditions hidden inside. For example some HepMC versions check // if number of hard collisions is equal to 0 and assign 'false' in that case, which is usually wrong. // This might be changed in the future if (event.genEvent()->heavy_ion()) { + const HepMCHeavyIon & hi = apply(event, "HepMC"); // Prepare centrality projection and value const CentralityProjection& centrProj = apply(event, (_mode == 1 ? "V0M" : "V0M_IMP")); double centr = centrProj(); // Veto event for too large centralities since those are not used in the analysis at all if ((centr < 0.) || (centr > 80.)) vetoEvent; // Fill the right PbPb histograms and add weights based on centrality value for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { if (inRange(centr, _centrRegions[ihist].first, _centrRegions[ihist].second)) { _counterSOW[PBPB][ihist]->fill(weight); - _counterNcoll[ihist]->fill(event.genEvent()->heavy_ion()->Ncoll, weight); + _counterNcoll[ihist]->fill(hi.Ncoll(), weight); foreach (const Particle& p, chargedParticles) { float pT = p.pT()/GeV; if (pT < 50.) { double pTAtBinCenter = _histNch[PBPB][ihist]->binAt(pT).xMid(); _histNch[PBPB][ihist]->fill(pT, weight/pTAtBinCenter); } } } } } else { // Fill all pp histograms and add weights for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { _counterSOW[PP][ihist]->fill(weight); foreach (const Particle& p, chargedParticles) { float pT = p.pT()/GeV; if (pT < 50.) { double pTAtBinCenter = _histNch[PP][ihist]->binAt(pT).xMid(); _histNch[PP][ihist]->fill(pT, weight/pTAtBinCenter); } } } } } /// Normalise histograms etc., after the run void finalize() { // Right scaling of the histograms with their individual weights. for (size_t itype = 0; itype < EVENT_TYPES; ++itype ) { for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { if (_counterSOW[itype][ihist]->sumW() > 0.) { scale(_histNch[itype][ihist], (1./_counterSOW[itype][ihist]->sumW() / 2. / M_PI)); } } } // Postprocessing of the histograms for (size_t ihist = 0; ihist < NHISTOS; ++ihist) { // If there are entires in histograms for both beam types if (_histNch[PP][ihist]->numEntries() > 0 && _histNch[PBPB][ihist]->numEntries() > 0) { // Initialize and fill R_AA histograms _histRAA[ihist] = bookScatter2D(ihist+16, 1, 1); divide(_histNch[PBPB][ihist], _histNch[PP][ihist], _histRAA[ihist]); // Scale by Ncoll. Unfortunately some generators does not provide Ncoll (eg. JEWEL), // so the following scaling will be done only if there are entries in the counters if (_counterNcoll[ihist]->sumW() > 1e-6 && _counterSOW[PBPB][ihist]->sumW() > 1e-6) { _histRAA[ihist]->scaleY(1. / (_counterNcoll[ihist]->sumW() / _counterSOW[PBPB][ihist]->sumW())); } } } } private: static const int NHISTOS = 15; static const int EVENT_TYPES = 2; static const int PP = 0; static const int PBPB = 1; Histo1DPtr _histNch[EVENT_TYPES][NHISTOS]; CounterPtr _counterSOW[EVENT_TYPES][NHISTOS]; CounterPtr _counterNcoll[NHISTOS]; Scatter2DPtr _histRAA[NHISTOS]; std::vector> _centrRegions; int _mode; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2012_I1127497); } diff --git a/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc b/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc --- a/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc +++ b/analyses/pluginCMS/CMS_2013_I1224539_ZJET.cc @@ -1,201 +1,200 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/ZFinder.hh" #include "fastjet/tools/Filter.hh" #include "fastjet/tools/Pruner.hh" namespace Rivet { class CMS_2013_I1224539_ZJET : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor CMS_2013_I1224539_ZJET() : Analysis("CMS_2013_I1224539_ZJET"), _filter(fastjet::Filter(fastjet::JetDefinition(fastjet::cambridge_algorithm, 0.3), fastjet::SelectorNHardest(3))), _trimmer(fastjet::Filter(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.03))), _pruner(fastjet::Pruner(fastjet::cambridge_algorithm, 0.1, 0.5)) { } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs(Cuts::abseta < 2.4); declare(fs, "FS"); // Find Zs with pT > 120 GeV ZFinder zfinder(fs, Cuts::abseta < 2.4 && Cuts::pT > 30*GeV, PID::ELECTRON, 80*GeV, 100*GeV, 0.2, ZFinder::CLUSTERNODECAY, ZFinder::TRACK); declare(zfinder, "ZFinder"); // Z+jet jet collections declare(FastJets(zfinder.remainingFinalState(), FastJets::ANTIKT, 0.7), "JetsAK7_zj"); declare(FastJets(zfinder.remainingFinalState(), FastJets::CAM, 0.8), "JetsCA8_zj"); declare(FastJets(zfinder.remainingFinalState(), FastJets::CAM, 1.2), "JetsCA12_zj"); // Histograms /// @note These are 2D histos rendered into slices const int zjetsOffset = 28; for (size_t i = 0; i < N_PT_BINS_vj; ++i ) { _h_ungroomedJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+0*N_PT_BINS_vj, 1, 1); _h_filteredJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+1*N_PT_BINS_vj,1,1); _h_trimmedJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+2*N_PT_BINS_vj,1,1); _h_prunedJetMass_AK7_zj[i] = bookHisto1D(zjetsOffset+i+1+3*N_PT_BINS_vj,1,1); _h_prunedJetMass_CA8_zj[i] = bookHisto1D(zjetsOffset+i+1+4*N_PT_BINS_vj,1,1); if (i > 0) _h_filteredJetMass_CA12_zj[i] = bookHisto1D(zjetsOffset+i+5*N_PT_BINS_vj,1,1); } } bool isBackToBack_zj(const ZFinder& zf, const fastjet::PseudoJet& psjet) { const FourMomentum& z = zf.bosons()[0].momentum(); const FourMomentum& l1 = zf.constituents()[0].momentum(); const FourMomentum& l2 = zf.constituents()[1].momentum(); /// @todo We should make FourMomentum know how to construct itself from a PseudoJet const FourMomentum jmom(psjet.e(), psjet.px(), psjet.py(), psjet.pz()); return (deltaPhi(z, jmom) > 2.0 && deltaR(l1, jmom) > 1.0 && deltaR(l2, jmom) > 1.0); } // Find the pT histogram bin index for value pt (in GeV), to hack a 2D histogram equivalent /// @todo Use a YODA axis/finder alg when available size_t findPtBin(double ptJ) { const double ptBins_vj[N_PT_BINS_vj+1] = { 125.0, 150.0, 220.0, 300.0, 450.0 }; for (size_t ibin = 0; ibin < N_PT_BINS_vj; ++ibin) { if (inRange(ptJ, ptBins_vj[ibin], ptBins_vj[ibin+1])) return ibin; } return N_PT_BINS_vj; } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); // Get the Z const ZFinder& zfinder = apply(event, "ZFinder"); if (zfinder.bosons().size() != 1) vetoEvent; const Particle& z = zfinder.bosons()[0]; - if ( zfinder.constituents().size() < 2 ) { - cerr << "ZFinder problems: #constinuents: " - << zfinder.constituents().size() << endl - << " #Z-constituents: " - << z.constituents().size() << endl; + if (z.constituents().size() < 2) { + MSG_WARNING("Found a Z with less than 2 constituents."); + vetoEvent; } const Particle l1 = zfinder.constituents()[0]; const Particle l2 = zfinder.constituents()[1]; + MSG_DEBUG(l1.pT() << " " << l2.pT()); // Require a high-pT Z (and constituents) if (l1.pT() < 30*GeV || l2.pT() < 30*GeV || z.pT() < 120*GeV) vetoEvent; // AK7 jets const PseudoJets& psjetsAK7_zj = apply(event, "JetsAK7_zj").pseudoJetsByPt(50.0*GeV); if (!psjetsAK7_zj.empty()) { // Get the leading jet and make sure it's back-to-back with the Z const fastjet::PseudoJet& j0 = psjetsAK7_zj[0]; if (isBackToBack_zj(zfinder, j0)) { const size_t njetBin = findPtBin(j0.pt()/GeV); if (njetBin < N_PT_BINS_vj) { fastjet::PseudoJet filtered0 = _filter(j0); fastjet::PseudoJet trimmed0 = _trimmer(j0); fastjet::PseudoJet pruned0 = _pruner(j0); _h_ungroomedJetMass_AK7_zj[njetBin]->fill(j0.m()/GeV, weight); _h_filteredJetMass_AK7_zj[njetBin]->fill(filtered0.m()/GeV, weight); _h_trimmedJetMass_AK7_zj[njetBin]->fill(trimmed0.m()/GeV, weight); _h_prunedJetMass_AK7_zj[njetBin]->fill(pruned0.m()/GeV, weight); } } } // CA8 jets const PseudoJets& psjetsCA8_zj = apply(event, "JetsCA8_zj").pseudoJetsByPt(50.0*GeV); if (!psjetsCA8_zj.empty()) { // Get the leading jet and make sure it's back-to-back with the Z const fastjet::PseudoJet& j0 = psjetsCA8_zj[0]; if (isBackToBack_zj(zfinder, j0)) { const size_t njetBin = findPtBin(j0.pt()/GeV); if (njetBin < N_PT_BINS_vj) { fastjet::PseudoJet pruned0 = _pruner(j0); _h_prunedJetMass_CA8_zj[njetBin]->fill(pruned0.m()/GeV, weight); } } } // CA12 jets const PseudoJets& psjetsCA12_zj = apply(event, "JetsCA12_zj").pseudoJetsByPt(50.0*GeV); if (!psjetsCA12_zj.empty()) { // Get the leading jet and make sure it's back-to-back with the Z const fastjet::PseudoJet& j0 = psjetsCA12_zj[0]; if (isBackToBack_zj(zfinder, j0)) { const size_t njetBin = findPtBin(j0.pt()/GeV); if (njetBin>0 && njetBin < N_PT_BINS_vj) { fastjet::PseudoJet filtered0 = _filter(j0); _h_filteredJetMass_CA12_zj[njetBin]->fill( filtered0.m() / GeV, weight); } } } } /// Normalise histograms etc., after the run void finalize() { const double normalizationVal = 1000; for (size_t i = 0; i < N_PT_BINS_vj; ++i ) { normalize( _h_ungroomedJetMass_AK7_zj[i], normalizationVal); normalize( _h_filteredJetMass_AK7_zj[i], normalizationVal); normalize( _h_trimmedJetMass_AK7_zj[i], normalizationVal); normalize( _h_prunedJetMass_AK7_zj[i], normalizationVal); normalize( _h_prunedJetMass_CA8_zj[i], normalizationVal); if (i > 0) normalize( _h_filteredJetMass_CA12_zj[i], normalizationVal); } } //@} private: /// @name FastJet grooming tools (configured in constructor init list) //@{ const fastjet::Filter _filter; const fastjet::Filter _trimmer; const fastjet::Pruner _pruner; //@} /// @name Histograms //@{ enum BINS_vj { PT_125_150_vj=0, PT_150_220_vj, PT_220_300_vj, PT_300_450_vj, N_PT_BINS_vj }; Histo1DPtr _h_ungroomedJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_filteredJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_trimmedJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_prunedJetMass_AK7_zj[N_PT_BINS_vj]; Histo1DPtr _h_prunedJetMass_CA8_zj[N_PT_BINS_vj]; Histo1DPtr _h_filteredJetMass_CA12_zj[N_PT_BINS_vj]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2013_I1224539_ZJET); } diff --git a/analyses/pluginLHCf/LHCF_2015_I1351909.cc b/analyses/pluginLHCf/LHCF_2015_I1351909.cc --- a/analyses/pluginLHCf/LHCF_2015_I1351909.cc +++ b/analyses/pluginLHCf/LHCF_2015_I1351909.cc @@ -1,311 +1,311 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/Beam.hh" namespace Rivet { /// @brief Add a short analysis description here class LHCF_2015_I1351909 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(LHCF_2015_I1351909); static constexpr bool lhcf_like = true; static constexpr int ndecay = 1; static constexpr int nbeam = 2; static constexpr double D1_begin = 82000.; //mm 60000.; //mm static constexpr double D1_end = 82000; //mm 90000.; //mm static constexpr double IPtoLHCf = 141050.; //mm /// @name Analysis methods bool isParticleFromCollision(Particle p, vector parents, const Beam & beams) { bool beam[nbeam]={false}; if(parents.size()==nbeam) { for ( int ipar=0; ipar < nbeam; ++ipar ) { if ( parents[ipar].genParticle() == beams.beams().first.genParticle() || parents[ipar].genParticle() == beams.beams().second.genParticle() ) beam[ipar] = true; } if(beam[0] && beam[1]) return true; } return false; } bool isParticleFromDecay(Particle p, vector parents) { if(parents.size()==ndecay) return true; else return false; } bool isDeviated(Particle p, Particle parent) { //Select/Remove particles decayed between IP and LHCf ConstGenVertexPtr pv = p.genParticle()->production_vertex(); assert(pv != NULL); const double decay_vertex = pv->position().z()/mm; const double parent_charge = PID::charge(parent.pid()); const double descendant_charge = PID::charge(p.pid()); if(parent_charge == 0) { //Particles produced by neutral parent decay if(descendant_charge == 0) { return false; } else { if(decay_vertex >= D1_end) return false; else return true; //Remove charged descendants produced from decay before end of D1 } } else { //Particles produced by charged parent decay if(decay_vertex <= D1_begin) { if(descendant_charge == 0) return false; else return true; //Remove charged descendants produced from decay before end of D1 } else { return true; //Remove particles produced by charged parent decay after begin of D1 } } return false; } bool isSameParticle(Particle p1, Particle p2) { if(p1.pid() == p2.pid() && mom(p1).t() == mom(p2).t() && mom(p1).x() == mom(p2).x() && mom(p1).y() == mom(p2).y() && mom(p1).z() == mom(p2).z()) return true; else return false; } bool isAlreadyProcessed(Particle p, vector list) { for(unsigned int ipar=0; iparproduction_vertex(); const double x0 = pv->position().x()/mm; const double y0 = pv->position().y()/mm; const double z0 = pv->position().z()/mm; const double px = p.px()/MeV; const double py = p.py()/MeV; const double pz = abs(p.pz()/MeV); const double dist_to_lhcf = IPtoLHCf - z0; const double x1 = x0 + (dist_to_lhcf * px/pz); const double y1 = y0 + (dist_to_lhcf * py/pz); const double r = sqrt(pow(x1, 2.)+pow(y1, 2.)); const double theta = atan(abs(r / IPtoLHCf)); const double pseudorapidity = - log (tan (theta/2.) ); return pseudorapidity; } /// Book histograms and initialise projections before the run void init() { // Initialise and register projections // declare(FinalState("FS"); - addProjection(FinalState(), "Beam"); + addProjection(FinalState(), "FS"); addProjection(Beam(), "Beams"); // Book histograms _h_n_en_eta1 = bookHisto1D(1, 1, 1); _h_n_en_eta2 = bookHisto1D(1, 1, 2); _h_n_en_eta3 = bookHisto1D(1, 1, 3); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = event.weight(); const FinalState &fs = applyProjection (event, "FS"); Particles fs_particles = fs.particles(); const Beam & beams = applyProjection (event, "Beam"); vector processed_parents; processed_parents.clear(); for (Particle& p: fs_particles ) { if(p.pz()/GeV<0.) continue; double eta = 0.; double en = 0.; if(lhcf_like) { //====================================================================== //========== LHCf-like analysis ======================================== //====================================================================== vector parents = p.parents(); if(isParticleFromCollision(p, parents, beams)) { //Particles directly produced in collisions if(!PID::isHadron(p.pid())) continue; //Remove non-hadron particles if(PID::charge(p.pid()) != 0) continue; //Remove charged particles eta = p.eta(); en = p.E()/GeV; } else if(isParticleFromDecay(p, parents)) { //Particles produced from decay ConstGenVertexPtr pv = p.genParticle()->production_vertex(); assert(pv != NULL); const double decay_vertex = pv->position().z()/mm; Particle parent = parents[0]; if(decay_vertex < IPtoLHCf) { //If decay happens before LHCf we consider descendants if(!PID::isHadron(p.pid())) continue; //Remove non-hadron descendants if(isDeviated(p, parent)) continue; //Remove descendants deviated by D1 eta = RecomputeEta(p); en = p.E()/GeV; } else {//If decay happens after LHCf we consider parents vector ancestors; ancestors.clear(); int ngeneration=0; bool isValid=true; bool isEnded=false; while(!isEnded) //Loop over all generations in the decay { vector temp_part; temp_part.clear(); if(ngeneration==0) { parent = parents[0]; temp_part = parent.parents(); } else { parent = ancestors[0]; temp_part = parent.parents(); } ancestors.clear(); ancestors = temp_part; Particle ancestor = ancestors[0]; if(isParticleFromCollision(parent, ancestors, beams)) { //if we found first particles produced in collisions we consider them isEnded=true; if(!PID::isHadron(parent.pid())) isValid=false; //Remove non-hadron ancestors/parents if(PID::charge(parent.pid()) != 0) isValid=false; //Remove charged ancestors/parents if(isAlreadyProcessed(parent, processed_parents)) isValid=false; //Remove already processed ancestors/parents when looping other descendants else processed_parents.push_back(parent); //Fill ancestors/parents in the list eta = parent.eta(); en = parent.E()/GeV; } else if (isParticleFromDecay(parent, ancestors)) { //if we found first particles produced entering LHCf we consider them ConstGenVertexPtr pv_prev = parent.genParticle()->production_vertex(); assert(pv_prev != NULL); const double previous_decay_vertex = pv_prev->position().z()/mm; if(previous_decay_vertex < IPtoLHCf) { isEnded=true; if(!PID::isHadron(parent.pid())) isValid=false; //Remove non-hadron ancestors/parents if(isDeviated(parent, ancestor)) isValid=false; //Remove ancestors/parents deviated by D1 if(isAlreadyProcessed(parent, processed_parents)) isValid=false; //Remove already processed ancestors/parents when looping other descendants else processed_parents.push_back(parent); //Fill ancestors/parents in the list eta = RecomputeEta(parent); en = parent.E()/GeV; } } else { //This condition should never happen cout << "Looping over particles generation ended without match : Exit..." << endl; exit(EXIT_FAILURE); } ++ngeneration; } if(!isValid) continue; } } else { //This condition should never happen cout << "Particle seems not to be produced in collision or decay : Exit..." << endl; exit(EXIT_FAILURE); } } else { //====================================================================== //========== Only neutrons at IP ======================================= //====================================================================== vector parents = p.parents(); //if(isParticleFromCollision(p, parents)) { //Particles directly produced in collisions if(p.pid() != 2112 ) continue; eta = p.eta(); en = p.E()/GeV; //} } // Fill histograms if( eta > 10.76 ){ _h_n_en_eta1->fill( en , weight ); }else if(eta > 8.99 && eta < 9.22){ _h_n_en_eta2->fill( en , weight ); }else if(eta > 8.81 && eta < 8.99){ _h_n_en_eta3->fill( en , weight ); } } } /// Normalise histograms etc., after the run void finalize() { scale(_h_n_en_eta1, crossSection()/millibarn/sumOfWeights()); // norm to cross section scale(_h_n_en_eta2, crossSection()/millibarn/sumOfWeights()); // norm to cross section scale(_h_n_en_eta3, crossSection()/millibarn/sumOfWeights()); // norm to cross section } //@} private: /// @name Histograms //@{ Histo1DPtr _h_n_en_eta1; Histo1DPtr _h_n_en_eta2; Histo1DPtr _h_n_en_eta3; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(LHCF_2015_I1351909); } diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh --- a/include/Rivet/Particle.hh +++ b/include/Rivet/Particle.hh @@ -1,711 +1,711 @@ // -*- C++ -*- #ifndef RIVET_Particle_HH #define RIVET_Particle_HH #include "Rivet/Particle.fhh" #include "Rivet/ParticleBase.hh" #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Tools/Utils.hh" #include "Rivet/Math/LorentzTrans.hh" // NOTE: Rivet/Tools/ParticleUtils.hh included at the end #include "fastjet/PseudoJet.hh" namespace Rivet { /// Particle representation, either from a HepMC::GenEvent or reconstructed. class Particle : public ParticleBase { public: /// @name Constructors //@{ /// Default constructor. /// @note A particle without info is useless. This only exists to keep STL containers happy. Particle() : ParticleBase(), _original(nullptr), _id(PID::ANY), _isDirect{false,false} { } /// Constructor without GenParticle. Particle(PdgId pid, const FourMomentum& mom, const FourVector& pos=FourVector()) : ParticleBase(), _original(nullptr), _id(pid), _momentum(mom), _origin(pos), _isDirect{false,false} { } /// Constructor from a HepMC GenParticle pointer. Particle(ConstGenParticlePtr gp) : ParticleBase(), _original(gp), _id(gp->pdg_id()), _momentum(gp->momentum()), _isDirect{false,false} { ConstGenVertexPtr vprod = gp->production_vertex(); if (vprod != nullptr) { setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z()); } } /// Constructor from a HepMC GenParticle reference. Particle(const RivetHepMC::GenParticle& gp) - : Particle(gp.shared_from_this()) + : Particle(HepMCUtils::getParticlePtr(gp)) { } //@} /// @name Kinematic properties //@{ /// The momentum. const FourMomentum& momentum() const { return _momentum; } /// Set the momentum. Particle& setMomentum(const FourMomentum& momentum) { _momentum = momentum; return *this; } /// Set the momentum via components. Particle& setMomentum(double E, double px, double py, double pz) { _momentum = FourMomentum(E, px, py, pz); return *this; } /// Apply an active Lorentz transform to this particle Particle& transformBy(const LorentzTransform& lt); //@ /// @name Positional properties //@{ /// The origin position. const FourVector& origin() const { return _origin; } /// Set the origin position. Particle& setOrigin(const FourVector& position) { _origin = position; return *this; } /// Set the origin position via components. Particle& setOrigin(double t, double x, double y, double z) { _origin = FourMomentum(t, x, y, z); return *this; } //@} /// @name Other representations and implicit casts to momentum-like objects //@{ /// Converter to FastJet3 PseudoJet virtual fastjet::PseudoJet pseudojet() const { return fastjet::PseudoJet(mom().px(), mom().py(), mom().pz(), mom().E()); } /// Cast operator to FastJet3 PseudoJet operator PseudoJet () const { return pseudojet(); } /// Get a const pointer to the original GenParticle ConstGenParticlePtr genParticle() const { return _original; } /// Cast operator for conversion to GenParticle* operator ConstGenParticlePtr () const { return genParticle(); } //@} /// @name Particle ID code accessors //@{ /// This Particle's PDG ID code. PdgId pid() const { return _id; } /// Absolute value of the PDG ID code. PdgId abspid() const { return std::abs(_id); } /// This Particle's PDG ID code (alias). /// @deprecated Prefer the pid/abspid form PdgId pdgId() const { return _id; } //@} /// @name Charge //@{ /// The charge of this Particle. double charge() const { return PID::charge(pid()); } /// The absolute charge of this Particle. double abscharge() const { return PID::abscharge(pid()); } /// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge). int charge3() const { return PID::charge3(pid()); } /// Alias for charge3 /// @deprecated Use charge3 int threeCharge() const { return PID::threeCharge(pid()); } /// Three times the absolute charge of this Particle (i.e. integer multiple of smallest quark charge). int abscharge3() const { return PID::abscharge3(pid()); } /// Is this Particle charged? bool isCharged() const { return charge3() != 0; } //@} /// @name Particle species //@{ /// Is this a hadron? bool isHadron() const { return PID::isHadron(pid()); } /// Is this a meson? bool isMeson() const { return PID::isMeson(pid()); } /// Is this a baryon? bool isBaryon() const { return PID::isBaryon(pid()); } /// Is this a lepton? bool isLepton() const { return PID::isLepton(pid()); } /// Is this a charged lepton? bool isChargedLepton() const { return PID::isChargedLepton(pid()); } /// Is this a neutrino? bool isNeutrino() const { return PID::isNeutrino(pid()); } /// Does this (hadron) contain a b quark? bool hasBottom() const { return PID::hasBottom(pid()); } /// Does this (hadron) contain a c quark? bool hasCharm() const { return PID::hasCharm(pid()); } // /// Does this (hadron) contain an s quark? // bool hasStrange() const { return PID::hasStrange(pid()); } /// Is this particle potentially visible in a detector? bool isVisible() const; /// Is this a parton? (Hopefully not very often... fiducial FTW) bool isParton() const { return PID::isParton(pid()); } //@} /// @name Constituents (for composite particles) //@{ /// Set direct constituents of this particle virtual void setConstituents(const Particles& cs, bool setmom=false); /// Add a single direct constituent to this particle virtual void addConstituent(const Particle& c, bool addmom=false); /// Add direct constituents to this particle virtual void addConstituents(const Particles& cs, bool addmom=false); /// Determine if this Particle is a composite of other Rivet Particles bool isComposite() const { return !constituents().empty(); } /// @brief Direct constituents of this particle, returned by reference /// /// The returned vector will be empty if this particle is non-composite, /// and its entries may themselves be composites. const Particles& constituents() const { return _constituents; } /// @brief Direct constituents of this particle, sorted by a functor /// @note Returns a copy, thanks to the sorting const Particles constituents(const ParticleSorter& sorter) const { return sortBy(constituents(), sorter); } /// @brief Direct constituents of this particle, filtered by a Cut /// @note Returns a copy, thanks to the filtering const Particles constituents(const Cut& c) const { return filter_select(constituents(), c); } /// @brief Direct constituents of this particle, sorted by a functor /// @note Returns a copy, thanks to the filtering and sorting const Particles constituents(const Cut& c, const ParticleSorter& sorter) const { return sortBy(constituents(c), sorter); } /// @brief Direct constituents of this particle, filtered by a selection functor /// @note Returns a copy, thanks to the filtering const Particles constituents(const ParticleSelector& selector) const { return filter_select(constituents(), selector); } /// @brief Direct constituents of this particle, filtered and sorted by functors /// @note Returns a copy, thanks to the filtering and sorting const Particles constituents(const ParticleSelector& selector, const ParticleSorter& sorter) const { return sortBy(constituents(selector), sorter); } /// @brief Fundamental constituents of this particle /// @note Returns {{*this}} if this particle is non-composite. Particles rawConstituents() const; /// @brief Fundamental constituents of this particle, sorted by a functor /// @note Returns a copy, thanks to the sorting const Particles rawConstituents(const ParticleSorter& sorter) const { return sortBy(rawConstituents(), sorter); } /// @brief Fundamental constituents of this particle, filtered by a Cut /// @note Returns a copy, thanks to the filtering const Particles rawConstituents(const Cut& c) const { return filter_select(rawConstituents(), c); } /// @brief Fundamental constituents of this particle, sorted by a functor /// @note Returns a copy, thanks to the filtering and sorting const Particles rawConstituents(const Cut& c, const ParticleSorter& sorter) const { return sortBy(rawConstituents(c), sorter); } /// @brief Fundamental constituents of this particle, filtered by a selection functor /// @note Returns a copy, thanks to the filtering const Particles rawConstituents(const ParticleSelector& selector) const { return filter_select(rawConstituents(), selector); } /// @brief Fundamental constituents of this particle, filtered and sorted by functors /// @note Returns a copy, thanks to the filtering and sorting const Particles rawConstituents(const ParticleSelector& selector, const ParticleSorter& sorter) const { return sortBy(rawConstituents(selector), sorter); } //@} /// @name Ancestry (for fundamental particles with a HepMC link) //@{ /// Get a list of the direct parents of the current particle (with optional selection Cut) /// /// @note This is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! Particles parents(const Cut& c=Cuts::OPEN) const; /// Get a list of the direct parents of the current particle (with selector function) /// /// @note This is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! Particles parents(const ParticleSelector& f) const { return filter_select(parents(), f); } /// Check whether any particle in the particle's parent list has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasParentWith(const ParticleSelector& f) const { return !parents(f).empty(); } /// Check whether any particle in the particle's parent list has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasParentWith(const Cut& c) const; /// Check whether any particle in the particle's parent list does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasParentWithout(const ParticleSelector& f) const { return hasParentWith([&](const Particle& p){ return !f(p); }); } /// Check whether any particle in the particle's parent list does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasParentWithout(const Cut& c) const; /// Check whether a given PID is found in the particle's parent list /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! /// /// @deprecated Prefer e.g. hasParentWith(Cut::pid == 123) //DEPRECATED("Prefer e.g. hasParentWith(Cut::pid == 123)"); bool hasParent(PdgId pid) const; /// Get a list of the ancestors of the current particle (with optional selection Cut) /// /// @note By default only physical ancestors, with status=2, are returned. /// /// @note This is valid in MC, but may not be answerable experimentally -- /// use this function with care when replicating experimental analyses! Particles ancestors(const Cut& c=Cuts::OPEN, bool only_physical=true) const; /// Get a list of the direct parents of the current particle (with selector function) /// /// @note By default only physical ancestors, with status=2, are returned. /// /// @note This is valid in MC, but may not be answerable experimentally -- /// use this function with care when replicating experimental analyses! Particles ancestors(const ParticleSelector& f, bool only_physical=true) const { return filter_select(ancestors(Cuts::OPEN, only_physical), f); } /// Check whether any particle in the particle's ancestor list has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasAncestorWith(const ParticleSelector& f, bool only_physical=true) const { return !ancestors(f, only_physical).empty(); } /// Check whether any particle in the particle's ancestor list has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasAncestorWith(const Cut& c, bool only_physical=true) const; /// Check whether any particle in the particle's ancestor list does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasAncestorWithout(const ParticleSelector& f, bool only_physical=true) const { return hasAncestorWith([&](const Particle& p){ return !f(p); }, only_physical); } /// Check whether any particle in the particle's ancestor list does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasAncestorWithout(const Cut& c, bool only_physical=true) const; /// Check whether a given PID is found in the particle's ancestor list /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! /// /// @deprecated Prefer hasAncestorWith(Cuts::pid == pid) etc. //DEPRECATED("Prefer e.g. hasAncestorWith(Cut::pid == 123)"); bool hasAncestor(PdgId pid, bool only_physical=true) const; /// @brief Determine whether the particle is from a b-hadron decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromBottom() const; /// @brief Determine whether the particle is from a c-hadron decay /// /// @note If a hadron contains b and c quarks it is considered a bottom /// hadron and NOT a charm hadron. /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromCharm() const; // /// @brief Determine whether the particle is from a s-hadron decay // /// // /// @note If a hadron contains b or c quarks as well as strange it is // /// considered a b or c hadron, but NOT a strange hadron. // /// // /// @note This question is valid in MC, but may not be perfectly answerable // /// experimentally -- use this function with care when replicating // /// experimental analyses! // bool fromStrange() const; /// @brief Determine whether the particle is from a hadron decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromHadron() const; /// @brief Determine whether the particle is from a tau decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromTau(bool prompt_taus_only=false) const; /// @brief Determine whether the particle is from a prompt tau decay /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromPromptTau() const { return fromTau(true); } /// @brief Determine whether the particle is from a tau which decayed hadronically /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool fromHadronicTau(bool prompt_taus_only=false) const; /// @brief Determine whether the particle is from a hadron or tau decay /// /// Specifically, walk up the ancestor chain until a status 2 hadron or /// tau is found, if at all. /// /// @note This question is valid in MC, but may not be perfectly answerable /// experimentally -- use this function with care when replicating /// experimental analyses! /// /// @deprecated Too vague: use fromHadron or fromHadronicTau bool fromDecay() const { return fromHadron() || fromPromptTau(); } /// @brief Shorthand definition of 'promptness' based on set definition flags /// /// A "direct" particle is one directly connected to the hard process. It is a /// preferred alias for "prompt", since it has no confusing implications about /// distinguishability by timing information. /// /// The boolean arguments allow a decay lepton to be considered direct if /// its parent was a "real" direct lepton. /// /// @note This one doesn't make any judgements about final-stateness bool isDirect(bool allow_from_direct_tau=false, bool allow_from_direct_mu=false) const; /// Alias for isDirect bool isPrompt(bool allow_from_prompt_tau=false, bool allow_from_prompt_mu=false) const { return isDirect(allow_from_prompt_tau, allow_from_prompt_mu); } //@} /// @name Decay info //@{ /// Whether this particle is stable according to the generator bool isStable() const; /// @todo isDecayed? How to restrict to physical particles? /// Get a list of the direct descendants from the current particle (with optional selection Cut) Particles children(const Cut& c=Cuts::OPEN) const; /// Get a list of the direct descendants from the current particle (with selector function) Particles children(const ParticleSelector& f) const { return filter_select(children(), f); } /// Check whether any direct child of this particle has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasChildWith(const ParticleSelector& f) const { return !children(f).empty(); } /// Check whether any direct child of this particle has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasChildWith(const Cut& c) const; /// Check whether any direct child of this particle does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasChildWithout(const ParticleSelector& f) const { return hasChildWith([&](const Particle& p){ return !f(p); }); } /// Check whether any direct child of this particle does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasChildWithout(const Cut& c) const; /// Get a list of all the descendants from the current particle (with optional selection Cut) Particles allDescendants(const Cut& c=Cuts::OPEN, bool remove_duplicates=true) const; /// Get a list of all the descendants from the current particle (with selector function) Particles allDescendants(const ParticleSelector& f, bool remove_duplicates=true) const { return filter_select(allDescendants(Cuts::OPEN, remove_duplicates), f); } /// Check whether any descendant of this particle has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasDescendantWith(const ParticleSelector& f, bool remove_duplicates=true) const { return !allDescendants(f, remove_duplicates).empty(); } /// Check whether any descendant of this particle has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasDescendantWith(const Cut& c, bool remove_duplicates=true) const; /// Check whether any descendant of this particle does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasDescendantWithout(const ParticleSelector& f, bool remove_duplicates=true) const { return hasDescendantWith([&](const Particle& p){ return !f(p); }, remove_duplicates); } /// Check whether any descendant of this particle does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasDescendantWithout(const Cut& c, bool remove_duplicates=true) const; /// Get a list of all the stable descendants from the current particle (with optional selection Cut) /// /// @todo Use recursion through replica-avoiding MCUtils functions to avoid bookkeeping duplicates /// @todo Insist that the current particle is post-hadronization, otherwise throw an exception? Particles stableDescendants(const Cut& c=Cuts::OPEN) const; /// Get a list of all the stable descendants from the current particle (with selector function) Particles stableDescendants(const ParticleSelector& f) const { return filter_select(stableDescendants(), f); } /// Check whether any stable descendant of this particle has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasStableDescendantWith(const ParticleSelector& f) const { return !stableDescendants(f).empty(); } /// Check whether any stable descendant of this particle has the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasStableDescendantWith(const Cut& c) const; /// Check whether any stable descendant of this particle does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasStableDescendantWithout(const ParticleSelector& f) const { return hasStableDescendantWith([&](const Particle& p){ return !f(p); }); } /// Check whether any stable descendant of this particle does not have the requested property /// /// @note This question is valid in MC, but may not be answerable /// experimentally -- use this function with care when replicating /// experimental analyses! bool hasStableDescendantWithout(const Cut& c) const; /// Flight length (divide by mm or cm to get the appropriate units) double flightLength() const; //@} /// @name Duplicate testing //@{ /// @brief Determine whether a particle is the first in a decay chain to meet the function requirement inline bool isFirstWith(const ParticleSelector& f) const { if (!f(*this)) return false; //< This doesn't even meet f, let alone being the last to do so if (any(parents(), f)) return false; //< If a direct parent has this property, this isn't the first return true; } /// @brief Determine whether a particle is the first in a decay chain not to meet the function requirement inline bool isFirstWithout(const ParticleSelector& f) const { return isFirstWith([&](const Particle& p){ return !f(p); }); } /// @brief Determine whether a particle is the last in a decay chain to meet the function requirement inline bool isLastWith(const ParticleSelector& f) const { if (!f(*this)) return false; //< This doesn't even meet f, let alone being the last to do so if (any(children(), f)) return false; //< If a child has this property, this isn't the last return true; } /// @brief Determine whether a particle is the last in a decay chain not to meet the function requirement inline bool isLastWithout(const ParticleSelector& f) const { return isLastWith([&](const Particle& p){ return !f(p); }); } //@} protected: /// A pointer to the original GenParticle from which this Particle is projected. ConstGenParticlePtr _original; /// Constituent particles if this is a composite (may be empty) Particles _constituents; /// The PDG ID code for this Particle. PdgId _id; /// The momentum of this particle. FourMomentum _momentum; /// The creation position of this particle. FourVector _origin; /// Cached computation of directness, via ancestry. Second element is cache status /// @todo Replace this awkward caching with C++17 std::optional mutable std::pair _isDirect; }; /// @name String representation and streaming support //@{ /// Allow a Particle to be passed to an ostream. std::ostream& operator << (std::ostream& os, const Particle& p); /// Allow ParticlePair to be passed to an ostream. std::ostream& operator << (std::ostream& os, const ParticlePair& pp); //@} } #include "Rivet/Tools/ParticleUtils.hh" #endif diff --git a/include/Rivet/Projections/HepMCHeavyIon.hh b/include/Rivet/Projections/HepMCHeavyIon.hh --- a/include/Rivet/Projections/HepMCHeavyIon.hh +++ b/include/Rivet/Projections/HepMCHeavyIon.hh @@ -1,177 +1,180 @@ // -*- C++ -*- #ifndef RIVET_HepMCHeavyIon_HH #define RIVET_HepMCHeavyIon_HH #include "Rivet/Projection.hh" #include "Rivet/Tools/RivetHepMC.hh" #include "Rivet/Event.hh" namespace Rivet { class HepMCHeavyIon : public Projection { public: /// @name Constructors etc. //@{ /// Constructor HepMCHeavyIon(); /// Clone on the heap. DEFAULT_RIVET_PROJ_CLONE(HepMCHeavyIon); //@} protected: /// Perform the projection on the Event void project(const Event& e); /// Compare with other projections //int compare(const Projection& p) const; // Taken from Thrust.hh int compare(const Projection& p) const { return 0; } public: + /// Check that there were at all any heavy ion info in HepMC + bool ok() const { return _hi != nullptr; } + /// @brief the number of hard nucleon-nucleon collisions. /// /// Model-dependent. Usually the number of nucleon-nucleon /// collisions containing a special signal process. A negative /// value means that the information is not available. int Ncoll_hard() const; /// @brief the number of participating nucleons in the projectile. /// /// The number of nucleons in the projectile participating in an /// inelastic collision (see Ncoll). A negative value means that /// the information is not available. int Npart_proj() const; /// @brief the number of participating nucleons in the target. /// /// The number of nucleons in the target participating in an /// inelastic collision (see Ncoll). A negative value means that /// the information is not available. int Npart_targ() const; /// @brief the number of inelastic nucleon-nucleon collisions. /// /// Note that a one participating nucleon can be involved in many /// inelastic collisions, and that inelastic also includes /// diffractive excitation. A negative value means that the /// information is not available. /// int Ncoll() const; /// @brief Collisions with a diffractively excited target nucleon. /// /// The number of single diffractive nucleon-nucleon collisions /// where the target nucleon is excited. A negative value means /// that the information is not available. int N_Nwounded_collisions() const; /// @brief Collisions with a diffractively excited projectile nucleon. /// /// The number of single diffractive nucleon-nucleon collisions /// where the projectile nucleon is excited. A negative value /// means that the information is not available. int Nwounded_N_collisions() const; /// @brief Non-diffractive or doubly diffractive collisions. /// /// The number of nucleon-nucleon collisions where both projectile /// and target nucleons are wounded. A negative value means that /// the information is not available. int Nwounded_Nwounded_collisions() const; /// @brief The impact parameter. /// /// The impact parameter given in units of femtometer. A negative /// value means that the information is not available. double impact_parameter() const; /// @brief The event plane angle. /// /// The angle wrt. the x-axix of the impact parameter vector /// (pointing frm the target to the projectile). A positive number /// between 0 and two pi. A negative value means that the /// information is not available. double event_plane_angle() const; /// @brief The assumed inelastic nucleon-nucleon cross section /// /// in units of millibarn. As used in a Glauber calculation to /// simulate the distribution in Ncoll. A negative value means /// that the information is not available. double sigma_inel_NN() const; /// @brief The centrality. /// /// The generated centrality in percentiles, where 0 is the /// maximally central and 100 is the minimally central. A negative /// value means that the information is not available. double centrality() const; /// @brief A user defined centrality estimator. /// /// This variable may contain anything a generator feels is /// reasonable for estimating centrality. The value should be /// non-negative, and a low value corresponds to a low /// centrality. A negative value indicatess that the information /// is not available. double user_cent_estimate() const; /// @brief The number of spectator neutrons in the projectile /// /// ie. those that thave not participated in any inelastic /// nucleon-nucleon collision. A negative value indicatess that /// the information is not available. int Nspec_proj_n() const; /// @brief The number of spectator neutrons in the target /// /// ie. those that thave not participated in any inelastic /// nucleon-nucleon collision. A negative value indicatess that /// the information is not available. int Nspec_targ_n() const; /// @brief The number of spectator protons in the projectile /// /// ie. those that thave not participated in any inelastic /// nucleon-nucleon collision. A negative value indicatess that /// the information is not available. int Nspec_proj_p() const; /// @brief The number of spectator protons in the target /// /// ie. those that thave not participated in any inelastic /// nucleon-nucleon collision. A negative value indicatess that /// the information is not available. int Nspec_targ_p() const; /// @brief Participant plane angles /// /// calculated to different orders. The key of the map specifies /// the order, and the value gives to the angle wrt. the /// event plane. map participant_plane_angles() const; /// @brief Eccentricities /// /// Calculated to different orders. The key of the map specifies /// the order, and the value gives the corresponding eccentricity. map eccentricities() const; private: /// A pointer to the actual heavy ion object - RivetHepMC::ConstGenHeavyIonPtr hi; + RivetHepMC::ConstGenHeavyIonPtr _hi; }; } #endif diff --git a/include/Rivet/Tools/RivetHepMC.hh b/include/Rivet/Tools/RivetHepMC.hh --- a/include/Rivet/Tools/RivetHepMC.hh +++ b/include/Rivet/Tools/RivetHepMC.hh @@ -1,96 +1,97 @@ // -*- C++ -*- #ifndef RIVET_RivetHepMC_HH #define RIVET_RivetHepMC_HH #ifdef ENABLE_HEPMC_3 #include "HepMC3/HepMC3.h" #include "HepMC3/Relatives.h" #include "HepMC3/Reader.h" namespace Rivet{ namespace RivetHepMC = HepMC3; using RivetHepMC::ConstGenParticlePtr; using RivetHepMC::ConstGenVertexPtr; using RivetHepMC::Relatives; using RivetHepMC::ConstGenHeavyIonPtr; using HepMC_IO_type = RivetHepMC::Reader; using PdfInfo = RivetHepMC::GenPdfInfo; } #else #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/GenVertex.h" #include "HepMC/Version.h" #include "HepMC/GenRanges.h" #include "HepMC/IO_GenEvent.h" namespace Rivet{ namespace RivetHepMC = HepMC; // HepMC 2.07 provides its own #defines #define ConstGenParticlePtr const HepMC::GenParticle* #define ConstGenVertexPtr const HepMC::GenVertex* #define ConstGenHeavyIonPtr const HepMC::HeavyIon* /// @brief Replicated the HepMC3 Relatives syntax using HepMC2 IteratorRanges /// This is necessary mainly because of capitalisation differences class Relatives{ public: constexpr Relatives(HepMC::IteratorRange relo): _internal(relo){} constexpr HepMC::IteratorRange operator()() const {return _internal;} operator HepMC::IteratorRange() const {return _internal;} const static Relatives PARENTS; const static Relatives CHILDREN; const static Relatives ANCESTORS; const static Relatives DESCENDANTS; private: const HepMC::IteratorRange _internal; }; using HepMC_IO_type = HepMC::IO_GenEvent; using PdfInfo = RivetHepMC::PdfInfo; } #endif #include "Rivet/Tools/RivetSTL.hh" #include "Rivet/Tools/Exceptions.hh" namespace Rivet { using RivetHepMC::GenEvent; using ConstGenEventPtr = std::shared_ptr; /// @todo Use mcutils? namespace HepMCUtils{ - + + ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp); std::vector particles(ConstGenEventPtr ge); std::vector particles(const GenEvent *ge); std::vector vertices(ConstGenEventPtr ge); std::vector vertices(const GenEvent *ge); std::vector particles(ConstGenVertexPtr gv, const Relatives &relo); std::vector particles(ConstGenParticlePtr gp, const Relatives &relo); int uniqueId(ConstGenParticlePtr gp); int particles_size(ConstGenEventPtr ge); int particles_size(const GenEvent *ge); std::pair beams(const GenEvent *ge); std::shared_ptr makeReader(std::istream &istr, std::string * errm = 0); bool readEvent(std::shared_ptr io, std::shared_ptr evt); } } #endif diff --git a/src/Projections/HepMCHeavyIon.cc b/src/Projections/HepMCHeavyIon.cc --- a/src/Projections/HepMCHeavyIon.cc +++ b/src/Projections/HepMCHeavyIon.cc @@ -1,104 +1,104 @@ // -*- C++ -*- #include "Rivet/Projections/HepMCHeavyIon.hh" #ifdef ENABLE_HEPMC_3 #define IMPLEMENTATION(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ - return hi? hi->functionname: defret; \ + return _hi? _hi->functionname: defret; \ } #define IMPLEMENTATION_ALT_HEPMC2(rettype, functionname, hepmc2name, defret) \ rettype HepMCHeavyIon::functionname() const { \ - return hi? hi->functionname: defret; \ + return _hi? _hi->functionname: defret; \ } #define IMPLEMENTATION_NO_HEPMC2(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ - return hi? hi->functionname: defret; \ + return _hi? _hi->functionname: defret; \ } #else #define IMPLEMENTATION(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ - return hi? hi->functionname(): defret; \ + return _hi? _hi->functionname(): defret; \ } #define IMPLEMENTATION_ALT_HEPMC2(rettype, functionname, hepmc2name, defret) \ rettype HepMCHeavyIon::functionname() const { \ - return hi? hi->, hepmc2name(): defret; \ + return _hi? _hi->, hepmc2name(): defret; \ } #define IMPLEMENTATION_NO_HEPMC2(rettype, functionname, defret) \ rettype HepMCHeavyIon::functionname() const { \ MSG_WARNING("HeavyIon::" #functionname " is only avialable in HepMC3"); \ return defret; \ } #endif namespace Rivet { HepMCHeavyIon::HepMCHeavyIon() { setName("HepMCHeavyIon"); } void HepMCHeavyIon::project(const Event& e) { - hi = e.genEvent()->heavy_ion(); - if ( !hi ) + _hi = e.genEvent()->heavy_ion(); + if ( !_hi ) MSG_WARNING("Could not find the HepMC HeavyIon object"); } IMPLEMENTATION(int, Ncoll_hard, -1) IMPLEMENTATION(int, Npart_proj, -1) IMPLEMENTATION(int, Npart_targ, -1) IMPLEMENTATION(int, Ncoll, -1) IMPLEMENTATION(int, N_Nwounded_collisions, -1) IMPLEMENTATION(int, Nwounded_N_collisions, -1) IMPLEMENTATION(int, Nwounded_Nwounded_collisions, -1) IMPLEMENTATION(double, impact_parameter, -1.0) IMPLEMENTATION(double, event_plane_angle, -1.0) IMPLEMENTATION(double, sigma_inel_NN, -1.0) IMPLEMENTATION_NO_HEPMC2(double, centrality, -1.0) IMPLEMENTATION_NO_HEPMC2(double, user_cent_estimate, -1.0) IMPLEMENTATION_NO_HEPMC2(int, Nspec_proj_n, -1) IMPLEMENTATION_NO_HEPMC2(int, Nspec_targ_n, -1) IMPLEMENTATION_NO_HEPMC2(int, Nspec_proj_p, -1) IMPLEMENTATION_NO_HEPMC2(int, Nspec_targ_p, -1) map HepMCHeavyIon::participant_plane_angles() const { #ifdef ENABLE_HEPMC_3 - return hi? hi->participant_plane_angles: map(); + return _hi? _hi->participant_plane_angles: map(); #else MSG_WARNING("HeavyIon::participant_plane_angles is only avialable in HepMC3"); return map(); #endif } map HepMCHeavyIon::eccentricities() const { #ifdef ENABLE_HEPMC_3 - return hi? hi->eccentricities: map(); + return _hi? _hi->eccentricities: map(); #else MSG_WARNING("HeavyIon::eccentricities is only avialable in HepMC3"); return map(); #endif } } diff --git a/src/Tools/RivetHepMC_2.cc b/src/Tools/RivetHepMC_2.cc --- a/src/Tools/RivetHepMC_2.cc +++ b/src/Tools/RivetHepMC_2.cc @@ -1,116 +1,119 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" namespace Rivet{ const Relatives Relatives::PARENTS = HepMC::parents; const Relatives Relatives::CHILDREN = HepMC::children; const Relatives Relatives::ANCESTORS = HepMC::ancestors; const Relatives Relatives::DESCENDANTS = HepMC::descendants; namespace HepMCUtils{ + ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) { + return &gp; + } std::vector particles(ConstGenEventPtr ge){ std::vector result; for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){ result.push_back(*pi); } return result; } std::vector particles(const GenEvent *ge){ std::vector result; for(GenEvent::particle_const_iterator pi = ge->particles_begin(); pi != ge->particles_end(); ++pi){ result.push_back(*pi); } return result; } std::vector vertices(ConstGenEventPtr ge){ std::vector result; for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){ result.push_back(*vi); } return result; } std::vector vertices(const GenEvent *ge){ std::vector result; for(GenEvent::vertex_const_iterator vi = ge->vertices_begin(); vi != ge->vertices_end(); ++vi){ result.push_back(*vi); } return result; } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ std::vector result; /// @todo A particle_const_iterator on GenVertex would be nice... // Before HepMC 2.7.0 there were no GV::particles_const_iterators and constness consistency was all screwed up :-/ #if HEPMC_VERSION_CODE >= 2007000 for (HepMC::GenVertex::particle_iterator pi = gv->particles_begin(relo); pi != gv->particles_end(relo); ++pi) result.push_back(*pi); #else HepMC::GenVertex* gv2 = const_cast(gv); for (HepMC::GenVertex::particle_iterator pi = gv2->particles_begin(relo); pi != gv2->particles_end(relo); ++pi) result.push_back(const_cast(*pi)); #endif return result; } std::vector particles(ConstGenParticlePtr gp, const Relatives &relo){ ConstGenVertexPtr vtx; switch(relo){ case HepMC::parents: case HepMC::ancestors: vtx = gp->production_vertex(); break; case HepMC::children: case HepMC::descendants: vtx = gp->end_vertex(); break; default: throw std::runtime_error("Not implemented!"); break; } return particles(vtx, relo); } int uniqueId(ConstGenParticlePtr gp){ return gp->barcode(); } int particles_size(ConstGenEventPtr ge){ return ge->particles_size(); } int particles_size(const GenEvent *ge){ return ge->particles_size(); } std::pair beams(const GenEvent *ge){ return ge->beam_particles(); } std::shared_ptr makeReader(std::istream &istr, std::string *) { return make_shared(istr); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ if(io->rdstate() != 0) return false; if(!io->fill_next_event(evt.get())) return false; return true; } } } diff --git a/src/Tools/RivetHepMC_3.cc b/src/Tools/RivetHepMC_3.cc --- a/src/Tools/RivetHepMC_3.cc +++ b/src/Tools/RivetHepMC_3.cc @@ -1,109 +1,113 @@ // -*- C++ -*- #include "Rivet/Tools/RivetHepMC.hh" #include "Rivet/Tools/ReaderCompressedAscii.hh" #include "HepMC3/ReaderAscii.h" #include "HepMC3/ReaderAsciiHepMC2.h" namespace Rivet{ namespace HepMCUtils{ - + + ConstGenParticlePtr getParticlePtr(const RivetHepMC::GenParticle & gp) { + return gp.shared_from_this(); + } + std::vector particles(ConstGenEventPtr ge){ return ge->particles(); } std::vector particles(const GenEvent *ge){ assert(ge); return ge->particles(); } std::vector vertices(ConstGenEventPtr ge){ return ge->vertices(); } std::vector vertices(const GenEvent *ge){ assert(ge); return ge->vertices(); } std::vector particles(ConstGenVertexPtr gv, const Relatives &relo){ return relo(gv); } std::vector particles(ConstGenParticlePtr gp, const Relatives &relo){ return relo(gp); } int particles_size(ConstGenEventPtr ge){ return particles(ge).size(); } int particles_size(const GenEvent *ge){ return particles(ge).size(); } int uniqueId(ConstGenParticlePtr gp){ return gp->id(); } std::pair beams(const GenEvent *ge) { std::vector beamlist = ge->beams(); if ( beamlist.size() < 2 ) { std::cerr << "CANNOT FIND ANY BEAMS!" << std::endl; return std::pair(); } return std::make_pair(beamlist[0], beamlist[1]); } bool readEvent(std::shared_ptr io, std::shared_ptr evt){ return io->read_event(*evt) && !io->failed(); } shared_ptr makeReader(std::istream & istr, std::string * errm) { shared_ptr ret; // First scan forward and check if there is some hint as to what // kind of file we are looking att. int ntry = 10; std::string header; int filetype = -1; while ( ntry ) { std::getline(istr, header); if ( header.empty() ) continue; if ( header.substr(0, 34) == "HepMC::Asciiv3-START_EVENT_LISTING" ) { filetype = 3; break; } if ( header.substr(0, 44) == "HepMC::CompressedAsciiv3-START_EVENT_LISTING" ) { filetype = 4; break; } if ( header.substr(0, 38) == "HepMC::IO_GenEvent-START_EVENT_LISTING" ) { filetype = 2; break; } --ntry; } if ( filetype == 3 ) ret = make_shared(istr); else if ( filetype == 4 ) ret = make_shared(istr); else ret = make_shared(istr); if ( filetype == 0 && errm ) *errm += "Could not determine file type. Assuming HepMC2 file. "; // Check that everything was ok. if ( ret->failed() ) { if ( errm ) *errm = "Problems reading from HepMC file."; ret = shared_ptr(); } return ret; } } }