diff --git a/analyses/pluginCMS/CMS_2015_I1397174.cc b/analyses/pluginCMS/CMS_2015_I1397174.cc --- a/analyses/pluginCMS/CMS_2015_I1397174.cc +++ b/analyses/pluginCMS/CMS_2015_I1397174.cc @@ -1,386 +1,386 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/PartonicTops.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// Fully leptonic partonic ttbar analysis class CMS_2015_I1397174 : public Analysis { public: /// Minimal constructor CMS_2015_I1397174() : Analysis("CMS_2015_I1397174") { } /// @name Analysis methods //@{ /// Set up projections and book histograms void init() { // Parton level top quarks addProjection(PartonicTops(PartonicTops::E_MU, false), "PartonTops"); // Find jets not related to the top/W decays VetoedFinalState vfs; vfs.addDecayProductsVeto(PID::WPLUSBOSON); vfs.addDecayProductsVeto(PID::WMINUSBOSON); FastJets fj(vfs, FastJets::ANTIKT, 0.5, JetAlg::ALL_MUONS, JetAlg::ALL_INVISIBLES); addProjection(fj, "Jets"); // Book histograms _hVis_nJet30_abs = bookHisto1D( 1, 1, 1); _hVis_nJet30 = bookHisto1D( 2, 1, 1); _hVis_nJet60_abs = bookHisto1D( 3, 1, 1); _hVis_nJet60 = bookHisto1D( 4, 1, 1); _hVis_nJet100_abs = bookHisto1D( 5, 1, 1); _hVis_nJet100 = bookHisto1D( 6, 1, 1); _hVis_addJet1Pt_abs = bookHisto1D( 7, 1, 1); _hVis_addJet1Pt = bookHisto1D( 8, 1, 1); _hVis_addJet1Eta_abs = bookHisto1D( 9, 1, 1); _hVis_addJet1Eta = bookHisto1D(10, 1, 1); _hVis_addJet2Pt_abs = bookHisto1D(11, 1, 1); _hVis_addJet2Pt = bookHisto1D(12, 1, 1); _hVis_addJet2Eta_abs = bookHisto1D(13, 1, 1); _hVis_addJet2Eta = bookHisto1D(14, 1, 1); _hVis_addJJMass_abs = bookHisto1D(15, 1, 1); _hVis_addJJMass = bookHisto1D(16, 1, 1); _hVis_addJJDR_abs = bookHisto1D(17, 1, 1); _hVis_addJJDR = bookHisto1D(18, 1, 1); _hVis_addJJHT_abs = bookHisto1D(19, 1, 1); _hVis_addJJHT = bookHisto1D(20, 1, 1); _hFull_addJet1Pt_abs = bookHisto1D(21, 1, 1); _hFull_addJet1Pt = bookHisto1D(22, 1, 1); _hFull_addJet1Eta_abs = bookHisto1D(23, 1, 1); _hFull_addJet1Eta = bookHisto1D(24, 1, 1); _hFull_addJet2Pt_abs = bookHisto1D(25, 1, 1); _hFull_addJet2Pt = bookHisto1D(26, 1, 1); _hFull_addJet2Eta_abs = bookHisto1D(27, 1, 1); _hFull_addJet2Eta = bookHisto1D(28, 1, 1); _hFull_addJJMass_abs = bookHisto1D(29, 1, 1); _hFull_addJJMass = bookHisto1D(30, 1, 1); _hFull_addJJDR_abs = bookHisto1D(31, 1, 1); _hFull_addJJDR = bookHisto1D(32, 1, 1); _hFull_addJJHT_abs = bookHisto1D(33, 1, 1); _hFull_addJJHT = bookHisto1D(34, 1, 1); _hVis_addBJet1Pt_abs = bookHisto1D(35, 1, 1); _hVis_addBJet1Pt = bookHisto1D(36, 1, 1); _hVis_addBJet1Eta_abs = bookHisto1D(37, 1, 1); _hVis_addBJet1Eta = bookHisto1D(38, 1, 1); _hVis_addBJet2Pt_abs = bookHisto1D(39, 1, 1); _hVis_addBJet2Pt = bookHisto1D(40, 1, 1); _hVis_addBJet2Eta_abs = bookHisto1D(41, 1, 1); _hVis_addBJet2Eta = bookHisto1D(42, 1, 1); _hVis_addBBMass_abs = bookHisto1D(43, 1, 1); _hVis_addBBMass = bookHisto1D(44, 1, 1); _hVis_addBBDR_abs = bookHisto1D(45, 1, 1); _hVis_addBBDR = bookHisto1D(46, 1, 1); _hFull_addBJet1Pt_abs = bookHisto1D(47, 1, 1); _hFull_addBJet1Pt = bookHisto1D(48, 1, 1); _hFull_addBJet1Eta_abs = bookHisto1D(49, 1, 1); _hFull_addBJet1Eta = bookHisto1D(50, 1, 1); _hFull_addBJet2Pt_abs = bookHisto1D(51, 1, 1); _hFull_addBJet2Pt = bookHisto1D(52, 1, 1); _hFull_addBJet2Eta_abs = bookHisto1D(53, 1, 1); _hFull_addBJet2Eta = bookHisto1D(54, 1, 1); _hFull_addBBMass_abs = bookHisto1D(55, 1, 1); _hFull_addBBMass = bookHisto1D(56, 1, 1); _hFull_addBBDR_abs = bookHisto1D(57, 1, 1); _hFull_addBBDR = bookHisto1D(58, 1, 1); _h_gap_addJet1Pt = bookProfile1D(59, 1, 1); _h_gap_addJet1Pt_eta0 = bookProfile1D(60, 1, 1); _h_gap_addJet1Pt_eta1 = bookProfile1D(61, 1, 1); _h_gap_addJet1Pt_eta2 = bookProfile1D(62, 1, 1); _h_gap_addJet2Pt = bookProfile1D(63, 1, 1); _h_gap_addJet2Pt_eta0 = bookProfile1D(64, 1, 1); _h_gap_addJet2Pt_eta1 = bookProfile1D(65, 1, 1); _h_gap_addJet2Pt_eta2 = bookProfile1D(66, 1, 1); _h_gap_addJetHT = bookProfile1D(67, 1, 1); _h_gap_addJetHT_eta0 = bookProfile1D(68, 1, 1); _h_gap_addJetHT_eta1 = bookProfile1D(69, 1, 1); _h_gap_addJetHT_eta2 = bookProfile1D(70, 1, 1); } void analyze(const Event& event) { // The objects used in the PAPER 12-041 are defined as follows (see p.16 for details): // // * Leptons : from the W boson decays after FSR // * Jets : anti-kT R=0.5 to all stable particles // exclude W->enu, munu, taunu // * B jet : B-Ghost matched // * B from top : B hadron from top->b decay // // Visible phase space definition: // // * Leptons : pT > 20, |eta| < 2.4 // * B jets from top : pT > 30, |eta| < 2.4 // Additional jets : pT > 20, |eta| < 2.4 // * // Full phase space definition: // // * Correction to dilepton BR from W boson BR // * No cut on top decay products // * Additional jets : pT > 20, |eta| < 2.4 // Do the analysis only for the ttbar full leptonic channel, removing tau decays const Particles partontops = apply(event, "PartonTops").particlesByPt(); if (partontops.size() != 2) vetoEvent; const Particle& t1 = partontops[0]; const Particle& t2 = partontops[1]; // Apply acceptance cuts on top-decay leptons (existence should be guaranteed) const auto isPromptChLepton = [](const Particle& p){return isChargedLepton(p) && !fromDecay(p);}; const Particle lep1 = t1.allDescendants(lastParticleWith(isPromptChLepton)).front(); const Particle lep2 = t2.allDescendants(lastParticleWith(isPromptChLepton)).front(); if (lep1.pT() < 1e-9*GeV || lep2.pT() < 1e-9*GeV) vetoEvent; // sanity check? const Jets jets = apply(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.4); int nJet30 = 0, nJet60 = 0, nJet100 = 0; Jets topBJets, addJets, addBJets, addJets_eta0, addJets_eta1, addJets_eta2; for (const Jet& jet : jets) { if (jet.pT() > 30*GeV) nJet30 += 1; if (jet.pT() > 60*GeV) nJet60 += 1; if (jet.pT() > 100*GeV) nJet100 += 1; const bool isBtagged = jet.bTagged(); - const bool isBFromTop = any(jet.bTags(), hasParticleAncestorWith(Cuts::abspid == PID::TQUARK)); + const bool isBFromTop = any(jet.bTags(), hasParticleAncestorWith(Cuts::abspid == PID::TQUARK, false)); if (isBFromTop) { if (jet.pT() > 30*GeV) topBJets.push_back(jet); } else { addJets.push_back(jet); if (isBtagged) addBJets.push_back(jet); if (jet.abseta() < 0.8 ) addJets_eta0.push_back(jet); else if (jet.abseta() < 1.5 ) addJets_eta1.push_back(jet); else if (jet.abseta() < 2.4 ) addJets_eta2.push_back(jet); } } const bool isVisiblePS = topBJets.size() >= 2 && lep1.pT() > 20*GeV && lep1.abseta() < 2.4 && lep2.pT() > 20*GeV && lep2.abseta() < 2.4; MSG_DEBUG(isVisiblePS << ": #b(top) = " << topBJets.size() << "; l1 = " << lep1.pT() << ", " << lep1.abseta() << "; l2 = " << lep2.pT() << ", " << lep2.abseta()); const double weight = event.weight(); if (isVisiblePS) { fillWithOF(_hVis_nJet30_abs, nJet30, weight); fillWithOF(_hVis_nJet30, nJet30, weight); fillWithOF(_hVis_nJet60_abs, nJet60, weight); fillWithOF(_hVis_nJet60, nJet60, weight); fillWithOF(_hVis_nJet100_abs, nJet100, weight); fillWithOF(_hVis_nJet100, nJet100, weight); fillGapFractions(addJets, _h_gap_addJet1Pt, _h_gap_addJet2Pt, _h_gap_addJetHT, weight); fillGapFractions(addJets_eta0, _h_gap_addJet1Pt_eta0, _h_gap_addJet2Pt_eta0, _h_gap_addJetHT_eta0, weight); fillGapFractions(addJets_eta1, _h_gap_addJet1Pt_eta1, _h_gap_addJet2Pt_eta1, _h_gap_addJetHT_eta1, weight); fillGapFractions(addJets_eta2, _h_gap_addJet1Pt_eta2, _h_gap_addJet2Pt_eta2, _h_gap_addJetHT_eta2, weight); } // Plots with two additional jets if (addJets.size() >= 1) { const double ht = sum(addJets, Kin::pT, 0.0); _hFull_addJJHT_abs->fill(ht/GeV, weight); _hFull_addJJHT ->fill(ht/GeV, weight); if (isVisiblePS) { _hVis_addJJHT_abs->fill(ht/GeV, weight); _hVis_addJJHT ->fill(ht/GeV, weight); } const Jet& j1 = addJets[0]; _hFull_addJet1Pt_abs ->fill(j1.pT()/GeV, weight); _hFull_addJet1Pt ->fill(j1.pT()/GeV, weight); _hFull_addJet1Eta_abs->fill(j1.abseta(), weight); _hFull_addJet1Eta ->fill(j1.abseta(), weight); if (isVisiblePS) { _hVis_addJet1Pt_abs ->fill(j1.pT()/GeV, weight); _hVis_addJet1Pt ->fill(j1.pT()/GeV, weight); _hVis_addJet1Eta_abs->fill(j1.abseta(), weight); _hVis_addJet1Eta ->fill(j1.abseta(), weight); } if (addJets.size() >= 2) { const Jet& j2 = addJets[1]; _hFull_addJet2Pt_abs ->fill(j2.pT()/GeV, weight); _hFull_addJet2Pt ->fill(j2.pT()/GeV, weight); _hFull_addJet2Eta_abs->fill(j2.abseta(), weight); _hFull_addJet2Eta ->fill(j2.abseta(), weight); if (isVisiblePS) { _hVis_addJet2Pt_abs ->fill(j2.pT()/GeV, weight); _hVis_addJet2Pt ->fill(j2.pT()/GeV, weight); _hVis_addJet2Eta_abs->fill(j2.abseta(), weight); _hVis_addJet2Eta ->fill(j2.abseta(), weight); } const double jjmass = (j1.mom() + j2.mom()).mass(); const double jjdR = deltaR(j1, j2); _hFull_addJJMass_abs->fill(jjmass/GeV, weight); _hFull_addJJMass ->fill(jjmass/GeV, weight); _hFull_addJJDR_abs ->fill(jjdR, weight); _hFull_addJJDR ->fill(jjdR, weight); if (isVisiblePS) { _hVis_addJJMass_abs->fill(jjmass/GeV, weight); _hVis_addJJMass ->fill(jjmass/GeV, weight); _hVis_addJJDR_abs ->fill(jjdR, weight); _hVis_addJJDR ->fill(jjdR, weight); } } } // Same set of plots if there are additional b-jets if (addBJets.size() >= 1) { const Jet& b1 = addBJets[0]; _hFull_addBJet1Pt_abs ->fill(b1.pT()/GeV, weight); _hFull_addBJet1Pt ->fill(b1.pT()/GeV, weight); _hFull_addBJet1Eta_abs->fill(b1.abseta(), weight); _hFull_addBJet1Eta ->fill(b1.abseta(), weight); if (isVisiblePS) { _hVis_addBJet1Pt_abs ->fill(b1.pT()/GeV, weight); _hVis_addBJet1Pt ->fill(b1.pT()/GeV, weight); _hVis_addBJet1Eta_abs->fill(b1.abseta(), weight); _hVis_addBJet1Eta ->fill(b1.abseta(), weight); } if (addBJets.size() >= 2) { const Jet& b2 = addBJets[1]; _hFull_addBJet2Pt_abs ->fill(b2.pT()/GeV, weight); _hFull_addBJet2Pt ->fill(b2.pT()/GeV, weight); _hFull_addBJet2Eta_abs->fill(b2.abseta(), weight); _hFull_addBJet2Eta ->fill(b2.abseta(), weight); if (isVisiblePS) { _hVis_addBJet2Pt_abs ->fill(b2.pT()/GeV, weight); _hVis_addBJet2Pt ->fill(b2.pT()/GeV, weight); _hVis_addBJet2Eta_abs->fill(b2.abseta(), weight); _hVis_addBJet2Eta ->fill(b2.abseta(), weight); } const double bbmass = (b1.mom() + b2.mom()).mass(); const double bbdR = deltaR(b1, b2); _hFull_addBBMass_abs->fill(bbmass/GeV, weight); _hFull_addBBMass ->fill(bbmass/GeV, weight); _hFull_addBBDR_abs ->fill(bbdR, weight); _hFull_addBBDR ->fill(bbdR, weight); if (isVisiblePS) { _hVis_addBBMass_abs->fill(bbmass/GeV, weight); _hVis_addBBMass ->fill(bbmass/GeV, weight); _hVis_addBBDR_abs ->fill(bbdR, weight); _hVis_addBBDR ->fill(bbdR, weight); } } } } void finalize() { const double ttbarXS = !std::isnan(crossSectionPerEvent()) ? crossSection() : 252.89*picobarn; if (std::isnan(crossSectionPerEvent())) MSG_INFO("No valid cross-section given, using NNLO (arXiv:1303.6254; sqrt(s)=8 TeV, m_t=172.5 GeV): " << ttbarXS/picobarn << " pb"); normalize({_hVis_nJet30,_hVis_nJet60, _hVis_nJet100, _hVis_addJet1Pt, _hVis_addJet1Eta, _hVis_addJet2Pt, _hVis_addJet2Eta, _hVis_addJJMass, _hVis_addJJDR, _hVis_addJJHT, _hFull_addJet1Pt, _hFull_addJet1Eta, _hFull_addJet2Pt, _hFull_addJet2Eta, _hFull_addJJMass, _hFull_addJJDR, _hFull_addJJHT, _hVis_addBJet1Pt, _hVis_addBJet1Eta, _hVis_addBJet2Pt, _hVis_addBJet2Eta, _hVis_addBBMass, _hVis_addBBDR, _hFull_addBJet1Pt, _hFull_addBJet1Eta, _hFull_addBJet2Pt, _hFull_addBJet2Eta, _hFull_addBBMass, _hFull_addBBDR}); const double xsPerWeight = ttbarXS/picobarn / sumOfWeights(); scale({_hVis_nJet30_abs, _hVis_nJet60_abs, _hVis_nJet100_abs, _hVis_addJet1Pt_abs, _hVis_addJet1Eta_abs, _hVis_addJet2Pt_abs, _hVis_addJet2Eta_abs, _hVis_addJJMass_abs, _hVis_addJJDR_abs, _hVis_addJJHT_abs, _hVis_addBJet1Pt_abs, _hVis_addBJet1Eta_abs, _hVis_addBJet2Pt_abs, _hVis_addBJet2Eta_abs, _hVis_addBBMass_abs, _hVis_addBBDR_abs}, xsPerWeight); const double sfull = xsPerWeight / 0.0454; //< correct for dilepton branching fraction scale({_hFull_addJet1Pt_abs, _hFull_addJet1Eta_abs, _hFull_addJet2Pt_abs, _hFull_addJet2Eta_abs, _hFull_addJJMass_abs, _hFull_addJJDR_abs, _hFull_addJJHT_abs, _hFull_addBJet1Pt_abs, _hFull_addBJet1Eta_abs, _hFull_addBJet2Pt_abs, _hFull_addBJet2Eta_abs, _hFull_addBBMass_abs, _hFull_addBBDR_abs}, sfull); } //@} void fillWithOF(Histo1DPtr h, double x, double w) { h->fill(std::min(x, h->xMax()-1e-9), w); } void fillGapFractions(const Jets& addJets, Profile1DPtr h_gap_addJet1Pt, Profile1DPtr h_gap_addJet2Pt, Profile1DPtr h_gap_addJetHT, double weight) { const double j1pt = (addJets.size() > 0) ? addJets[0].pT() : 0; for (size_t i = 0; i < h_gap_addJet1Pt->numBins(); ++i) { const double binCenter = h_gap_addJet1Pt->bin(i).xMid(); h_gap_addJet1Pt->fillBin(i, int(j1pt/GeV < binCenter), weight); } const double j2pt = (addJets.size() > 1) ? addJets[1].pT() : 0; for (size_t i = 0; i < h_gap_addJet2Pt->numBins(); ++i) { const double binCenter = h_gap_addJet2Pt->bin(i).xMid(); h_gap_addJet2Pt->fillBin(i, int(j2pt/GeV < binCenter), weight); } const double ht = sum(addJets, Kin::pT, 0.); for (size_t i = 0; i < h_gap_addJetHT->numBins(); ++i) { const double binCenter = h_gap_addJetHT->bin(i).xMid(); h_gap_addJetHT->fillBin(i, int(ht/GeV < binCenter) , weight); } } // @name Histogram data members //@{ Histo1DPtr _hVis_nJet30_abs, _hVis_nJet60_abs, _hVis_nJet100_abs; Histo1DPtr _hVis_addJet1Pt_abs, _hVis_addJet1Eta_abs, _hVis_addJet2Pt_abs, _hVis_addJet2Eta_abs; Histo1DPtr _hVis_addJJMass_abs, _hVis_addJJDR_abs, _hVis_addJJHT_abs; Histo1DPtr _hFull_addJet1Pt_abs, _hFull_addJet1Eta_abs, _hFull_addJet2Pt_abs, _hFull_addJet2Eta_abs; Histo1DPtr _hFull_addJJMass_abs, _hFull_addJJDR_abs, _hFull_addJJHT_abs; Histo1DPtr _hVis_addBJet1Pt_abs, _hVis_addBJet1Eta_abs, _hVis_addBJet2Pt_abs, _hVis_addBJet2Eta_abs; Histo1DPtr _hVis_addBBMass_abs, _hVis_addBBDR_abs; Histo1DPtr _hFull_addBJet1Pt_abs, _hFull_addBJet1Eta_abs, _hFull_addBJet2Pt_abs, _hFull_addBJet2Eta_abs; Histo1DPtr _hFull_addBBMass_abs, _hFull_addBBDR_abs; Histo1DPtr _hVis_nJet30, _hVis_nJet60, _hVis_nJet100; Histo1DPtr _hVis_addJet1Pt, _hVis_addJet1Eta, _hVis_addJet2Pt, _hVis_addJet2Eta; Histo1DPtr _hVis_addJJMass, _hVis_addJJDR, _hVis_addJJHT; Histo1DPtr _hFull_addJet1Pt, _hFull_addJet1Eta, _hFull_addJet2Pt, _hFull_addJet2Eta; Histo1DPtr _hFull_addJJMass, _hFull_addJJDR, _hFull_addJJHT; Histo1DPtr _hVis_addBJet1Pt, _hVis_addBJet1Eta, _hVis_addBJet2Pt, _hVis_addBJet2Eta; Histo1DPtr _hVis_addBBMass, _hVis_addBBDR; Histo1DPtr _hFull_addBJet1Pt, _hFull_addBJet1Eta, _hFull_addBJet2Pt, _hFull_addBJet2Eta; Histo1DPtr _hFull_addBBMass, _hFull_addBBDR; Profile1DPtr _h_gap_addJet1Pt, _h_gap_addJet1Pt_eta0, _h_gap_addJet1Pt_eta1, _h_gap_addJet1Pt_eta2; Profile1DPtr _h_gap_addJet2Pt, _h_gap_addJet2Pt_eta0, _h_gap_addJet2Pt_eta1, _h_gap_addJet2Pt_eta2; Profile1DPtr _h_gap_addJetHT, _h_gap_addJetHT_eta0, _h_gap_addJetHT_eta1, _h_gap_addJetHT_eta2; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2015_I1397174); } diff --git a/include/Rivet/Tools/ParticleUtils.hh b/include/Rivet/Tools/ParticleUtils.hh --- a/include/Rivet/Tools/ParticleUtils.hh +++ b/include/Rivet/Tools/ParticleUtils.hh @@ -1,778 +1,782 @@ #ifndef RIVET_PARTICLEUTILS_HH #define RIVET_PARTICLEUTILS_HH #include "Rivet/Particle.hh" #include "Rivet/Tools/ParticleBaseUtils.hh" #include "Rivet/Tools/ParticleIdUtils.hh" // Macros to map Rivet::Particle functions to PID:: functions of the same name #define PARTICLE_TO_PID_BOOLFN(fname) inline bool fname (const Particle& p) { return PID:: fname (p.pid()); } #define PARTICLE_TO_PID_INTFN(fname) inline int fname (const Particle& p) { return PID:: fname (p.pid()); } #define PARTICLE_TO_PID_DBLFN(fname) inline double fname (const Particle& p) { return PID:: fname (p.pid()); } namespace Rivet { /// @name Particle classifier functions //@{ /// Unbound function access to PID code inline int pid(const Particle& p) { return p.pid(); } /// Unbound function access to abs PID code inline int abspid(const Particle& p) { return p.abspid(); } /// Is this particle species charged? PARTICLE_TO_PID_BOOLFN(isCharged) /// Is this particle species neutral? PARTICLE_TO_PID_BOOLFN(isNeutral) /// Is this a neutrino? PARTICLE_TO_PID_BOOLFN(isNeutrino) /// Determine if the PID is that of a charged lepton PARTICLE_TO_PID_BOOLFN(isChargedLepton) PARTICLE_TO_PID_BOOLFN(isChLepton) /// Determine if the PID is that of a lepton (charged or neutral) PARTICLE_TO_PID_BOOLFN(isLepton) /// Determine if the PID is that of a photon PARTICLE_TO_PID_BOOLFN(isPhoton) /// Determine if the PID is that of an electron or positron PARTICLE_TO_PID_BOOLFN(isElectron) /// Determine if the PID is that of an muon or antimuon PARTICLE_TO_PID_BOOLFN(isMuon) /// Determine if the PID is that of an tau or antitau PARTICLE_TO_PID_BOOLFN(isTau) /// Determine if the PID is that of a hadron PARTICLE_TO_PID_BOOLFN(isHadron) /// Determine if the PID is that of a meson PARTICLE_TO_PID_BOOLFN(isMeson) /// Determine if the PID is that of a baryon PARTICLE_TO_PID_BOOLFN(isBaryon) /// Determine if the PID is that of a quark PARTICLE_TO_PID_BOOLFN(isQuark) /// Determine if the PID is that of a parton (quark or gluon) PARTICLE_TO_PID_BOOLFN(isParton) /// Determine if the PID is that of a W+ PARTICLE_TO_PID_BOOLFN(isWplus) /// Determine if the PID is that of a W- PARTICLE_TO_PID_BOOLFN(isWminus) /// Determine if the PID is that of a W+- PARTICLE_TO_PID_BOOLFN(isW) /// Determine if the PID is that of a Z0 PARTICLE_TO_PID_BOOLFN(isZ) /// Determine if the PID is that of an SM/lightest SUSY Higgs PARTICLE_TO_PID_BOOLFN(isHiggs) /// Determine if the PID is that of an s/sbar PARTICLE_TO_PID_BOOLFN(isStrange) /// Determine if the PID is that of a c/cbar PARTICLE_TO_PID_BOOLFN(isCharm) /// Determine if the PID is that of a b/bbar PARTICLE_TO_PID_BOOLFN(isBottom) /// Determine if the PID is that of a t/tbar PARTICLE_TO_PID_BOOLFN(isTop) /// Determine if the particle is a heavy flavour hadron or parton PARTICLE_TO_PID_BOOLFN(isHeavyFlavour) /// Determine if the PID is that of a heavy parton (c,b,t) PARTICLE_TO_PID_BOOLFN(isHeavyParton) /// Determine if the PID is that of a light parton (u,d,s) PARTICLE_TO_PID_BOOLFN(isLightParton) /// Determine if the PID is that of a heavy flavour (b or c) meson PARTICLE_TO_PID_BOOLFN(isHeavyMeson) /// Determine if the PID is that of a heavy flavour (b or c) baryon PARTICLE_TO_PID_BOOLFN(isHeavyBaryon) /// Determine if the PID is that of a heavy flavour (b or c) hadron PARTICLE_TO_PID_BOOLFN(isHeavyHadron) /// Determine if the PID is that of a light flavour (not b or c) meson PARTICLE_TO_PID_BOOLFN(isLightMeson) /// Determine if the PID is that of a light flavour (not b or c) baryon PARTICLE_TO_PID_BOOLFN(isLightBaryon) /// Determine if the PID is that of a light flavour (not b or c) hadron PARTICLE_TO_PID_BOOLFN(isLightHadron) /// Determine if the PID is that of a b-meson. PARTICLE_TO_PID_BOOLFN(isBottomMeson) /// Determine if the PID is that of a b-baryon. PARTICLE_TO_PID_BOOLFN(isBottomBaryon) /// Determine if the PID is that of a b-hadron. PARTICLE_TO_PID_BOOLFN(isBottomHadron) /// @brief Determine if the PID is that of a c-meson. /// /// Specifically, the _heaviest_ quark is a c: a B_c is a b-meson and NOT a c-meson. /// Charmonia (closed charm) are counted as c-mesons here. PARTICLE_TO_PID_BOOLFN(isCharmMeson) /// @brief Determine if the PID is that of a c-baryon. /// /// Specifically, the _heaviest_ quark is a c: a baryon containing a b & c /// is a b-baryon and NOT a c-baryon. To test for the simpler case, just use /// a combination of hasCharm() and isBaryon(). PARTICLE_TO_PID_BOOLFN(isCharmBaryon) /// Determine if the PID is that of a c-hadron. PARTICLE_TO_PID_BOOLFN(isCharmHadron) // /// Determine if the PID is that of a strange meson // PARTICLE_TO_PID_BOOLFN(isStrangeMeson) // /// Determine if the PID is that of a strange baryon // PARTICLE_TO_PID_BOOLFN(isStrangeBaryon) // /// Determine if the PID is that of a strange hadron // PARTICLE_TO_PID_BOOLFN(isStrangeHadron) /// Is this a pomeron, odderon, or generic reggeon? PARTICLE_TO_PID_BOOLFN(isReggeon) /// Determine if the PID is that of a diquark (used in hadronization models) PARTICLE_TO_PID_BOOLFN(isDiquark) /// Determine if the PID is that of a pentaquark (hypothetical hadron) PARTICLE_TO_PID_BOOLFN(isPentaquark) /// Is this a fundamental SUSY particle? PARTICLE_TO_PID_BOOLFN(isSUSY) /// Is this an R-hadron? PARTICLE_TO_PID_BOOLFN(isRhadron) /// Is this a technicolor particle? PARTICLE_TO_PID_BOOLFN(isTechnicolor) /// Is this an excited (composite) quark or lepton? PARTICLE_TO_PID_BOOLFN(isExcited) /// Is this a Kaluza-Klein excitation? PARTICLE_TO_PID_BOOLFN(isKK) /// Is this a graviton? PARTICLE_TO_PID_BOOLFN(isGraviton) /// Is this a BSM particle (including graviton)? PARTICLE_TO_PID_BOOLFN(isBSM) /// Determine if the PID is in the generator-specific range PARTICLE_TO_PID_BOOLFN(isGenSpecific) /// Determine if the PID is that of an EW scale resonance PARTICLE_TO_PID_BOOLFN(isResonance) /// Check the PID for usability in transport codes like Geant4 PARTICLE_TO_PID_BOOLFN(isTransportable) /// Does this particle contain an up quark? PARTICLE_TO_PID_BOOLFN(hasUp) /// Does this particle contain a down quark? PARTICLE_TO_PID_BOOLFN(hasDown) /// Does this particle contain a strange quark? PARTICLE_TO_PID_BOOLFN(hasStrange) /// Does this particle contain a charm quark? PARTICLE_TO_PID_BOOLFN(hasCharm) /// Does this particle contain a bottom quark? PARTICLE_TO_PID_BOOLFN(hasBottom) /// Does this particle contain a top quark? PARTICLE_TO_PID_BOOLFN(hasTop) /// jSpin returns 2J+1, where J is the total spin PARTICLE_TO_PID_INTFN(jSpin) /// sSpin returns 2S+1, where S is the spin PARTICLE_TO_PID_INTFN(sSpin) /// lSpin returns 2L+1, where L is the orbital angular momentum PARTICLE_TO_PID_INTFN(lSpin) /// Return the charge PARTICLE_TO_PID_DBLFN(charge) /// Return 3 times the charge (3 x quark charge is an int) PARTICLE_TO_PID_INTFN(charge3) /// Return the absolute charge PARTICLE_TO_PID_DBLFN(abscharge) /// Return 3 times the abs charge (3 x quark charge is an int) PARTICLE_TO_PID_INTFN(abscharge3) /// Alias for charge3 /// @deprecated Use charge3 PARTICLE_TO_PID_INTFN(threeCharge) /// Get the atomic number (number of protons) in a nucleus/ion PARTICLE_TO_PID_INTFN(nuclZ) /// Get the atomic weight (number of nucleons) in a nucleus/ion PARTICLE_TO_PID_INTFN(nuclA) /// If this is a nucleus (ion), get nLambda PARTICLE_TO_PID_INTFN(nuclNlambda) //@} /// @name Particle pair classifiers /// @todo Make versions that work on ParticlePair? //@{ inline bool isSameSign(const Particle& a, const Particle& b) { return PID::isSameSign(a.pid(), b.pid()); } inline bool isOppSign(const Particle& a, const Particle& b) { return PID::isOppSign(a.pid(), b.pid()); } inline bool isSameFlav(const Particle& a, const Particle& b) { return PID::isSameFlav(a.pid(), b.pid()); } inline bool isOppFlav(const Particle& a, const Particle& b) { return PID::isOppFlav(a.pid(), b.pid()); } inline bool isOSSF(const Particle& a, const Particle& b) { return PID::isOSSF(a.pid(), b.pid()); } inline bool isSSSF(const Particle& a, const Particle& b) { return PID::isSSSF(a.pid(), b.pid()); } inline bool isOSOF(const Particle& a, const Particle& b) { return PID::isOSOF(a.pid(), b.pid()); } inline bool isSSOF(const Particle& a, const Particle& b) { return PID::isSSOF(a.pid(), b.pid()); } //@} /// @name Particle charge/sign comparison functions //@{ /// @brief Return true if Particles @a a and @a b have the opposite charge sign /// @note Two neutrals returns false inline bool oppSign(const Particle& a, const Particle& b) { return sign(a.charge3()) == -sign(b.charge3()) && sign(a.charge3()) != ZERO; } /// Return true if Particles @a a and @a b have the same charge sign /// @note Two neutrals returns true inline bool sameSign(const Particle& a, const Particle& b) { return sign(a.charge3()) == sign(b.charge3()); } /// Return true if Particles @a a and @a b have the exactly opposite charge /// @note Two neutrals returns false inline bool oppCharge(const Particle& a, const Particle& b) { return a.charge3() == -b.charge3() && a.charge3() != 0; } /// Return true if Particles @a a and @a b have the same charge (including neutral) /// @note Two neutrals returns true inline bool sameCharge(const Particle& a, const Particle& b) { return a.charge3() == b.charge3(); } /// Return true if Particles @a a and @a b have a different (not necessarily opposite) charge inline bool diffCharge(const Particle& a, const Particle& b) { return a.charge3() != b.charge3(); } //@} ////////////////////////////////////// /// @name Non-PID particle properties, via unbound functions //@{ /// @brief Determine whether a particle is the first in a decay chain to meet the function requirement inline bool isFirstWith(const Particle& p, const ParticleSelector& f) { return p.isFirstWith(f); } /// @brief Determine whether a particle is the first in a decay chain not to meet the function requirement inline bool isFirstWithout(const Particle& p, const ParticleSelector& f) { return p.isFirstWithout(f); } /// @brief Determine whether a particle is the last in a decay chain to meet the function requirement inline bool isLastWith(const Particle& p, const ParticleSelector& f) { return p.isLastWith(f); } /// @brief Determine whether a particle is the last in a decay chain not to meet the function requirement inline bool isLastWithout(const Particle& p, const ParticleSelector& f) { return p.isLastWithout(f); } /// @brief Determine whether a particle has an ancestor which meets the function requirement inline bool hasAncestorWith(const Particle& p, const ParticleSelector& f, bool only_physical=true) { return p.hasAncestorWith(f, only_physical); } /// @brief Determine whether a particle has an ancestor which doesn't meet the function requirement inline bool hasAncestorWithout(const Particle& p, const ParticleSelector& f, bool only_physical=true) { return p.hasAncestorWithout(f, only_physical); } /// @brief Determine whether a particle has a parent which meets the function requirement inline bool hasParentWith(const Particle& p, const ParticleSelector& f) { return p.hasParentWith(f); } /// @brief Determine whether a particle has a parent which doesn't meet the function requirement inline bool hasParentWithout(const Particle& p, const ParticleSelector& f) { return p.hasParentWithout(f); } /// @brief Determine whether a particle has a child which meets the function requirement inline bool hasChildWith(const Particle& p, const ParticleSelector& f) { return p.hasChildWith(f); } /// @brief Determine whether a particle has a child which doesn't meet the function requirement inline bool hasChildWithout(const Particle& p, const ParticleSelector& f) { return p.hasChildWithout(f); } /// @brief Determine whether a particle has a descendant which meets the function requirement inline bool hasDescendantWith(const Particle& p, const ParticleSelector& f, bool remove_duplicates=true) { return p.hasDescendantWith(f, remove_duplicates); } /// @brief Determine whether a particle has a descendant which doesn't meet the function requirement inline bool hasDescendantWithout(const Particle& p, const ParticleSelector& f, bool remove_duplicates=true) { return p.hasDescendantWithout(f, remove_duplicates); } /// @brief Determine whether a particle has a stable descendant which meets the function requirement inline bool hasStableDescendantWith(const Particle& p, const ParticleSelector& f) { return p.hasStableDescendantWith(f); } /// @brief Determine whether a particle has a stable descendant which doesn't meet the function requirement inline bool hasStableDescendantWithout(const Particle& p, const ParticleSelector& f) { return p.hasStableDescendantWithout(f); } /// Is this particle potentially visible in a detector? inline bool isVisible(const Particle& p) { return p.isVisible(); } /// @brief Decide if a given particle is direct, via Particle::isDirect() /// /// 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. inline bool isDirect(const Particle& p, bool allow_from_direct_tau=false, bool allow_from_direct_mu=false) { return p.isDirect(allow_from_direct_tau, allow_from_direct_mu); } /// @brief Decide if a given particle is prompt, via Particle::isPrompt() /// /// The boolean arguments allow a decay lepton to be considered prompt if /// its parent was a "real" prompt lepton. inline bool isPrompt(const Particle& p, bool allow_from_prompt_tau=false, bool allow_from_prompt_mu=false) { return p.isPrompt(allow_from_prompt_tau, allow_from_prompt_mu); } /// Decide if a given particle is stable, via Particle::isStable() inline bool isStable(const Particle& p) { return p.isStable(); } /// Decide if a given particle decays hadronically inline bool hasHadronicDecay(const Particle& p) { if (p.isStable()) return false; if (p.hasChildWith(isHadron)) return true; return false; } /// Decide if a given particle decays leptonically (decays, and no hadrons) inline bool hasLeptonicDecay(const Particle& p) { if (p.isStable()) return false; if (p.hasChildWith(isHadron)) return false; return true; } /// Check whether a given PID is found in the particle's ancestor list /// @deprecated Prefer hasAncestorWith inline bool hasAncestor(const Particle& p, PdgId pid) { return p.hasAncestor(pid); } /// Determine whether the particle is from a b-hadron decay inline bool fromBottom(const Particle& p) { return p.fromBottom(); } /// @brief Determine whether the particle is from a c-hadron decay inline bool fromCharm(const Particle& p) { return p.fromCharm(); } /// @brief Determine whether the particle is from a hadron decay inline bool fromHadron(const Particle& p) { return p.fromHadron(); } /// @brief Determine whether the particle is from a tau decay inline bool fromTau(const Particle& p, bool prompt_taus_only=false) { return p.fromTau(prompt_taus_only); } /// @brief Determine whether the particle is from a prompt tau decay inline bool fromPromptTau(const Particle& p) { return p.fromPromptTau(); } /// @brief Determine whether the particle is from a hadron or tau decay /// @deprecated Too vague: use fromHadron or fromHadronicTau inline bool fromDecay(const Particle& p) { return p.fromDecay(); } //@} /// @name Particle classifier -> bool functors /// /// To be passed to any() or all() e.g. any(p.children(), HasPID(PID::MUON)) //@{ /// Base type for Particle -> bool functors struct BoolParticleFunctor { virtual bool operator()(const Particle& p) const = 0; virtual ~BoolParticleFunctor() {} }; /// Functor for and-combination of selector logic struct BoolParticleAND : public BoolParticleFunctor { BoolParticleAND(const std::vector& sels) : selectors(sels) {} BoolParticleAND(const ParticleSelector& a, const ParticleSelector& b) : selectors({a,b}) {} BoolParticleAND(const ParticleSelector& a, const ParticleSelector& b, const ParticleSelector& c) : selectors({a,b,c}) {} bool operator()(const Particle& p) const { for (const ParticleSelector& sel : selectors) if (!sel(p)) return false; return true; } std::vector selectors; }; /// Operator syntactic sugar for AND construction inline BoolParticleAND operator && (const ParticleSelector& a, const ParticleSelector& b) { return BoolParticleAND(a, b); } /// Functor for or-combination of selector logic struct BoolParticleOR : public BoolParticleFunctor { BoolParticleOR(const std::vector& sels) : selectors(sels) {} BoolParticleOR(const ParticleSelector& a, const ParticleSelector& b) : selectors({a,b}) {} BoolParticleOR(const ParticleSelector& a, const ParticleSelector& b, const ParticleSelector& c) : selectors({a,b,c}) {} bool operator()(const Particle& p) const { for (const ParticleSelector& sel : selectors) if (sel(p)) return true; return false; } std::vector selectors; }; /// Operator syntactic sugar for OR construction inline BoolParticleOR operator || (const ParticleSelector& a, const ParticleSelector& b) { return BoolParticleOR(a, b); } /// Functor for inverting selector logic struct BoolParticleNOT : public BoolParticleFunctor { BoolParticleNOT(const ParticleSelector& sel) : selector(sel) {} bool operator()(const Particle& p) const { return !selector(p); } ParticleSelector selector; }; /// Operator syntactic sugar for NOT construction inline BoolParticleNOT operator ! (const ParticleSelector& a) { return BoolParticleNOT(a); } /// PID matching functor struct HasPID : public BoolParticleFunctor { HasPID(PdgId pid) : targetpids{pid} { } HasPID(vector pids) : targetpids{pids} { } HasPID(initializer_list pids) : targetpids{pids} { } bool operator()(const Particle& p) const { return contains(targetpids, p.pid()); } vector targetpids; }; using hasPID = HasPID; /// |PID| matching functor struct HasAbsPID : public BoolParticleFunctor { HasAbsPID(PdgId pid) : targetapids{abs(pid)} { } HasAbsPID(vector pids) { for (PdgId pid : pids) targetapids.push_back(abs(pid)); } HasAbsPID(initializer_list pids) { for (PdgId pid : pids) targetapids.push_back(abs(pid)); } bool operator()(const Particle& p) const { return contains(targetapids, p.abspid()); } vector targetapids; }; using hasAbsPID = HasAbsPID; /// Determine whether a particle is the first in a decay chain to meet the cut/function struct FirstParticleWith : public BoolParticleFunctor { FirstParticleWith(const ParticleSelector& f) : fn(f) { } FirstParticleWith(const Cut& c); bool operator()(const Particle& p) const { return isFirstWith(p, fn); } ParticleSelector fn; }; using firstParticleWith = FirstParticleWith; /// Determine whether a particle is the first in a decay chain not to meet the cut/function struct FirstParticleWithout : public BoolParticleFunctor { FirstParticleWithout(const ParticleSelector& f) : fn(f) { } FirstParticleWithout(const Cut& c); bool operator()(const Particle& p) const { return isFirstWithout(p, fn); } ParticleSelector fn; }; using firstParticleWithout = FirstParticleWithout; /// Determine whether a particle is the last in a decay chain to meet the cut/function struct LastParticleWith : public BoolParticleFunctor { template LastParticleWith(const FN& f) : fn(f) { } LastParticleWith(const Cut& c); bool operator()(const Particle& p) const { return isLastWith(p, fn); } std::function fn; }; using lastParticleWith = LastParticleWith; /// Determine whether a particle is the last in a decay chain not to meet the cut/function struct LastParticleWithout : public BoolParticleFunctor { LastParticleWithout(const ParticleSelector& f) : fn(f) { } LastParticleWithout(const Cut& c); bool operator()(const Particle& p) const { return isLastWithout(p, fn); } ParticleSelector fn; }; using lastParticleWithout = LastParticleWithout; /// Determine whether a particle has an ancestor which meets the cut/function struct HasParticleAncestorWith : public BoolParticleFunctor { - HasParticleAncestorWith(const ParticleSelector& f) : fn(f) { } - HasParticleAncestorWith(const Cut& c); - bool operator()(const Particle& p) const { return hasAncestorWith(p, fn); } + HasParticleAncestorWith(const ParticleSelector& f, bool only_physical=true) : fn(f), onlyphysical(only_physical) { } + HasParticleAncestorWith(const Cut& c, bool only_physical=true); + bool operator()(const Particle& p) const { return hasAncestorWith(p, fn, onlyphysical); } ParticleSelector fn; + bool onlyphysical; }; using hasParticleAncestorWith = HasParticleAncestorWith; /// Determine whether a particle has an ancestor which doesn't meet the cut/function struct HasParticleAncestorWithout : public BoolParticleFunctor { - HasParticleAncestorWithout(const ParticleSelector& f) : fn(f) { } - HasParticleAncestorWithout(const Cut& c); - bool operator()(const Particle& p) const { return hasAncestorWithout(p, fn); } + HasParticleAncestorWithout(const ParticleSelector& f, bool only_physical=true) : fn(f), onlyphysical(only_physical) { } + HasParticleAncestorWithout(const Cut& c, bool only_physical=true); + bool operator()(const Particle& p) const { return hasAncestorWithout(p, fn, onlyphysical); } ParticleSelector fn; + bool onlyphysical; }; using hasParticleAncestorWithout = HasParticleAncestorWithout; /// Determine whether a particle has an parent which meets the cut/function struct HasParticleParentWith : public BoolParticleFunctor { HasParticleParentWith(const ParticleSelector& f) : fn(f) { } HasParticleParentWith(const Cut& c); bool operator()(const Particle& p) const { return hasParentWith(p, fn); } ParticleSelector fn; }; using hasParticleParentWith = HasParticleParentWith; /// Determine whether a particle has an parent which doesn't meet the cut/function struct HasParticleParentWithout : public BoolParticleFunctor { HasParticleParentWithout(const ParticleSelector& f) : fn(f) { } HasParticleParentWithout(const Cut& c); bool operator()(const Particle& p) const { return hasParentWithout(p, fn); } ParticleSelector fn; }; using hasParticleParentWithout = HasParticleParentWithout; /// Determine whether a particle has a child which meets the cut/function struct HasParticleChildWith : public BoolParticleFunctor { HasParticleChildWith(const ParticleSelector& f) : fn(f) { } HasParticleChildWith(const Cut& c); bool operator()(const Particle& p) const { return hasChildWith(p, fn); } ParticleSelector fn; }; using hasParticleChildWith = HasParticleChildWith; /// Determine whether a particle has a child which doesn't meet the cut/function struct HasParticleChildWithout : public BoolParticleFunctor { HasParticleChildWithout(const ParticleSelector& f) : fn(f) { } HasParticleChildWithout(const Cut& c); bool operator()(const Particle& p) const { return hasChildWithout(p, fn); } ParticleSelector fn; }; using hasParticleChildWithout = HasParticleChildWithout; /// Determine whether a particle has a descendant which meets the cut/function struct HasParticleDescendantWith : public BoolParticleFunctor { - HasParticleDescendantWith(const ParticleSelector& f) : fn(f) { } - HasParticleDescendantWith(const Cut& c); - bool operator()(const Particle& p) const { return hasDescendantWith(p, fn); } + HasParticleDescendantWith(const ParticleSelector& f, bool remove_duplicates=true) : fn(f), rmduplicates(remove_duplicates) { } + HasParticleDescendantWith(const Cut& c, bool remove_duplicates=true); + bool operator()(const Particle& p) const { return hasDescendantWith(p, fn, rmduplicates); } ParticleSelector fn; + bool rmduplicates; }; using hasParticleDescendantWith = HasParticleDescendantWith; /// Determine whether a particle has a descendant which doesn't meet the cut/function struct HasParticleDescendantWithout : public BoolParticleFunctor { - HasParticleDescendantWithout(const ParticleSelector& f) : fn(f) { } - HasParticleDescendantWithout(const Cut& c); - bool operator()(const Particle& p) const { return hasDescendantWithout(p, fn); } + HasParticleDescendantWithout(const ParticleSelector& f, bool remove_duplicates=true) : fn(f), rmduplicates(remove_duplicates) { } + HasParticleDescendantWithout(const Cut& c, bool remove_duplicates=true); + bool operator()(const Particle& p) const { return hasDescendantWithout(p, fn, rmduplicates); } ParticleSelector fn; + bool rmduplicates; }; using hasParticleDescendantWithout = HasParticleDescendantWithout; //@} /// @name Unbound functions for filtering particles //@{ /// Filter a particle collection in-place to the subset that passes the supplied Cut Particles& ifilter_select(Particles& particles, const Cut& c); /// Alias for ifilter_select /// @deprecated Use ifilter_select inline Particles& ifilterBy(Particles& particles, const Cut& c) { return ifilter_select(particles, c); } /// Filter a particle collection in-place to the subset that passes the supplied Cut inline Particles filter_select(const Particles& particles, const Cut& c) { Particles rtn = particles; return ifilter_select(rtn, c); } /// Alias for ifilter_select /// @deprecated Use filter_select inline Particles filterBy(const Particles& particles, const Cut& c) { return filter_select(particles, c); } /// Filter a particle collection in-place to the subset that passes the supplied Cut inline Particles filter_select(const Particles& particles, const Cut& c, Particles& out) { out = filter_select(particles, c); return out; } /// Alias for ifilter_select /// @deprecated Use filter_select inline Particles filterBy(const Particles& particles, const Cut& c, Particles& out) { return filter_select(particles, c, out); } /// Filter a particle collection in-place to the subset that fails the supplied Cut Particles& ifilter_discard(Particles& particles, const Cut& c); /// Filter a particle collection in-place to the subset that fails the supplied Cut inline Particles filter_discard(const Particles& particles, const Cut& c) { Particles rtn = particles; return ifilter_discard(rtn, c); } /// Filter a particle collection in-place to the subset that fails the supplied Cut inline Particles filter_discard(const Particles& particles, const Cut& c, Particles& out) { out = filter_discard(particles, c); return out; } // inline void ifilterIsolateDeltaR(Particles& particles, const FourMomenta& vecs) { // ifilter_discard(particles, // } // inline Particles filterIsolateDeltaR(const Particles& particles, const FourMomenta& vecs) { // } //@} /// @name Particle pair functions //@{ /// Get the PDG ID codes of a ParticlePair /// @todo Make ParticlePair a custom class instead? inline PdgIdPair pids(const ParticlePair& pp) { return make_pair(pp.first.pid(), pp.second.pid()); } //@} /// @name Operations on collections of Particle /// @note This can't be done on generic collections of ParticleBase -- thanks, C++ :-/ //@{ namespace Kin { inline double sumPt(const Particles& ps) { return sum(ps, pT, 0.0); } inline FourMomentum sumP4(const Particles& ps) { return sum(ps, p4, FourMomentum()); } inline Vector3 sumP3(const Particles& ps) { return sum(ps, p3, Vector3()); } /// @todo Min dPhi, min dR? /// @todo Isolation routines? } //@} // Import Kin namespace into Rivet using namespace Kin; /// Check Particle equivalence inline bool isSame(const Particle& a, const Particle& b) { return a.isSame(b); } } #endif diff --git a/src/Tools/ParticleUtils.cc b/src/Tools/ParticleUtils.cc --- a/src/Tools/ParticleUtils.cc +++ b/src/Tools/ParticleUtils.cc @@ -1,58 +1,62 @@ #include "Rivet/Tools/ParticleUtils.hh" #include "Rivet/Tools/Cuts.hh" namespace Rivet { FirstParticleWith::FirstParticleWith(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } FirstParticleWithout::FirstParticleWithout(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } LastParticleWith::LastParticleWith(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } LastParticleWithout::LastParticleWithout(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } - HasParticleAncestorWith::HasParticleAncestorWith(const Cut& c) - : fn([&](const Particle& p){ return c->accept(p); }) { } + HasParticleAncestorWith::HasParticleAncestorWith(const Cut& c, bool only_physical) + : fn([&](const Particle& p){ return c->accept(p); }), + onlyphysical(only_physical) { } - HasParticleAncestorWithout::HasParticleAncestorWithout(const Cut& c) - : fn([&](const Particle& p){ return c->accept(p); }) { } + HasParticleAncestorWithout::HasParticleAncestorWithout(const Cut& c, bool only_physical) + : fn([&](const Particle& p){ return c->accept(p); }), + onlyphysical(only_physical) { } HasParticleParentWith::HasParticleParentWith(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } HasParticleParentWithout::HasParticleParentWithout(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } HasParticleChildWith::HasParticleChildWith(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } HasParticleChildWithout::HasParticleChildWithout(const Cut& c) : fn([&](const Particle& p){ return c->accept(p); }) { } - HasParticleDescendantWith::HasParticleDescendantWith(const Cut& c) - : fn([&](const Particle& p){ return c->accept(p); }) { } + HasParticleDescendantWith::HasParticleDescendantWith(const Cut& c, bool remove_duplicates) + : fn([&](const Particle& p){ return c->accept(p); }), + rmduplicates(remove_duplicates) { } - HasParticleDescendantWithout::HasParticleDescendantWithout(const Cut& c) - : fn([&](const Particle& p){ return c->accept(p); }) { } + HasParticleDescendantWithout::HasParticleDescendantWithout(const Cut& c, bool remove_duplicates) + : fn([&](const Particle& p){ return c->accept(p); }), + rmduplicates(remove_duplicates) { } Particles& ifilter_select(Particles& particles, const Cut& c) { if (c == Cuts::OPEN) return particles; // return ifilter_select(particles, *c); return ifilter_select(particles, [&](const Particle& p){return c->accept(p);}); } Particles& ifilter_discard(Particles& particles, const Cut& c) { if (c == Cuts::OPEN) { particles.clear(); return particles; } // return ifilter_discard(particles, *c); return ifilter_discard(particles, [&](const Particle& p){return c->accept(p);}); } }