diff --git a/analyses/pluginATLAS/ATLAS_2017_I1591327.cc b/analyses/pluginATLAS/ATLAS_2017_I1591327.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1591327.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1591327.cc @@ -1,186 +1,183 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// Isolated diphoton + X differential cross-sections class ATLAS_2017_I1591327 : public Analysis { public: // Constructor - ATLAS_2017_I1591327() : Analysis("ATLAS_2017_I1591327") { - - } + DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1591327); // Book histograms and initialise projections before the run void init() { const FinalState fs; declare(fs, "FS"); FastJets fj(fs, FastJets::KT, 0.5); _area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec()); fj.useJetArea(_area_def); declare(fj, "KtJetsD05"); IdentifiedFinalState photonfs(Cuts::abseta < 2.37 && Cuts::pT > 30*GeV); photonfs.acceptId(PID::PHOTON); declare(photonfs, "Photon"); // Histograms book(_h_M , 2, 1, 1); book(_h_pT , 3, 1, 1); book(_h_at , 4, 1, 1); book(_h_phistar, 5, 1, 1); book(_h_costh , 6, 1, 1); book(_h_dPhi , 7, 1, 1); } // Perform the per-event analysis void analyze(const Event& event) { // Require at least 2 photons in final state const Particles photons = apply(event, "Photon").particlesByPt(); if (photons.size() < 2) vetoEvent; // Compute the median energy density _ptDensity.clear(); _sigma.clear(); _Njets.clear(); vector > ptDensities; vector emptyVec; ptDensities.assign(ETA_BINS.size()-1, emptyVec); // Get jets, and corresponding jet areas const shared_ptr clust_seq_area = applyProjection(event, "KtJetsD05").clusterSeqArea(); for (const fastjet::PseudoJet& jet : apply(event, "KtJetsD05").pseudoJets(0.0*GeV)) { const double aeta = fabs(jet.eta()); const double pt = jet.perp(); const double area = clust_seq_area->area(jet); if (area < 1e-3) continue; const int ieta = binIndex(aeta, ETA_BINS); if (ieta != -1) ptDensities[ieta].push_back(pt/area); } // Compute median jet properties over the jets in the event for (size_t b = 0; b < ETA_BINS.size()-1; ++b) { double median = 0.0, sigma = 0.0; int Njets = 0; if (ptDensities[b].size() > 0) { std::sort(ptDensities[b].begin(), ptDensities[b].end()); int nDens = ptDensities[b].size(); median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2]; sigma = ptDensities[b][(int)(.15865*nDens)]; Njets = nDens; } _ptDensity.push_back(median); _sigma.push_back(sigma); _Njets.push_back(Njets); } // Loop over photons and fill vector of isolated ones Particles isolated_photons; for (const Particle& photon : photons) { // Check if it's a prompt photon (needed for SHERPA 2->5 sample, otherwise I also get photons from hadron decays in jets) - if (photon.fromDecay()) continue; + if (!photon.isPrompt()) continue; // Remove photons in ECAL crack region if (inRange(photon.abseta(), 1.37, 1.56)) continue; const double eta_P = photon.eta(); const double phi_P = photon.phi(); // Compute isolation via particles within an R=0.4 cone of the photon const Particles fs = apply(event, "FS").particles(); FourMomentum mom_in_EtCone; for (const Particle& p : fs) { // Reject if not in cone if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue; // Reject if in the 5x7 cell central core if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 && fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5) continue; // Sum momentum mom_in_EtCone += p.momentum(); } // Now figure out the correction (area*density) const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle const double correction = _ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area; // Discard the photon if there is more than 11 GeV of cone activity // NOTE: Shouldn't need to subtract photon itself (it's in the central core) if (mom_in_EtCone.Et() - correction > 11*GeV) continue; // Add isolated photon to list isolated_photons.push_back(photon); } // Require at least two isolated photons if (isolated_photons.size() < 2) vetoEvent; // Select leading pT pair sortByPt(isolated_photons); const FourMomentum y1 = isolated_photons[0]; const FourMomentum y2 = isolated_photons[1]; // Leading photon should have pT > 40 GeV, subleading > 30 GeV if (y1.pT() < 40.*GeV) vetoEvent; if (y2.pT() < 30.*GeV) vetoEvent; // Require the two photons to be separated (dR>0.4) if (deltaR(y1,y2) < 0.4) vetoEvent; const FourMomentum yy = y1 + y2; const double Myy = yy.mass(); const double pTyy = yy.pT(); const double dPhiyy = mapAngle0ToPi(y1.phi() - y2.phi()); // phi* const double costhetastar_ = fabs(tanh(( y1.eta() - y2.eta() ) / 2.)); const double sinthetastar_ = sqrt(1. - pow(costhetastar_, 2)); const double phistar = tan(0.5 * (PI - dPhiyy)) * sinthetastar_; // a_t const Vector3 t_hat(y1.x()-y2.x(), y1.y()-y2.y(), 0.); const double factor = t_hat.mod(); const Vector3 t_hatx(t_hat.x()/factor, t_hat.y()/factor, t_hat.z()/factor); const Vector3 At(y1.x()+y2.x(), y1.y()+y2.y(), 0.); // Compute a_t transverse component with respect to t_hat const double at = At.cross(t_hatx).mod(); // Fill histograms _h_M->fill(Myy); _h_pT->fill(pTyy); _h_dPhi->fill(dPhiyy); _h_costh->fill(costhetastar_); _h_phistar->fill(phistar); _h_at->fill(at); } // Normalise histograms etc., after the run void finalize() { - const double sf = crossSection() / (femtobarn * sumOfWeights()); + const double sf = crossSection()/femtobarn / sumOfWeights(); scale({_h_M, _h_pT, _h_dPhi, _h_costh, _h_phistar, _h_at}, sf); } private: Histo1DPtr _h_M, _h_pT, _h_dPhi, _h_costh, _h_phistar, _h_at; fastjet::AreaDefinition* _area_def; const vector ETA_BINS = {0.0, 1.5, 3.0}; vector _ptDensity, _sigma, _Njets; }; - // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1591327); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1604029.cc b/analyses/pluginATLAS/ATLAS_2017_I1604029.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1604029.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1604029.cc @@ -1,149 +1,153 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/UnstableFinalState.hh" namespace Rivet { - ///@brief: ttbar + gamma at 8 TeV + + /// @brief ttbar + gamma at 8 TeV class ATLAS_2017_I1604029 : public Analysis { public: // Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1604029); // Book histograms and initialise projections before the run void init() { const FinalState fs; // signal photons PromptFinalState prompt_ph(Cuts::abspid == PID::PHOTON && Cuts::pT > 15*GeV && Cuts::abseta < 2.37); declare(prompt_ph, "photons"); - // bare leptons + // bare leptons Cut base_cuts = (Cuts::abseta < 2.7) && (Cuts::pT > 10*GeV); IdentifiedFinalState bare_leps(base_cuts); bare_leps.acceptIdPair(PID::MUON); bare_leps.acceptIdPair(PID::ELECTRON); declare(bare_leps, "bare_leptons"); // dressed leptons Cut dressed_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV); PromptFinalState prompt_mu(base_cuts && Cuts::abspid == PID::MUON); PromptFinalState prompt_el(base_cuts && Cuts::abspid == PID::ELECTRON); IdentifiedFinalState all_photons(fs, PID::PHOTON); DressedLeptons elecs(all_photons, prompt_el, 0.1, dressed_cuts); declare(elecs, "elecs"); DressedLeptons muons(all_photons, prompt_mu, 0.1, dressed_cuts); declare(muons, "muons"); // auxiliary projections for 'single-lepton ttbar filter' PromptFinalState prompt_lep(Cuts::abspid == PID::MUON || Cuts::abspid == PID::ELECTRON); declare(prompt_lep, "prompt_leps"); declare(UnstableFinalState(), "ufs"); // jets - FastJets jets(fs, FastJets::ANTIKT, 0.4, JetAlg::NO_MUONS, JetAlg::NO_INVISIBLES); + FastJets jets(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); declare(jets, "jets"); + // BOOK HISTOGRAMS - _h["pt"] = bookHisto1D(2,1,1); - _h["eta"] = bookHisto1D(3,1,1); + book(_h["pt"], 2,1,1); + book(_h["eta"], 3,1,1); } + // Perform the per-event analysis void analyze(const Event& event) { - - const double weight = event.weight(); // analysis extrapolated to 1-lepton-plus-jets channel, where "lepton" cannot be a tau - // (i.e. contribution from dileptonic ttbar where one of the leptons is outside + // (i.e. contribution from dileptonic ttbar where one of the leptons is outside // the detector acceptance has been subtracted as a background) if (applyProjection(event, "prompt_leps").particles().size() != 1) vetoEvent; for (const auto& p : apply(event, "ufs").particles()) { if (p.fromPromptTau()) vetoEvent; } // photon selection Particles photons = applyProjection(event, "photons").particlesByPt(); Particles bare_leps = apply(event, "bare_leptons").particles(); for (const Particle& lep : bare_leps) ifilter_discard(photons, deltaRLess(lep, 0.1)); if (photons.size() != 1) vetoEvent; const Particle& photon = photons[0]; // jet selection Jets jets = apply(event, "jets").jetsByPt(Cuts::abseta < 2.5 && Cuts::pT > 25*GeV); // lepton selection const vector& elecs = apply(event, "elecs").dressedLeptons(); const vector& all_muons = apply(event, "muons").dressedLeptons(); // jet photon/electron overlap removal for (const DressedLepton& e : elecs) ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY)); for (const Particle& ph : photons) ifilter_discard(jets, deltaRLess(ph, 0.1, RAPIDITY)); if (jets.size() < 4) vetoEvent; // photon-jet minimum deltaR double mindR_phjet = 999.; for (Jet jet : jets) { const double dR_phjet = deltaR(photon, jet); if (dR_phjet < mindR_phjet) mindR_phjet = dR_phjet; } if (mindR_phjet < 0.5) vetoEvent; // muon jet overlap removal vector muons; - foreach (DressedLepton mu, all_muons) { + for (const DressedLepton& mu : all_muons) { bool overlaps = false; - foreach (Jet jet, jets) { + for (const Jet& jet : jets) { if (deltaR(mu, jet) < 0.4) { overlaps = true; break; } } if (overlaps) continue; muons.push_back(mu); } - // one electron XOR one muon + // one electron XOR one muon bool isEl = elecs.size() == 1 && muons.size() == 0; bool isMu = muons.size() == 1 && elecs.size() == 0; if (!isEl && !isMu) vetoEvent; // photon-lepton deltaR double mindR_phlep = deltaR(photon, isEl? elecs[0] : muons[0]); if (mindR_phlep < 0.7) vetoEvent; // b-tagging Jets bjets; - foreach (Jet jet, jets) { - if (jet.bTagged(Cuts::pT > 5*GeV)) bjets +=jet; + for (const Jet& jet : jets) { + if (jet.bTagged(Cuts::pT > 5*GeV)) bjets +=jet; } if (bjets.empty()) vetoEvent; - _h["pt"]->fill(photon.pT()/GeV, weight); - _h["eta"]->fill(photon.abseta(), weight); + _h["pt"]->fill(photon.pT()/GeV); + _h["eta"]->fill(photon.abseta()); } + // Normalise histograms etc., after the run void finalize() { const double normto(crossSection() / femtobarn / sumOfWeights()); for (auto &hist : _h) { scale(hist.second, normto); } } + private: map _h; + }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1604029); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1626105.cc b/analyses/pluginATLAS/ATLAS_2017_I1626105.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1626105.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1626105.cc @@ -1,130 +1,136 @@ +// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { + /// @brief lepton differential ttbar analysis at 8 TeV class ATLAS_2017_I1626105 : public Analysis { public: - + DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1626105); - + void init() { Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV; // All final state particles const FinalState fs; // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons PromptFinalState prompt_el(Cuts::abspid == PID::ELECTRON, true); DressedLeptons elecs(photons, prompt_el, 0.1, (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV)); DressedLeptons veto_elecs(photons, prompt_el, 0.1, eta_full, false); declare(elecs, "elecs"); // Projection to find the muons PromptFinalState prompt_mu(Cuts::abspid == PID::MUON, true); DressedLeptons muons(photons, prompt_mu, 0.1, (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV)); DressedLeptons veto_muons(photons, prompt_mu, 0.1, eta_full, false); declare(muons, "muons"); // Jet clustering. VetoedFinalState vfs; vfs.addVetoOnThisFinalState(veto_elecs); vfs.addVetoOnThisFinalState(veto_muons); FastJets jets(vfs, FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); // Book histograms bookHistos("lep_pt", 1); bookHistos("lep_eta", 3); bookHistos("dilep_pt", 5); bookHistos("dilep_mass", 7); bookHistos("dilep_rap", 9); bookHistos("dilep_dphi", 11); bookHistos("dilep_sumpt", 13); bookHistos("dilep_sumE", 15); } + void analyze(const Event& event) { vector elecs = sortByPt(apply(event, "elecs").dressedLeptons()); - vector muons = sortByPt(apply(event, "muons").dressedLeptons()); + vector muons = sortByPt(apply(event, "muons").dressedLeptons()); Jets jets = apply(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); - + // Check overlap of jets/leptons. for (const Jet& jet : jets) { for (const DressedLepton& el : elecs) { if (deltaR(jet, el) < 0.4) delete el; } for (const DressedLepton& mu : muons) { if (deltaR(jet, mu) < 0.4) delete mu; } } if (elecs.empty() || muons.empty()) vetoEvent; - - if (elecs[0].charge() == muons[0].charge()) vetoEvent; - + + if (elecs[0].charge() == muons[0].charge()) vetoEvent; + FourMomentum el = elecs[0].momentum(); FourMomentum mu = muons[0].momentum(); FourMomentum ll = elecs[0].momentum() + muons[0].momentum(); - + // Fill histograms - const double weight = event.weight(); - fillHistos("lep_pt", el.pT()/GeV, weight); - fillHistos("lep_pt", mu.pT()/GeV, weight); - fillHistos("lep_eta", el.abseta(), weight); - fillHistos("lep_eta", mu.abseta(), weight); - fillHistos("dilep_pt", ll.pT()/GeV, weight); - fillHistos("dilep_mass", ll.mass()/GeV, weight); - fillHistos("dilep_rap", ll.absrap(), weight); - fillHistos("dilep_dphi", deltaPhi(el, mu), weight); - fillHistos("dilep_sumpt", (el.pT() + mu.pT())/GeV, weight); - fillHistos("dilep_sumE", (el.E() + mu.E())/GeV, weight); + fillHistos("lep_pt", el.pT()/GeV); + fillHistos("lep_pt", mu.pT()/GeV); + fillHistos("lep_eta", el.abseta()); + fillHistos("lep_eta", mu.abseta()); + fillHistos("dilep_pt", ll.pT()/GeV); + fillHistos("dilep_mass", ll.mass()/GeV); + fillHistos("dilep_rap", ll.absrap()); + fillHistos("dilep_dphi", deltaPhi(el, mu)); + fillHistos("dilep_sumpt", (el.pT() + mu.pT())/GeV); + fillHistos("dilep_sumE", (el.E() + mu.E())/GeV); } + void finalize() { // Normalize to cross-section const double sf = crossSection()/femtobarn/sumOfWeights(); for (auto& hist : _h) { const double norm = 1.0 / hist.second->integral(); - // add overflow to last bin + // add overflow to last bin double overflow = hist.second->overflow().effNumEntries(); hist.second->fillBin(hist.second->numBins() - 1, overflow); // histogram normalisation if (hist.first.find("norm") != string::npos) scale(hist.second, norm); else scale(hist.second, sf); } } + private: /// @name Histogram helper functions //@{ void bookHistos(const std::string name, unsigned int index) { - _h[name] = bookHisto1D(index, 1, 1); - _h["norm_" + name] = bookHisto1D(index + 1, 1, 1); + book(_h[name], index, 1, 1); + book(_h["norm_" + name], index + 1, 1, 1); } - void fillHistos(const std::string name, double value, double weight) { - _h[name]->fill(value, weight); - _h["norm_" + name]->fill(value, weight); + void fillHistos(const std::string name, double value) { + _h[name]->fill(value); + _h["norm_" + name]->fill(value); } map _h; //@} }; + // Declare the class as a hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1626105); + } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1644367.cc b/analyses/pluginATLAS/ATLAS_2017_I1644367.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1644367.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1644367.cc @@ -1,185 +1,187 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { + /// @brief Isolated triphotons at 8 TeV class ATLAS_2017_I1644367 : public Analysis { public: // Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1644367); - + // Book histograms and initialise projections before the run void init() { - + const FinalState fs; - declare(fs, "FS"); + declare(fs, "FS"); FastJets fj(fs, FastJets::KT, 0.5); - fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec())); - declare(fj, "KtJetsD05"); + fj.useJetArea(new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec())); + declare(fj, "KtJetsD05"); - IdentifiedFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 15*GeV); + IdentifiedFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 15*GeV); declare(photonfs, "Photon"); - - // Histograms - _h["etg1"] = bookHisto1D( 1, 1, 1); - _h["etg2"] = bookHisto1D( 2, 1, 1); - _h["etg3"] = bookHisto1D( 3, 1, 1); - _h["dphig1g2"] = bookHisto1D( 4, 1, 1); - _h["dphig1g3"] = bookHisto1D( 5, 1, 1); - _h["dphig2g3"] = bookHisto1D( 6, 1, 1); - _h["detag1g2"] = bookHisto1D( 7, 1, 1); - _h["detag1g3"] = bookHisto1D( 8, 1, 1); - _h["detag2g3"] = bookHisto1D( 9, 1, 1); - _h["mg1g2"] = bookHisto1D(10, 1, 1); - _h["mg1g3"] = bookHisto1D(11, 1, 1); - _h["mg2g3"] = bookHisto1D(12, 1, 1); - _h["mg1g2g3"] = bookHisto1D(13, 1, 1); + + // Histograms + book(_h["etg1"] , 1, 1, 1); + book(_h["etg2"] , 2, 1, 1); + book(_h["etg3"] , 3, 1, 1); + book(_h["dphig1g2"], 4, 1, 1); + book(_h["dphig1g3"], 5, 1, 1); + book(_h["dphig2g3"], 6, 1, 1); + book(_h["detag1g2"], 7, 1, 1); + book(_h["detag1g3"], 8, 1, 1); + book(_h["detag2g3"], 9, 1, 1); + book(_h["mg1g2"] , 10, 1, 1); + book(_h["mg1g3"] , 11, 1, 1); + book(_h["mg2g3"] , 12, 1, 1); + book(_h["mg1g2g3"] , 13, 1, 1); } - + + // Perform the per-event analysis void analyze(const Event& event) { - const double weight = event.weight(); - // Require at least 2 photons in final state - const Particles photons = apply(event, "Photon").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.5); - if (photons.size() < 3) vetoEvent; + // Require at least 2 photons in final state + const Particles photons = apply(event, "Photon").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.5); + if (photons.size() < 3) vetoEvent; - // Get jets, and corresponding jet areas - vector > ptDensities(ETA_BINS.size()-1); + // Get jets, and corresponding jet areas + vector > ptDensities(ETA_BINS.size()-1); FastJets fastjets = apply(event, "KtJetsD05"); - const auto clust_seq_area = fastjets.clusterSeqArea(); - for (const Jet& jet : fastjets.jets()) { - const double area = clust_seq_area->area(jet); - if (area < 1e-3) continue; - const int ieta = binIndex(jet.abseta(), ETA_BINS); - if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area); - } + const auto clust_seq_area = fastjets.clusterSeqArea(); + for (const Jet& jet : fastjets.jets()) { + const double area = clust_seq_area->area(jet); + if (area < 1e-3) continue; + const int ieta = binIndex(jet.abseta(), ETA_BINS); + if (ieta != -1) ptDensities[ieta].push_back(jet.pT()/area); + } - // Compute median jet properties over the jets in the event - vector ptDensity; - for (size_t b = 0; b < ETA_BINS.size()-1; ++b) { - double median = 0.0; - if (ptDensities[b].size() > 0) { - std::sort(ptDensities[b].begin(), ptDensities[b].end()); - int nDens = ptDensities[b].size(); - median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2]; + // Compute median jet properties over the jets in the event + vector ptDensity; + for (size_t b = 0; b < ETA_BINS.size()-1; ++b) { + double median = 0.0; + if (ptDensities[b].size() > 0) { + std::sort(ptDensities[b].begin(), ptDensities[b].end()); + int nDens = ptDensities[b].size(); + median = (nDens % 2 == 0) ? (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2 : ptDensities[b][(nDens-1)/2]; } ptDensity.push_back(median); } // Loop over photons and fill vector of isolated ones Particles isolated_photons; for (const Particle& photon : photons) { - if (photon.fromDecay()) continue; - - // Remove photons in ECAL crack region - const double eta_P = photon.eta(); - const double phi_P = photon.phi(); - - // Compute isolation via particles within an R=0.4 cone of the photon - const Particles fs = apply(event, "FS").particles(); - FourMomentum mom_in_EtCone; - for (const Particle& p : fs) { - // Reject if not in cone - if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue; - // Reject if in the 5x7 cell central core + if (!photon.isPrompt()) continue; + + // Remove photons in ECAL crack region + const double eta_P = photon.eta(); + const double phi_P = photon.phi(); + + // Compute isolation via particles within an R=0.4 cone of the photon + const Particles fs = apply(event, "FS").particles(); + FourMomentum mom_in_EtCone; + for (const Particle& p : fs) { + // Reject if not in cone + if (deltaR(photon.momentum(), p.momentum()) > 0.4) continue; + // Reject if in the 5x7 cell central core if (fabs(eta_P - p.eta()) < 0.025 * 5 * 0.5 && fabs(phi_P - p.phi()) < PI/128. * 7 * 0.5) continue; // Sum momentum mom_in_EtCone += p.momentum(); } // Now figure out the correction (area*density) const double EtCone_area = M_PI*sqr(0.4) - (7*.025)*(5*M_PI/128.); // cone area - central core rectangle const double correction = ptDensity[binIndex(fabs(eta_P), ETA_BINS)] * EtCone_area; // Discard the photon if there is more than 11 GeV of cone activity // NOTE: Shouldn't need to subtract photon itself (it's in the central core) if (mom_in_EtCone.Et() - correction > 10*GeV) continue; // Add isolated photon to list isolated_photons.push_back(photon); - }///loop over photons - - // Require at least two isolated photons - if (isolated_photons.size() < 3) vetoEvent; - - // Select leading pT pair - sortByPt(isolated_photons); - const FourMomentum y1 = isolated_photons[0]; - const FourMomentum y2 = isolated_photons[1]; - const FourMomentum y3 = isolated_photons[2]; - - // Leading photon should have pT > 40 GeV, subleading > 30 GeV - if (y1.pT() < 27*GeV) vetoEvent; - if (y2.pT() < 22*GeV) vetoEvent; - if (y3.pT() < 15*GeV) vetoEvent; - + }///loop over photons + + // Require at least two isolated photons + if (isolated_photons.size() < 3) vetoEvent; + + // Select leading pT pair + sortByPt(isolated_photons); + const FourMomentum y1 = isolated_photons[0]; + const FourMomentum y2 = isolated_photons[1]; + const FourMomentum y3 = isolated_photons[2]; + + // Leading photon should have pT > 40 GeV, subleading > 30 GeV + if (y1.pT() < 27*GeV) vetoEvent; + if (y2.pT() < 22*GeV) vetoEvent; + if (y3.pT() < 15*GeV) vetoEvent; + // Require the two photons to be separated (dR>0.4) - if (deltaR(y1,y2) < 0.45) vetoEvent; - if (deltaR(y1,y3) < 0.45) vetoEvent; - if (deltaR(y2,y3) < 0.45) vetoEvent; - - - const FourMomentum yyy = y1 + y2 + y3; - const FourMomentum y1y2 = y1 + y2; - const FourMomentum y1y3 = y1 + y3; - const FourMomentum y2y3 = y2 + y3; + if (deltaR(y1,y2) < 0.45) vetoEvent; + if (deltaR(y1,y3) < 0.45) vetoEvent; + if (deltaR(y2,y3) < 0.45) vetoEvent; - const double Myyy = yyy.mass() / GeV; - const double dPhiy1y2 = mapAngle0ToPi(deltaPhi(y1, y2)); - const double dPhiy1y3 = mapAngle0ToPi(deltaPhi(y1, y3)); - const double dPhiy2y3 = mapAngle0ToPi(deltaPhi(y2, y3)); + const FourMomentum yyy = y1 + y2 + y3; + const FourMomentum y1y2 = y1 + y2; + const FourMomentum y1y3 = y1 + y3; + const FourMomentum y2y3 = y2 + y3; - const double dEtay1y2 = fabs(y1.eta() - y2.eta()); - const double dEtay1y3 = fabs(y1.eta() - y3.eta()); - const double dEtay2y3 = fabs(y2.eta() - y3.eta()); + const double Myyy = yyy.mass() / GeV; + + const double dPhiy1y2 = mapAngle0ToPi(deltaPhi(y1, y2)); + const double dPhiy1y3 = mapAngle0ToPi(deltaPhi(y1, y3)); + const double dPhiy2y3 = mapAngle0ToPi(deltaPhi(y2, y3)); + + const double dEtay1y2 = fabs(y1.eta() - y2.eta()); + const double dEtay1y3 = fabs(y1.eta() - y3.eta()); + const double dEtay2y3 = fabs(y2.eta() - y3.eta()); if(Myyy < 50.) vetoEvent; // Fill histograms - _h["etg1"]->fill(y1.pT() / GeV, weight); - _h["etg2"]->fill(y2.pT() / GeV, weight); - _h["etg3"]->fill(y3.pT() / GeV, weight); + _h["etg1"]->fill(y1.pT() / GeV); + _h["etg2"]->fill(y2.pT() / GeV); + _h["etg3"]->fill(y3.pT() / GeV); - _h["dphig1g2"]->fill(dPhiy1y2); - _h["dphig1g3"]->fill(dPhiy1y3); - _h["dphig2g3"]->fill(dPhiy2y3); + _h["dphig1g2"]->fill(dPhiy1y2); + _h["dphig1g3"]->fill(dPhiy1y3); + _h["dphig2g3"]->fill(dPhiy2y3); - _h["detag1g2"]->fill(dEtay1y2); - _h["detag1g3"]->fill(dEtay1y3); - _h["detag2g3"]->fill(dEtay2y3); + _h["detag1g2"]->fill(dEtay1y2); + _h["detag1g3"]->fill(dEtay1y3); + _h["detag2g3"]->fill(dEtay2y3); - _h["mg1g2"]->fill(y1y2.mass() / GeV, weight); - _h["mg1g3"]->fill(y1y3.mass() / GeV, weight); - _h["mg2g3"]->fill(y2y3.mass() / GeV, weight); + _h["mg1g2"]->fill(y1y2.mass() / GeV); + _h["mg1g3"]->fill(y1y3.mass() / GeV); + _h["mg2g3"]->fill(y2y3.mass() / GeV); - _h["mg1g2g3"]->fill(Myyy, weight); - } - - - // Normalise histograms etc., after the run - void finalize() { - const double sf = crossSection() / (femtobarn * sumOfWeights()); + _h["mg1g2g3"]->fill(Myyy); + } + + + // Normalise histograms etc., after the run + void finalize() { + const double sf = crossSection() / (femtobarn * sumOfWeights()); for (auto &hist : _h) { scale(hist.second, sf); } - } - - - private: - - map _h; - const vector ETA_BINS = { 0.0, 1.5, 3.0 }; - - }; - - // The hook for the plugin system - DECLARE_RIVET_PLUGIN(ATLAS_2017_I1644367); + } + + + private: + + map _h; + const vector ETA_BINS = { 0.0, 1.5, 3.0 }; + + }; + + + DECLARE_RIVET_PLUGIN(ATLAS_2017_I1644367); + } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1645627.cc b/analyses/pluginATLAS/ATLAS_2017_I1645627.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1645627.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1645627.cc @@ -1,154 +1,153 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief Isolated photon + jets at 13 TeV class ATLAS_2017_I1645627 : public Analysis { public: // Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1645627); // Book histograms and initialise projections before the run void init() { const FinalState fs; // calorimeter particles VisibleFinalState visFS(fs); VetoedFinalState calo_fs(visFS); calo_fs.addVetoPairId(PID::MUON); declare(calo_fs, "calo"); // Voronoi eta-phi tessellation with KT jets, for ambient energy density calculation - FastJets fj(fs, FastJets::KT, 0.5, JetAlg::NO_MUONS, JetAlg::NO_INVISIBLES); // E-scheme used by default; + FastJets fj(fs, FastJets::KT, 0.5, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); // E-scheme used by default; fj.useJetArea(new fastjet::AreaDefinition(fastjet::voronoi_area, fastjet::VoronoiAreaSpec(1.0))); declare(fj, "KtJetsD05"); // photon PromptFinalState photonfs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 2.37 && Cuts::pT > 125*GeV); declare(photonfs, "photons"); // Jets - FastJets jetpro(fs, FastJets::ANTIKT, 0.4, JetAlg::NO_MUONS, JetAlg::NO_INVISIBLES); + FastJets jetpro(fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); declare(jetpro, "Jets"); - - _h_photon_pt = bookHisto1D(1, 1, 1); - _h_jet_pt = bookHisto1D(2, 1, 1); - _h_phjet_dphi = bookHisto1D(3, 1, 1); - _h_phjet_mass = bookHisto1D(4, 1, 1); - _h_phjet_costheta = bookHisto1D(5, 1, 1); + + // Histograms + book(_h_photon_pt , 1, 1, 1); + book(_h_jet_pt , 2, 1, 1); + book(_h_phjet_dphi , 3, 1, 1); + book(_h_phjet_mass , 4, 1, 1); + book(_h_phjet_costheta, 5, 1, 1); } - size_t getEtaBin(double eta) const { - const double abseta = fabs(eta); - return binIndex(abseta, _eta_bins_areaoffset); + size_t getEtaBin(double eta) const { + return binIndex(fabs(eta), _eta_bins_areaoffset); } - + // Perform the per-event analysis void analyze(const Event& event) { - const double weight = event.weight(); - // Get the photon const Particles& photons = apply(event, "photons").particlesByPt(Cuts::abseta < 1.37 || Cuts::abseta > 1.56); if (photons.empty()) vetoEvent; const FourMomentum photon = photons[0].momentum(); // Get the jet Jets jets = apply(event, "Jets").jetsByPt(Cuts::pT > 100*GeV && Cuts::absrap < 2.37); ifilter_discard(jets, deltaRLess(photon, 0.8)); if (jets.empty()) vetoEvent; FourMomentum leadingJet = jets[0].momentum(); - + // Compute the jet pT densities vector< vector > ptDensities(_eta_bins_areaoffset.size()-1); FastJets fastjets = apply(event, "KtJetsD05"); const auto clust_seq_area = fastjets.clusterSeqArea(); for (const Jet& jet : fastjets.jets()) { - const double area = clust_seq_area->area(jet); // Implicit call to pseudojet(). - //const double area2 = (clust_seq_area->area_4vector(jet)).perp(); // Area definition used in egammaTruthParticles. + const double area = clust_seq_area->area(jet); // Implicit call to pseudojet(). + //const double area2 = (clust_seq_area->area_4vector(jet)).perp(); // Area definition used in egammaTruthParticles. if (area > 1e-3 && jet.abseta() < _eta_bins_areaoffset.back()) { ptDensities.at(getEtaBin(jet.abseta())) += jet.pT()/area; } } // Compute the median event energy density vector ptDensity; for (size_t b = 0; b < _eta_bins_areaoffset.size()-1; ++b) { ptDensity += ptDensities[b].empty() ? 0 : median(ptDensities[b]); } // Compute photon isolation with a standard ET cone FourMomentum mom_in_EtCone; const Particles calo_fs = apply(event, "calo").particles(); const double iso_dr = 0.4; for (const Particle& p : calo_fs) { // Check if it's in the cone of .4 if (sqrt(2.0*(cosh(p.eta()-photon.eta()) - cos(p.phi()-photon.phi()))) >= iso_dr) continue; // Increment sum mom_in_EtCone += p.momentum(); } - + // Remove the photon energy from the isolation mom_in_EtCone -= photon; // Figure out the correction (area*density) const double etcone_area = PI*iso_dr*iso_dr; const double correction = ptDensity[getEtaBin(photon.abseta())] * etcone_area; // Require photon to be isolated if ((mom_in_EtCone.Et()-correction) > (0.0042*photon.pT() + 10*GeV)) vetoEvent; - + // Fill histos const double photon_pt = photon.pT()/GeV; const double jet_pt = leadingJet.pT()/GeV; const double phjet_dphi = deltaPhi(photon, leadingJet); const double photon_eta = photon.eta(); const double jet_y = leadingJet.rapidity(); - _h_photon_pt->fill(photon_pt, weight); - _h_jet_pt->fill(jet_pt, weight); - _h_phjet_dphi->fill(phjet_dphi, weight); + _h_photon_pt->fill(photon_pt); + _h_jet_pt->fill(jet_pt); + _h_phjet_dphi->fill(phjet_dphi); double dy = fabs(jet_y-photon_eta); double phjet_costheta = tanh(dy/2.); double phjet_mass= (photon+leadingJet).mass()/GeV; if (phjet_mass <= 450.) vetoEvent; if (fabs(photon_eta + jet_y) >= 2.37) vetoEvent; if (phjet_costheta >= 0.83) vetoEvent; - _h_phjet_costheta->fill(phjet_costheta, weight); - _h_phjet_mass->fill(phjet_mass, weight); + _h_phjet_costheta->fill(phjet_costheta); + _h_phjet_mass->fill(phjet_mass); } /// Normalise histograms etc., after the run void finalize() { const double sf = crossSection() / sumOfWeights(); scale(_h_photon_pt, sf); scale(_h_jet_pt, sf); scale(_h_phjet_dphi, sf); scale(_h_phjet_mass, sf); scale(_h_phjet_costheta, sf); } private: Histo1DPtr _h_photon_pt, _h_jet_pt; Histo1DPtr _h_phjet_dphi, _h_phjet_mass, _h_phjet_costheta; const vector _eta_bins_areaoffset = {0.0, 1.5, 3.0, 4.0, 5.0}; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1645627); + } diff --git a/analyses/pluginATLAS/ATLAS_2018_I1646686.cc b/analyses/pluginATLAS/ATLAS_2018_I1646686.cc --- a/analyses/pluginATLAS/ATLAS_2018_I1646686.cc +++ b/analyses/pluginATLAS/ATLAS_2018_I1646686.cc @@ -1,299 +1,298 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/PartonicTops.hh" #include "Rivet/Math/LorentzTrans.hh" namespace Rivet { /// @brief All-hadronic ttbar at 13 TeV class ATLAS_2018_I1646686 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2018_I1646686); /// Book histograms and initialise projections before the run void init() { //histogram booking book(_h["inclusive"],2,1,1); bookHistograms("t_y", 1, true); bookHistograms("t_pt", 2, true); bookHistograms("t1_pt", 3); bookHistograms("t1_y", 4); bookHistograms("t2_pt", 5); bookHistograms("t2_y", 6); bookHistograms("tt_m", 7); bookHistograms("tt_pt", 8); bookHistograms("tt_y", 9); bookHistograms("tt_chi", 10); bookHistograms("tt_yboost", 11); bookHistograms("tt_pout", 12); bookHistograms("tt_dPhi", 13); bookHistograms("tt_Ht", 14); bookHistograms("tt_cosThStar", 15); // Projections Cut dressed_lep = (Cuts::abseta < 2.5) && (Cuts::pT >= 25*GeV); Cut eta_full = (Cuts::abseta < 5.0); // All final state particles FinalState fs(eta_full); // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); electrons.acceptTauDecays(true); DressedLeptons dressedelectrons(photons, electrons, 0.1, dressed_lep); declare(dressedelectrons, "elecs"); DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); DressedLeptons dressedmuons(photons, muons, 0.1, dressed_lep); declare(dressedmuons, "muons"); DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full); // Projection to find neutrinos for vetoing in jets IdentifiedFinalState nu_id(fs); nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); // Jet clustering. VetoedFinalState lvfs(fs); lvfs.addVetoOnThisFinalState(mu_id); lvfs.addVetoOnThisFinalState(nu_id); VetoedFinalState vfs; vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); FastJets jets(vfs, FastJets::ANTIKT, 0.4); jets.useInvisibles(); declare(jets, "jets"); FastJets ljets(lvfs, FastJets::ANTIKT, 1.0); ljets.useInvisibles(); declare(ljets, "ljets" ); PartonicTops partonTops; declare(partonTops, "partonicTops"); } void analyze(const Event& event) { // Parton-level top quarks const Particles partonicTops = apply( event, "partonicTops").particlesByPt(); FourMomentum top, tbar; bool foundT = false, foundTBar = false; for (const Particle& ptop : partonicTops) { const int pid = ptop.pid(); if (pid == PID::TQUARK) { top = ptop.momentum(); foundT = true; } else if (pid == -PID::TQUARK) { tbar = ptop.momentum(); foundTBar = true; } } FourMomentum t1_parton, t2_parton, ttbar_parton; if ( foundT && foundTBar ) { t1_parton = top.pT2() > tbar.pT2() ? top : tbar; t2_parton = top.pT2() > tbar.pT2() ? tbar : top; ttbar_parton = t1_parton + t2_parton; if ( t1_parton.pT() > 500*GeV && t2_parton.pT() > 350*GeV) { const double chi_parton = calcChi(t1_parton, t2_parton); const double cosThetaStar_parton = abs(calcCosThetaStar(t1_parton, t2_parton)); const double pout_parton = abs(calcPout(t1_parton, t2_parton)); const double dPhi_parton = deltaPhi(t1_parton, t2_parton); const int randomChoice = rand() % 2; const FourMomentum& randomTopParton = (randomChoice == 0) ? t1_parton : t2_parton; fillParton("t_pt", randomTopParton.pT()/GeV); fillParton("t_y", randomTopParton.absrap()); fillParton("t1_pt", t1_parton.pT()/GeV); fillParton("t1_y", t1_parton.absrap()); fillParton("t2_pt", t2_parton.pT()/GeV); fillParton("t2_y", t2_parton.absrap()); fillParton("tt_m", ttbar_parton.mass()/GeV); fillParton("tt_pt", ttbar_parton.pT()/GeV); fillParton("tt_Ht", (t1_parton.pT() + t2_parton.pT())/GeV); fillParton("tt_y", ttbar_parton.absrap()); fillParton("tt_yboost", 0.5 * abs(t1_parton.rapidity() + t2_parton.rapidity())); fillParton("tt_chi", chi_parton); fillParton("tt_cosThStar", cosThetaStar_parton); fillParton("tt_pout", pout_parton/GeV); fillParton("tt_dPhi", dPhi_parton); } } // Get and veto on dressed leptons const vector dressedElectrons = apply(event, "elecs").dressedLeptons(); const vector dressedMuons = apply(event, "muons").dressedLeptons(); if (!dressedElectrons.empty()) vetoEvent; if (!dressedMuons.empty()) vetoEvent; // Get jets const Jets& all_jets = apply( event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); const FastJets& ljets_fj = apply( event, "ljets"); const Jets all_ljets = ljets_fj.jetsByPt(); // Trim the large-R jets Jets trimmedJets; fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.2), fastjet::SelectorPtFractionMin(0.05)); for (const Jet& jet : all_ljets) trimmedJets += ljets_fj.trimJet(jet, trimmer); trimmedJets = sortByPt(trimmedJets); // Check large-R jets Jets ljets; vector b_tagged; for (const Jet& jet : trimmedJets) { if (jet.pT() < 250 * GeV) continue; if (jet.pT() > 3000 * GeV) continue; if (jet.mass() > jet.pT()) continue; if (jet.abseta() > 2.0 ) continue; ljets += jet; b_tagged += jet.bTagged(); } if (all_jets.size() < 2) vetoEvent; if (ljets.size() < 2) vetoEvent; // Identify top and anti top, compute some event variables const FourMomentum ttbar = ljets[0].momentum() + ljets[1].momentum(); const FourMomentum t1 = ljets[0].momentum(); const FourMomentum t2 = ljets[1].momentum(); const double chi = calcChi(t1, t2); const double cosThetaStar = abs(calcCosThetaStar(t1, t2)); const double pout = abs(calcPout(t1, t2)); const double dPhi = deltaPhi(t1, t2); if ( t2.pT() < 350*GeV) vetoEvent; if ( t1.pT() < 500*GeV) vetoEvent; // b-tagging for particle done on large-R jets if (!(b_tagged[0] && b_tagged[1])) vetoEvent; // Continues with signal region cuts if ( abs(t1.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent; if ( abs(t2.mass() - 172.5 * GeV) > 50*GeV ) vetoEvent; _h["inclusive"]->fill(0); fillHistograms("t1_pt", t1.pT()/GeV); fillHistograms("t1_y", t1.absrap()); fillHistograms("t2_pt", t2.pT()/GeV); fillHistograms("t2_y", t2.absrap()); fillHistograms("tt_m", ttbar.mass()/TeV); fillHistograms("tt_pt", ttbar.pT()/GeV); fillHistograms("tt_Ht", (t1.pT() + t2.pT())/GeV); fillHistograms("tt_y", ttbar.absrap()); fillHistograms("tt_yboost", 0.5 * abs(t1.rapidity() + t2.rapidity())); fillHistograms("tt_chi", chi); fillHistograms("tt_cosThStar", cosThetaStar); fillHistograms("tt_pout", pout/GeV); fillHistograms("tt_dPhi", dPhi); } void finalize() { // Normalize histograms const double sf = crossSection() / sumOfWeights(); for (auto &hist : _h) { scale(hist.second, sf); if ((hist.first.find("_norm") != string::npos) && hist.second->integral(false)>0) hist.second->normalize(1.0, false); } } double calcChi(const FourMomentum& t1, const FourMomentum& t2) { double ystar = 0.5 * (t1.rapidity()-t2.rapidity()); double chi = exp( 2 * abs(ystar)); return chi; } double calcCosThetaStar(const FourMomentum& t1, const FourMomentum& t2) { FourMomentum ttbar = t1 + t2; LorentzTransform centreOfMassTrans; ttbar.setX(0); ttbar.setY(0); centreOfMassTrans.setBetaVec( -ttbar.betaVec() ); FourMomentum t1_star = centreOfMassTrans.transform(t1); double cosThetaStar = t1_star.pz()/t1_star.p3().mod(); return cosThetaStar; } double calcPout(const FourMomentum& t1, const FourMomentum& t2) { Vector3 t1V = t1.p3(); Vector3 t2V = t2.p3(); Vector3 zUnit = Vector3(0., 0., 1.); Vector3 vPerp = zUnit.cross(t1V); double pout = vPerp.dot(t2V)/vPerp.mod(); return pout; } private: map _h; //some functions for booking, filling and scaling the histograms void fillHistograms(std::string name, double value) { _h[name]->fill(value); _h[name + "_norm"]->fill(value); } void fillParton(std::string name, double value) { _h[name + "_parton"]->fill(value); _h[name + "_parton_norm"]->fill(value); } void bookHistograms(const std::string name, unsigned int index, bool onlyParton = false) { if (!onlyParton) { book(_h[name], index, 1, 1 ); book(_h[name + "_norm"], index + 13, 1, 1 ); } book(_h[name + "_parton"], index + 82, 1, 1 ); book(_h[name + "_parton_norm"], index + 97, 1, 1 ); } }; DECLARE_RIVET_PLUGIN(ATLAS_2018_I1646686); - } diff --git a/analyses/pluginCMS/CMS_2014_I1303894.cc b/analyses/pluginCMS/CMS_2014_I1303894.cc --- a/analyses/pluginCMS/CMS_2014_I1303894.cc +++ b/analyses/pluginCMS/CMS_2014_I1303894.cc @@ -1,236 +1,235 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/InvMassFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/DressedLeptons.hh" namespace Rivet { /// @brief Differential cross-section of W bosons + jets in pp collisions at sqrt(s)=7 TeV /// @author Darin Baumgartel (darinb@cern.ch) /// /// Based on Rivet analysis originally created by Anil Singh (anil@cern.ch), Lovedeep Saini (lovedeep@cern.ch) class CMS_2014_I1303894 : public Analysis { public: /// Constructor CMS_2014_I1303894() : Analysis("CMS_2014_I1303894") { } // Book histograms and initialise projections before the run void init() { // Prompt leptons only, no test on nu flavour. // Projections const FinalState fs; declare(fs, "FS"); MissingMomentum missing(fs); declare(missing, "MET"); PromptFinalState pfs(fs); IdentifiedFinalState bareMuons(pfs); bareMuons.acceptIdPair(PID::MUON); DressedLeptons muonClusters(fs, bareMuons, -1); //, Cuts::open(), false, false); declare(muonClusters, "muonClusters"); IdentifiedFinalState neutrinos(pfs); neutrinos.acceptIdPair(PID::NU_MU); declare(neutrinos, "neutrinos"); VetoedFinalState jetFS(fs); jetFS.addVetoOnThisFinalState(muonClusters); jetFS.addVetoOnThisFinalState(neutrinos); jetFS.vetoNeutrinos(); - FastJets JetProjection(jetFS, FastJets::ANTIKT, 0.5); - JetProjection.useInvisibles(false); - declare(JetProjection, "Jets"); + FastJets jetprojection(jetFS, FastJets::ANTIKT, 0.5); + declare(jetprojection, "Jets"); // Histograms book(_histDPhiMuJet1 ,1,1,1); book(_histDPhiMuJet2 ,2,1,1); book(_histDPhiMuJet3 ,3,1,1); book(_histDPhiMuJet4 ,4,1,1); book(_histEtaJet1 ,5,1,1); book(_histEtaJet2 ,6,1,1); book(_histEtaJet3 ,7,1,1); book(_histEtaJet4 ,8,1,1); book(_histHT1JetInc ,9,1,1); book(_histHT2JetInc ,10,1,1); book(_histHT3JetInc ,11,1,1); book(_histHT4JetInc ,12,1,1); book(_histJet30MultExc ,13,1,1); book(_histJet30MultInc ,14,1,1); book(_histPtJet1 ,15,1,1); book(_histPtJet2 ,16,1,1); book(_histPtJet3 ,17,1,1); book(_histPtJet4 ,18,1,1); // Counters book(_n_1jet, "n_1jet"); book(_n_2jet, "n_2jet"); book(_n_3jet, "n_3jet"); book(_n_4jet, "n_4jet"); book(_n_inclusivebinsummation, "n_inclusivebinsummation"); } void analyze(const Event& event) { // Get the dressed muon const DressedLeptons& muonClusters = apply(event, "muonClusters"); int nmu = muonClusters.dressedLeptons().size(); if (nmu < 1) vetoEvent; DressedLepton dressedmuon = muonClusters.dressedLeptons()[0]; if (dressedmuon.momentum().abseta() > 2.1) vetoEvent; if (dressedmuon.momentum().pT() < 25.0*GeV) vetoEvent; // Get the muon neutrino //const Particles& neutrinos = apply(event, "neutrinos").particlesByPt(); // Check that the muon and neutrino are not decay products of tau if (dressedmuon.constituentLepton().hasAncestor( PID::TAU)) vetoEvent; if (dressedmuon.constituentLepton().hasAncestor(-PID::TAU)) vetoEvent; // Get the missing momentum const MissingMomentum& met = apply(event, "MET"); const double ptmet = met.visibleMomentum().pT(); const double phimet = (-met.visibleMomentum()).phi(); // Calculate MET and MT(mu,MET), and remove events with MT < 50 GeV const double ptmuon = dressedmuon.pT(); const double phimuon = dressedmuon.phi(); const double mt_mumet = sqrt(2*ptmuon*ptmet*(1.0 - cos(phimet-phimuon))); // Remove events in MT < 50 region if (mt_mumet < 50*GeV) vetoEvent; // Loop over jets and fill pt/eta/phi quantities in vectors const Jets& jets_filtered = apply(event, "Jets").jetsByPt(0.0*GeV); vector finaljet_pT_list, finaljet_eta_list, finaljet_phi_list; double htjets = 0.0; for (size_t ii = 0; ii < jets_filtered.size(); ++ii) { // Jet pT/eta/phi double jet_pt = jets_filtered[ii].pT(); double jet_eta = jets_filtered[ii].eta(); double jet_phi = jets_filtered[ii].phi(); // Kinemetic cuts for jet acceptance if (fabs(jet_eta) > 2.4) continue; if (jet_pt < 30.0*GeV) continue; if (deltaR(dressedmuon, jets_filtered[ii]) < 0.5) continue; // Add jet to jet list and increases the HT variable finaljet_pT_list.push_back(jet_pt); finaljet_eta_list.push_back(jet_eta); finaljet_phi_list.push_back(jet_phi); htjets += fabs(jet_pt); } // Filling of histograms: // Fill as many jets as there are into the exclusive jet multiplicity if (!finaljet_pT_list.empty()) _histJet30MultExc->fill(finaljet_pT_list.size()); for (size_t ij = 0; ij < finaljet_pT_list.size(); ++ij) { _histJet30MultInc->fill(ij+1); _n_inclusivebinsummation->fill(); } if (finaljet_pT_list.size() >= 1) { _histPtJet1->fill(finaljet_pT_list[0]); _histEtaJet1->fill(fabs(finaljet_eta_list[0])); _histDPhiMuJet1->fill(deltaPhi(finaljet_phi_list[0], phimuon)); _histHT1JetInc->fill(htjets); _n_1jet->fill(); } if (finaljet_pT_list.size() >= 2) { _histPtJet2->fill(finaljet_pT_list[1]); _histEtaJet2->fill(fabs(finaljet_eta_list[1])); _histDPhiMuJet2->fill(deltaPhi(finaljet_phi_list[1], phimuon)); _histHT2JetInc->fill(htjets); _n_2jet->fill(); } if (finaljet_pT_list.size() >= 3) { _histPtJet3->fill(finaljet_pT_list[2]); _histEtaJet3->fill(fabs(finaljet_eta_list[2])); _histDPhiMuJet3->fill(deltaPhi(finaljet_phi_list[2], phimuon)); _histHT3JetInc->fill(htjets); _n_3jet->fill(); } if (finaljet_pT_list.size() >=4 ) { _histPtJet4->fill(finaljet_pT_list[3]); _histEtaJet4->fill(fabs(finaljet_eta_list[3])); _histDPhiMuJet4->fill(deltaPhi(finaljet_phi_list[3], phimuon)); _histHT4JetInc-> fill(htjets); _n_4jet->fill(); } } // Finalize the histograms. void finalize() { const double inclusive_cross_section = crossSection(); const double norm_1jet_histo = inclusive_cross_section*dbl(*_n_1jet)/sumOfWeights(); const double norm_2jet_histo = inclusive_cross_section*dbl(*_n_2jet)/sumOfWeights(); const double norm_3jet_histo = inclusive_cross_section*dbl(*_n_3jet)/sumOfWeights(); const double norm_4jet_histo = inclusive_cross_section*dbl(*_n_4jet)/sumOfWeights(); const double norm_incmultiplicity = inclusive_cross_section*dbl(*_n_inclusivebinsummation)/sumOfWeights(); normalize(_histJet30MultExc, norm_1jet_histo); normalize(_histJet30MultInc, norm_incmultiplicity); normalize(_histPtJet1, norm_1jet_histo); normalize(_histHT1JetInc, norm_1jet_histo); normalize(_histEtaJet1, norm_1jet_histo); normalize(_histDPhiMuJet1, norm_1jet_histo); normalize(_histPtJet2, norm_2jet_histo); normalize(_histHT2JetInc, norm_2jet_histo); normalize(_histEtaJet2, norm_2jet_histo); normalize(_histDPhiMuJet2, norm_2jet_histo); normalize(_histPtJet3, norm_3jet_histo); normalize(_histHT3JetInc, norm_3jet_histo); normalize(_histEtaJet3, norm_3jet_histo); normalize(_histDPhiMuJet3, norm_3jet_histo); normalize(_histPtJet4, norm_4jet_histo); normalize(_histHT4JetInc, norm_4jet_histo); normalize(_histEtaJet4, norm_4jet_histo); normalize(_histDPhiMuJet4, norm_4jet_histo); } private: Histo1DPtr _histJet30MultExc, _histJet30MultInc; Histo1DPtr _histPtJet1, _histPtJet2, _histPtJet3, _histPtJet4; Histo1DPtr _histEtaJet1, _histEtaJet2, _histEtaJet3, _histEtaJet4; Histo1DPtr _histDPhiMuJet1, _histDPhiMuJet2, _histDPhiMuJet3, _histDPhiMuJet4; Histo1DPtr _histHT1JetInc, _histHT2JetInc, _histHT3JetInc, _histHT4JetInc; CounterPtr _n_1jet, _n_2jet, _n_3jet, _n_4jet, _n_inclusivebinsummation; }; DECLARE_RIVET_PLUGIN(CMS_2014_I1303894); } diff --git a/analyses/pluginCMS/CMS_2015_I1346843.cc b/analyses/pluginCMS/CMS_2015_I1346843.cc --- a/analyses/pluginCMS/CMS_2015_I1346843.cc +++ b/analyses/pluginCMS/CMS_2015_I1346843.cc @@ -1,112 +1,113 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedLeptons.hh" #include "Rivet/Projections/NeutralFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" namespace Rivet { /// Differential cross-section of FSR photons in Z decays class CMS_2015_I1346843 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2015_I1346843); /// Book histograms and initialise projections before the run void init() { Cut c_photons = Cuts::pT >= 5.0*GeV && (Cuts::etaIn(-2.5, 1.4) || Cuts::etaIn(1.6, 2.5)); IdentifiedFinalState photons(c_photons); photons.acceptId(PID::PHOTON); declare(photons, "PHOTFS"); Cut c_muons = Cuts::pT > 9*GeV && Cuts::abseta < 2.4; IdentifiedFinalState muons(c_muons); muons.acceptIdPair(PID::MUON); declare(muons, "MUFS"); book(_hist_pho_et ,1, 1, 1); // photon transverse energy book(_hist_pho_et_wide ,1, 2, 1); // photon transverse energy (0.5 < dr < 3.0) book(_hist_pho_et_close ,1, 3, 1); // photon transverse energy (0.05 < dr < 0.5) book(_hist_pho_et_lqt ,1, 4, 1); // photon transverse energy (q_T < 10) book(_hist_pho_et_hqt ,1, 5, 1); // photon transverse energy (q_T > 50) book(_hist_pho_dr ,2, 1, 1); // delta_R book(_hist_pho_dr_lqt ,2, 2, 1); // delta_R (q_T < 10) book(_hist_pho_dr_hqt ,2, 3, 1); // delta_R (q_T > 50) } // Perform the per-event analysis void analyze(const Event& event) { const Particles muons = apply(event, "MUFS").particlesByPt(); if (muons.size() < 2) vetoEvent; if (muons[0].pT()/GeV < 31) vetoEvent; if (muons[0].charge()*muons[1].charge() > 0) vetoEvent; const double mZ = (muons[0].momentum() + muons[1].momentum()).mass(); if (!inRange(mZ, 30*GeV, 87*GeV)) vetoEvent; const Particles photons = apply(event, "PHOTFS").particlesByPt(); // We want the photon with the highest pT that does not come from a decay - for(const Particle& p : photons) { - if (p.fromDecay() || !p.isStable()) continue; + for (const Particle& p : photons) { + if (!p.isDirect()) continue; + if (!p.isStable()) continue; const double dR = std::min(deltaR(p, muons[0]), deltaR(p, muons[1]) ); if (!inRange(dR, 0.05, 3.0)) continue; // Calculate the three-body (mu,mu,gamma) transverse momentum const double qT = (muons[0].mom() + muons[1].mom() + p.mom()).pT(); // Fill the analysis histograms _hist_pho_et->fill(p.pT()/GeV, 1.0); _hist_pho_dr->fill(dR, 1.0); (dR <= 0.5 ? _hist_pho_et_close : _hist_pho_et_wide)->fill(p.pT()/GeV, 1.0); if (qT / GeV < 10.) { _hist_pho_et_lqt->fill(p.pT()/GeV, 1.0); _hist_pho_dr_lqt->fill(dR, 1.0); } if (qT / GeV > 50.) { _hist_pho_et_hqt->fill(p.pT()/GeV, 1.0); _hist_pho_dr_hqt->fill(dR, 1.0); } break; // Exit the loop since we found the highest pT lepton already } } /// Normalise histograms etc., after the run void finalize() { scale(_hist_pho_et, crossSection() / sumOfWeights()); scale(_hist_pho_et_wide, crossSection() / sumOfWeights()); scale(_hist_pho_et_close, crossSection() / sumOfWeights()); scale(_hist_pho_et_lqt, crossSection() / sumOfWeights()); scale(_hist_pho_et_hqt, crossSection() / sumOfWeights()); scale(_hist_pho_dr, crossSection() / sumOfWeights()); scale(_hist_pho_dr_lqt, crossSection() / sumOfWeights()); scale(_hist_pho_dr_hqt, crossSection() / sumOfWeights()); } private: Histo1DPtr _hist_pho_et; Histo1DPtr _hist_pho_et_wide, _hist_pho_et_close; Histo1DPtr _hist_pho_et_lqt, _hist_pho_et_hqt; Histo1DPtr _hist_pho_dr; Histo1DPtr _hist_pho_dr_lqt, _hist_pho_dr_hqt; }; DECLARE_RIVET_PLUGIN(CMS_2015_I1346843); } 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 declare(PartonicTops(PartonicTops::DecayMode::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::Muons::ALL, JetAlg::Invisibles::ALL); declare(fj, "Jets"); // Book histograms book(_hVis_nJet30_abs , 1, 1, 1); book(_hVis_nJet30 , 2, 1, 1); book(_hVis_nJet60_abs , 3, 1, 1); book(_hVis_nJet60 , 4, 1, 1); book(_hVis_nJet100_abs , 5, 1, 1); book(_hVis_nJet100 , 6, 1, 1); book(_hVis_addJet1Pt_abs , 7, 1, 1); book(_hVis_addJet1Pt , 8, 1, 1); book(_hVis_addJet1Eta_abs , 9, 1, 1); book(_hVis_addJet1Eta ,10, 1, 1); book(_hVis_addJet2Pt_abs ,11, 1, 1); book(_hVis_addJet2Pt ,12, 1, 1); book(_hVis_addJet2Eta_abs ,13, 1, 1); book(_hVis_addJet2Eta ,14, 1, 1); book(_hVis_addJJMass_abs ,15, 1, 1); book(_hVis_addJJMass ,16, 1, 1); book(_hVis_addJJDR_abs ,17, 1, 1); book(_hVis_addJJDR ,18, 1, 1); book(_hVis_addJJHT_abs ,19, 1, 1); book(_hVis_addJJHT ,20, 1, 1); book(_hFull_addJet1Pt_abs ,21, 1, 1); book(_hFull_addJet1Pt ,22, 1, 1); book(_hFull_addJet1Eta_abs ,23, 1, 1); book(_hFull_addJet1Eta ,24, 1, 1); book(_hFull_addJet2Pt_abs ,25, 1, 1); book(_hFull_addJet2Pt ,26, 1, 1); book(_hFull_addJet2Eta_abs ,27, 1, 1); book(_hFull_addJet2Eta ,28, 1, 1); book(_hFull_addJJMass_abs ,29, 1, 1); book(_hFull_addJJMass ,30, 1, 1); book(_hFull_addJJDR_abs ,31, 1, 1); book(_hFull_addJJDR ,32, 1, 1); book(_hFull_addJJHT_abs ,33, 1, 1); book(_hFull_addJJHT ,34, 1, 1); book(_hVis_addBJet1Pt_abs ,35, 1, 1); book(_hVis_addBJet1Pt ,36, 1, 1); book(_hVis_addBJet1Eta_abs ,37, 1, 1); book(_hVis_addBJet1Eta ,38, 1, 1); book(_hVis_addBJet2Pt_abs ,39, 1, 1); book(_hVis_addBJet2Pt ,40, 1, 1); book(_hVis_addBJet2Eta_abs ,41, 1, 1); book(_hVis_addBJet2Eta ,42, 1, 1); book(_hVis_addBBMass_abs ,43, 1, 1); book(_hVis_addBBMass ,44, 1, 1); book(_hVis_addBBDR_abs ,45, 1, 1); book(_hVis_addBBDR ,46, 1, 1); book(_hFull_addBJet1Pt_abs ,47, 1, 1); book(_hFull_addBJet1Pt ,48, 1, 1); book(_hFull_addBJet1Eta_abs ,49, 1, 1); book(_hFull_addBJet1Eta ,50, 1, 1); book(_hFull_addBJet2Pt_abs ,51, 1, 1); book(_hFull_addBJet2Pt ,52, 1, 1); book(_hFull_addBJet2Eta_abs ,53, 1, 1); book(_hFull_addBJet2Eta ,54, 1, 1); book(_hFull_addBBMass_abs ,55, 1, 1); book(_hFull_addBBMass ,56, 1, 1); book(_hFull_addBBDR_abs ,57, 1, 1); book(_hFull_addBBDR ,58, 1, 1); book(_h_gap_addJet1Pt ,59, 1, 1); book(_h_gap_addJet1Pt_eta0 ,60, 1, 1); book(_h_gap_addJet1Pt_eta1 ,61, 1, 1); book(_h_gap_addJet1Pt_eta2 ,62, 1, 1); book(_h_gap_addJet2Pt ,63, 1, 1); book(_h_gap_addJet2Pt_eta0 ,64, 1, 1); book(_h_gap_addJet2Pt_eta1 ,65, 1, 1); book(_h_gap_addJet2Pt_eta2 ,66, 1, 1); book(_h_gap_addJetHT ,67, 1, 1); book(_h_gap_addJetHT_eta0 ,68, 1, 1); book(_h_gap_addJetHT_eta1 ,69, 1, 1); book(_h_gap_addJetHT_eta2 ,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 auto isPromptChLepton = [](const Particle& p){return p.isPrompt() && isChargedLepton(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)); 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 = 1.0; 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, 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, 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/analyses/pluginMC/MC_DILEPTON.cc b/analyses/pluginMC/MC_DILEPTON.cc --- a/analyses/pluginMC/MC_DILEPTON.cc +++ b/analyses/pluginMC/MC_DILEPTON.cc @@ -1,134 +1,132 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/PromptFinalState.hh" namespace Rivet { /// @brief Study of prompt dilepton kinematics sensitive to boson polarizations class MC_DILEPTON : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(MC_DILEPTON); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections declare(PromptFinalState((Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON) && Cuts::abseta < 5 && Cuts::pT > 10*GeV), "Leptons"); // Book histograms - _h_pt_l1 = bookHisto1D("lep1_pt", logspace(40, 10, 400)); - _h_costheta_l1 = bookHisto1D("lep1_costheta", linspace(25, -1, 1)); - _h_ppara_l1 = bookHisto1D("lep1_ppara", linspace(40, -50, 350)); - _h_pperp_l1 = bookHisto1D("lep1_pperp", linspace(25, 0, 100)); + book(_h_pt_l1, "lep1_pt", logspace(40, 10, 400)); + book(_h_costheta_l1, "lep1_costheta", linspace(25, -1, 1)); + book(_h_ppara_l1, "lep1_ppara", linspace(40, -50, 350)); + book(_h_pperp_l1, "lep1_pperp", linspace(25, 0, 100)); // - _h_pt_l2 = bookHisto1D("lep2_pt", logspace(40, 10, 400)); - _h_costheta_l2 = bookHisto1D("lep2_costheta", linspace(25, -1, 1)); - _h_ppara_l2 = bookHisto1D("lep2_ppara", linspace(40, -50, 350)); - _h_pperp_l2 = bookHisto1D("lep2_pperp", linspace(25, 0, 100)); + book(_h_pt_l2, "lep2_pt", logspace(40, 10, 400)); + book(_h_costheta_l2, "lep2_costheta", linspace(25, -1, 1)); + book(_h_ppara_l2, "lep2_ppara", linspace(40, -50, 350)); + book(_h_pperp_l2, "lep2_pperp", linspace(25, 0, 100)); // - _h_costheta_com_l1 = bookHisto1D("com_costheta_l1", linspace(25, -1, 1)); - _h_costheta_com_l2 = bookHisto1D("com_costheta_l2", linspace(25, -1, 1)); - _h_ppara_com_l1 = bookHisto1D("com_ppara_l1", linspace(25, -50, 50)); - _h_ppara_com_l2 = bookHisto1D("com_ppara_l2", linspace(25, -50, 50)); + book(_h_costheta_com_l1, "com_costheta_l1", linspace(25, -1, 1)); + book(_h_costheta_com_l2, "com_costheta_l2", linspace(25, -1, 1)); + book(_h_ppara_com_l1, "com_ppara_l1", linspace(25, -50, 50)); + book(_h_ppara_com_l2, "com_ppara_l2", linspace(25, -50, 50)); // - _h_costheta_com = bookHisto1D("com_costheta", linspace(25, -1, 1)); - _h_ppara_com = bookHisto1D("com_ppara", linspace(25, -50, 50)); - _h_pperp_com = bookHisto1D("com_pperp", linspace(25, 0, 100)); + book(_h_costheta_com, "com_costheta", linspace(25, -1, 1)); + book(_h_ppara_com, "com_ppara", linspace(25, -50, 50)); + book(_h_pperp_com, "com_pperp", linspace(25, 0, 100)); } /// Perform the per-event analysis void analyze(const Event& event) { - const double weight = event.weight(); - const Particles& leptons = apply(event, "Leptons").particlesByPt(); if (leptons.size() != 2) vetoEvent; const Particle& l1 = leptons[0]; const Particle& l2 = leptons[1]; - _h_pt_l1->fill(l1.pT()/GeV, weight); - _h_pt_l2->fill(l2.pT()/GeV, weight); + _h_pt_l1->fill(l1.pT()/GeV); + _h_pt_l2->fill(l2.pT()/GeV); const FourMomentum pcom = l1.mom() + l2.mom(); const Vector3 betacom = pcom.betaVec(); const Vector3 unitboostvec = betacom.unit(); const LorentzTransform comboost = LorentzTransform::mkFrameTransformFromBeta(betacom); const double l1_costheta = cos(l1.p3().angle(unitboostvec)); const double l1_ppara = l2.p3().dot(unitboostvec); const double l1_pperp = l2.p3().cross(unitboostvec).mod(); - _h_costheta_l1->fill(l1_costheta, weight); - _h_ppara_l1->fill(l1_ppara, weight); - _h_pperp_l1->fill(l1_pperp, weight); + _h_costheta_l1->fill(l1_costheta); + _h_ppara_l1->fill(l1_ppara); + _h_pperp_l1->fill(l1_pperp); const double l2_costheta = cos(l2.p3().angle(unitboostvec)); const double l2_ppara = l2.p3().dot(unitboostvec); const double l2_pperp = l2.p3().cross(unitboostvec).mod(); - _h_costheta_l2->fill(l2_costheta, weight); - _h_ppara_l2->fill(l2_ppara, weight); - _h_pperp_l2->fill(l2_pperp, weight); + _h_costheta_l2->fill(l2_costheta); + _h_ppara_l2->fill(l2_ppara); + _h_pperp_l2->fill(l2_pperp); const FourMomentum p1com = comboost.transform(l1.mom()); const FourMomentum p2com = comboost.transform(l2.mom()); const double com_costheta1 = cos(p1com.p3().angle(unitboostvec)); const double com_costheta2 = cos(p2com.p3().angle(unitboostvec)); MSG_DEBUG("CoM cos(th)s: " << com_costheta1 << ", " << com_costheta2); // assert(com_costheta1 == com_costheta2 && "CoM cos(th)s differ"); const double com_ppara1 = p1com.p3().dot(unitboostvec); const double com_ppara2 = p2com.p3().dot(unitboostvec); MSG_DEBUG("CoM p_paras: " << com_ppara1 << ", " << com_ppara2); // assert(com_ppara1 == com_ppara2 && "CoM p_paras differ"); const double com_pperp1 = p1com.p3().cross(unitboostvec).mod(); const double com_pperp2 = p2com.p3().cross(unitboostvec).mod(); MSG_DEBUG("CoM p_pperps: " << com_pperp1 << ", " << com_pperp2); // assert(com_pperp1 == com_pperp2 && "CoM p_perps differ"); - _h_costheta_com_l1->fill(com_costheta1, weight); - _h_costheta_com_l2->fill(com_costheta2, weight); - _h_costheta_com->fill(com_costheta1, 0.5*weight); - _h_costheta_com->fill(com_costheta2, 0.5*weight); - _h_ppara_com_l1->fill(com_ppara1, weight); - _h_ppara_com_l2->fill(com_ppara2, weight); - _h_ppara_com->fill(com_ppara1, 0.5*weight); - _h_ppara_com->fill(com_ppara2, 0.5*weight); - _h_pperp_com->fill(com_pperp1, weight); + _h_costheta_com_l1->fill(com_costheta1); + _h_costheta_com_l2->fill(com_costheta2); + _h_costheta_com->fill(com_costheta1, 0.5); + _h_costheta_com->fill(com_costheta2, 0.5); + _h_ppara_com_l1->fill(com_ppara1); + _h_ppara_com_l2->fill(com_ppara2); + _h_ppara_com->fill(com_ppara1, 0.5); + _h_ppara_com->fill(com_ppara2, 0.5); + _h_pperp_com->fill(com_pperp1); } /// Normalise histograms etc., after the run void finalize() { normalize({_h_pt_l1, _h_pt_l2}); normalize({_h_ppara_com, _h_pperp_com, _h_costheta_com}); normalize({_h_ppara_com_l1, _h_ppara_com_l2, _h_costheta_com_l1, _h_costheta_com_l2}); normalize({_h_ppara_l1, _h_pperp_l1, _h_costheta_l1}); normalize({_h_ppara_l2, _h_pperp_l2, _h_costheta_l2}); } //@} /// @name Histograms //@{ Histo1DPtr _h_pt_l1, _h_pt_l2; Histo1DPtr _h_ppara_com, _h_pperp_com, _h_costheta_com; Histo1DPtr _h_ppara_com_l1, _h_ppara_com_l2, _h_costheta_com_l1, _h_costheta_com_l2; Histo1DPtr _h_ppara_l1, _h_pperp_l1, _h_costheta_l1; Histo1DPtr _h_ppara_l2, _h_pperp_l2, _h_costheta_l2; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(MC_DILEPTON); } diff --git a/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc b/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc --- a/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc +++ b/analyses/pluginPetra/JADE_OPAL_2000_S4300807.cc @@ -1,186 +1,177 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/FinalState.hh" namespace Rivet { /// @brief Jet rates in \f$ e^+ e^- \f$ at OPAL and JADE /// @author Frank Siegert class JADE_OPAL_2000_S4300807 : public Analysis { public: - /// @name Constructors etc. - //@{ - /// Constructor - JADE_OPAL_2000_S4300807() - : Analysis("JADE_OPAL_2000_S4300807") - { } - - //@} + DEFAULT_RIVET_ANALYSIS_CTOR(JADE_OPAL_2000_S4300807) /// @name Analysis methods //@{ void init() { // Projections const FinalState fs; declare(fs, "FS"); - FastJets jadeJets = FastJets(fs, FastJets::JADE, 0.7); - FastJets durhamJets = FastJets(fs, FastJets::DURHAM, 0.7); - jadeJets.useInvisibles(true); - durhamJets.useInvisibles(true); + FastJets jadeJets = FastJets(fs, FastJets::JADE, 0.7, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); + FastJets durhamJets = FastJets(fs, FastJets::DURHAM, 0.7, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jadeJets, "JadeJets"); declare(durhamJets, "DurhamJets"); // Histos int offset = 0; switch (int(sqrtS()/GeV + 0.5)) { case 35: offset = 7; break; case 44: offset = 8; break; case 91: offset = 9; break; case 133: offset = 10; break; case 161: offset = 11; break; case 172: offset = 12; break; case 183: offset = 13; break; case 189: offset = 14; break; default: break; } for (size_t i = 0; i < 5; ++i) { book(_h_R_Jade[i] ,offset, 1, i+1); book(_h_R_Durham[i] ,offset+9, 1, i+1); if (i < 4) book(_h_y_Durham[i] ,offset+17, 1, i+1); } } void analyze(const Event& e) { MSG_DEBUG("Num particles = " << apply(e, "FS").particles().size()); const FastJets& jadejet = apply(e, "JadeJets"); if (jadejet.clusterSeq()) { const double y_23 = jadejet.clusterSeq()->exclusive_ymerge_max(2); const double y_34 = jadejet.clusterSeq()->exclusive_ymerge_max(3); const double y_45 = jadejet.clusterSeq()->exclusive_ymerge_max(4); const double y_56 = jadejet.clusterSeq()->exclusive_ymerge_max(5); for (size_t i = 0; i < _h_R_Jade[0]->numBins(); ++i) { double ycut = _h_R_Jade[0]->bin(i).xMid(); double width = _h_R_Jade[0]->bin(i).xWidth(); if (y_23 < ycut) { _h_R_Jade[0]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Jade[1]->numBins(); ++i) { double ycut = _h_R_Jade[1]->bin(i).xMid(); double width = _h_R_Jade[1]->bin(i).xWidth(); if (y_34 < ycut && y_23 > ycut) { _h_R_Jade[1]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Jade[2]->numBins(); ++i) { double ycut = _h_R_Jade[2]->bin(i).xMid(); double width = _h_R_Jade[2]->bin(i).xWidth(); if (y_45 < ycut && y_34 > ycut) { _h_R_Jade[2]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Jade[3]->numBins(); ++i) { double ycut = _h_R_Jade[3]->bin(i).xMid(); double width = _h_R_Jade[3]->bin(i).xWidth(); if (y_56 < ycut && y_45 > ycut) { _h_R_Jade[3]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Jade[4]->numBins(); ++i) { double ycut = _h_R_Jade[4]->bin(i).xMid(); double width = _h_R_Jade[4]->bin(i).xWidth(); if (y_56 > ycut) { _h_R_Jade[4]->fillBin(i, width); } } } const FastJets& durjet = apply(e, "DurhamJets"); if (durjet.clusterSeq()) { const double y_23 = durjet.clusterSeq()->exclusive_ymerge_max(2); const double y_34 = durjet.clusterSeq()->exclusive_ymerge_max(3); const double y_45 = durjet.clusterSeq()->exclusive_ymerge_max(4); const double y_56 = durjet.clusterSeq()->exclusive_ymerge_max(5); _h_y_Durham[0]->fill(y_23); _h_y_Durham[1]->fill(y_34); _h_y_Durham[2]->fill(y_45); _h_y_Durham[3]->fill(y_56); for (size_t i = 0; i < _h_R_Durham[0]->numBins(); ++i) { double ycut = _h_R_Durham[0]->bin(i).xMid(); double width = _h_R_Durham[0]->bin(i).xWidth(); if (y_23 < ycut) { _h_R_Durham[0]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Durham[1]->numBins(); ++i) { double ycut = _h_R_Durham[1]->bin(i).xMid(); double width = _h_R_Durham[1]->bin(i).xWidth(); if (y_34 < ycut && y_23 > ycut) { _h_R_Durham[1]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Durham[2]->numBins(); ++i) { double ycut = _h_R_Durham[2]->bin(i).xMid(); double width = _h_R_Durham[2]->bin(i).xWidth(); if (y_45 < ycut && y_34 > ycut) { _h_R_Durham[2]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Durham[3]->numBins(); ++i) { double ycut = _h_R_Durham[3]->bin(i).xMid(); double width = _h_R_Durham[3]->bin(i).xWidth(); if (y_56 < ycut && y_45 > ycut) { _h_R_Durham[3]->fillBin(i, width); } } for (size_t i = 0; i < _h_R_Durham[4]->numBins(); ++i) { double ycut = _h_R_Durham[4]->bin(i).xMid(); double width = _h_R_Durham[4]->bin(i).xWidth(); if (y_56 > ycut) { _h_R_Durham[4]->fillBin(i, width); } } } } /// Finalize void finalize() { for (size_t n = 0; n < 4; ++n) normalize(_h_y_Durham[n]); for (size_t n = 0; n < 5; ++n) scale(_h_R_Jade[n], 100/sumOfWeights()); for (size_t n = 0; n < 5; ++n) scale(_h_R_Durham[n], 100/sumOfWeights()); } //@} private: /// @name Histograms //@{ Histo1DPtr _h_R_Jade[5]; Histo1DPtr _h_R_Durham[5]; Histo1DPtr _h_y_Durham[4]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(JADE_OPAL_2000_S4300807); } diff --git a/include/Rivet/Particle.hh b/include/Rivet/Particle.hh --- a/include/Rivet/Particle.hh +++ b/include/Rivet/Particle.hh @@ -1,700 +1,700 @@ // -*- 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) { } /// Constructor without GenParticle. Particle(PdgId pid, const FourMomentum& mom, const FourVector& pos=FourVector()) : ParticleBase(), _original(nullptr), _id(pid), _momentum(mom), _origin(pos) { } /// Constructor from a HepMC GenParticle pointer. Particle(const GenParticle* gp) : ParticleBase(), _original(gp), _id(gp->pdg_id()), _momentum(gp->momentum()) { const GenVertex* 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. Particle(const GenParticle& gp) : ParticleBase(), _original(&gp), _id(gp.pdg_id()), _momentum(gp.momentum()) { const GenVertex* vprod = gp.production_vertex(); if (vprod != nullptr) { setOrigin(vprod->position().t(), vprod->position().x(), vprod->position().y(), vprod->position().z()); } } //@} /// @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 const GenParticle* genParticle() const { return _original; } /// Cast operator for conversion to GenParticle* operator const GenParticle* () 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); } //@} /// @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()); } /// 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; //@} /// @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) 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. 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") + DEPRECATED("Too vague: use fromHadron() || fromPromptTau(), or isDirect()") 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 (may be null) const GenParticle* _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; }; /// @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/Tools/ParticleUtils.hh b/include/Rivet/Tools/ParticleUtils.hh --- a/include/Rivet/Tools/ParticleUtils.hh +++ b/include/Rivet/Tools/ParticleUtils.hh @@ -1,768 +1,768 @@ #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) /// 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) { return p.hasAncestorWith(f); } /// @brief Determine whether a particle has an ancestor which doesn't meet the function requirement inline bool hasAncestorWithout(const Particle& p, const ParticleSelector& f) { return p.hasAncestorWithout(f); } /// @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) { return p.hasDescendantWith(f); // return !p.allDescendants(f).empty(); } /// @brief Determine whether a particle has a descendant which doesn't meet the function requirement inline bool hasDescendantWithout(const Particle& p, const ParticleSelector& f) { return p.hasDescendantWithout(f); } /// @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(); } + // /// @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); } ParticleSelector fn; }; 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); } ParticleSelector fn; }; 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); } ParticleSelector fn; }; 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); } ParticleSelector fn; }; 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; } #endif diff --git a/src/AnalysisTools/MC_ParticleAnalysis.cc b/src/AnalysisTools/MC_ParticleAnalysis.cc --- a/src/AnalysisTools/MC_ParticleAnalysis.cc +++ b/src/AnalysisTools/MC_ParticleAnalysis.cc @@ -1,168 +1,168 @@ // -*- C++ -*- #include "Rivet/Analyses/MC_ParticleAnalysis.hh" namespace Rivet { - + MC_ParticleAnalysis::MC_ParticleAnalysis(const string& name, size_t nparticles, const string& particle_name) : Analysis(name), _nparts(nparticles), _pname(particle_name), _h_pt(nparticles), _h_eta(nparticles), _h_eta_plus(nparticles), _h_eta_minus(nparticles), _h_rap(nparticles), _h_rap_plus(nparticles), _h_rap_minus(nparticles), tmpeta(nparticles), tmprap(nparticles) { } // Book histograms void MC_ParticleAnalysis::init() { for (size_t i = 0; i < _nparts; ++i) { book(tmpeta[i], _pname + "_eta_pmratio_" + to_str(i+1)); book(tmprap[i], _pname + "_y_pmratio_" + to_str(i+1)); const string ptname = _pname + "_pt_" + to_str(i+1); const double ptmax = 1.0/(double(i)+2.0) * (sqrtS()>0.?sqrtS():14000.)/GeV/2.0; const int nbins_pt = 100/(i+1); book(_h_pt[i] ,ptname, logspace(nbins_pt, 1.0, ptmax)); const string etaname = _pname + "_eta_" + to_str(i+1); book(_h_eta[i] ,etaname, i > 1 ? 25 : 50, -5.0, 5.0); book(_h_eta_plus[i], "_" + etaname + "_plus", i > 1 ? 15 : 25, 0, 5); book(_h_eta_minus[i], "_" + etaname + "_minus", i > 1 ? 15 : 25, 0, 5); const string rapname = _pname + "_y_" + to_str(i+1); book(_h_rap[i] ,rapname, i > 1 ? 25 : 50, -5.0, 5.0); book(_h_rap_plus[i], "_" + rapname + "_plus", i > 1 ? 15 : 25, 0, 5); book(_h_rap_minus[i], "_" + rapname + "_minus", i > 1 ? 15 : 25, 0, 5); for (size_t j = i+1; j < min(size_t(3), _nparts); ++j) { const pair ij = std::make_pair(i, j); string detaname = _pname + "s_deta_" + to_str(i+1) + to_str(j+1); Histo1DPtr tmpeta; book(tmpeta, detaname, 25, -5.0, 5.0); _h_deta.insert(make_pair(ij, tmpeta)); string dphiname = _pname + "s_dphi_" + to_str(i+1) + to_str(j+1); Histo1DPtr tmpphi; book(tmpphi, dphiname, 25, 0.0, M_PI); _h_dphi.insert(make_pair(ij, tmpphi)); string dRname = _pname + "s_dR_" + to_str(i+1) + to_str(j+1); Histo1DPtr tmpR; book(tmpR, dRname, 25, 0.0, 5.0); _h_dR.insert(make_pair(ij, tmpR)); } } book(_h_multi_exclusive ,_pname + "_multi_exclusive", _nparts+3, -0.5, _nparts+3-0.5); book(_h_multi_inclusive ,_pname + "_multi_inclusive", _nparts+3, -0.5, _nparts+3-0.5); book(_h_multi_ratio, _pname + "_multi_ratio"); book(_h_multi_exclusive_prompt ,_pname + "_multi_exclusive_prompt", _nparts+3, -0.5, _nparts+3-0.5); book(_h_multi_inclusive_prompt ,_pname + "_multi_inclusive_prompt", _nparts+3, -0.5, _nparts+3-0.5); book(_h_multi_ratio_prompt, _pname + "_multi_ratio_prompt"); } // Do the analysis void MC_ParticleAnalysis::_analyze(const Event& event, const Particles& particles) { Particles promptparticles; for (const Particle& p : particles) - if (!p.fromDecay()) promptparticles += p; + if (p.isPrompt()) promptparticles += p; for (size_t i = 0; i < _nparts; ++i) { if (particles.size() < i+1) continue; _h_pt[i]->fill(particles[i].pt()/GeV); // Eta const double eta_i = particles[i].eta(); _h_eta[i]->fill(eta_i); (eta_i > 0.0 ? _h_eta_plus : _h_eta_minus)[i]->fill(fabs(eta_i)); // Rapidity const double rap_i = particles[i].rapidity(); _h_rap[i]->fill(rap_i); (rap_i > 0.0 ? _h_rap_plus : _h_rap_minus)[i]->fill(fabs(rap_i)); // Inter-particle properties for (size_t j = i+1; j < min(size_t(3),_nparts); ++j) { if (particles.size() < j+1) continue; std::pair ij = std::make_pair(i, j); double deta = particles[i].eta() - particles[j].eta(); double dphi = deltaPhi(particles[i].momentum(), particles[j].momentum()); double dR = deltaR(particles[i].momentum(), particles[j].momentum()); _h_deta[ij]->fill(deta); _h_dphi[ij]->fill(dphi); _h_dR[ij]->fill(dR); } } // Multiplicities _h_multi_exclusive->fill(particles.size()); _h_multi_exclusive_prompt->fill(promptparticles.size()); for (size_t i = 0; i < _nparts+2; ++i) { if (particles.size() >= i) _h_multi_inclusive->fill(i); if (promptparticles.size() >= i) _h_multi_inclusive_prompt->fill(i); } } // Finalize void MC_ParticleAnalysis::finalize() { const double scaling = crossSection()/sumOfWeights(); for (size_t i = 0; i < _nparts; ++i) { scale(_h_pt[i], scaling); scale(_h_eta[i], scaling); scale(_h_rap[i], scaling); // Create eta/rapidity ratio plots divide(_h_eta_plus[i], _h_eta_minus[i], tmpeta[i]); divide(_h_rap_plus[i], _h_rap_minus[i], tmprap[i]); } // Scale the d{eta,phi,R} histograms typedef map, Histo1DPtr> HistMap; for (HistMap::value_type& it : _h_deta) scale(it.second, scaling); for (HistMap::value_type& it : _h_dphi) scale(it.second, scaling); for (HistMap::value_type& it : _h_dR) scale(it.second, scaling); // Fill inclusive multi ratios for (size_t i = 0; i < _h_multi_inclusive->numBins()-1; ++i) { _h_multi_ratio->addPoint(i+1, 0, 0.5, 0); if (_h_multi_inclusive->bin(i).sumW() > 0.0) { const double ratio = _h_multi_inclusive->bin(i+1).sumW() / _h_multi_inclusive->bin(i).sumW(); const double relerr_i = _h_multi_inclusive->bin(i).relErr(); const double relerr_j = _h_multi_inclusive->bin(i+1).relErr(); const double err = ratio * (relerr_i + relerr_j); _h_multi_ratio->point(i).setY(ratio, err); } } for (size_t i = 0; i < _h_multi_inclusive_prompt->numBins()-1; ++i) { _h_multi_ratio_prompt->addPoint(i+1, 0, 0.5, 0); if (_h_multi_inclusive_prompt->bin(i).sumW() > 0.0) { const double ratio = _h_multi_inclusive_prompt->bin(i+1).sumW() / _h_multi_inclusive_prompt->bin(i).sumW(); const double relerr_i = _h_multi_inclusive_prompt->bin(i).relErr(); const double relerr_j = _h_multi_inclusive_prompt->bin(i+1).relErr(); const double err = ratio * (relerr_i + relerr_j); _h_multi_ratio_prompt->point(i).setY(ratio, err); } } scale(_h_multi_exclusive, scaling); scale(_h_multi_exclusive_prompt, scaling); scale(_h_multi_inclusive, scaling); scale(_h_multi_inclusive_prompt, scaling); } } diff --git a/src/Core/AnalysisHandler.cc b/src/Core/AnalysisHandler.cc --- a/src/Core/AnalysisHandler.cc +++ b/src/Core/AnalysisHandler.cc @@ -1,519 +1,520 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisHandler.hh" #include "Rivet/Analysis.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/BeamConstraint.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/Beam.hh" #include "YODA/ReaderYODA.h" #include "YODA/WriterYODA.h" #include #include using std::cout; using std::cerr; namespace { inline std::vector split(const std::string& input, const std::string& regex) { // passing -1 as the submatch index parameter performs splitting std::regex re(regex); std::sregex_token_iterator first{input.begin(), input.end(), re, -1}, last; return {first, last}; } } namespace Rivet { AnalysisHandler::AnalysisHandler(const string& runname) : _runname(runname), _eventCounter(vector(), Counter()), _xs(vector(), Scatter1D()), _initialised(false), _ignoreBeams(false), _defaultWeightIdx(0) {} AnalysisHandler::~AnalysisHandler() {} Log& AnalysisHandler::getLog() const { return Log::getLog("Rivet.Analysis.Handler"); } /// http://stackoverflow.com/questions/4654636/how-to-determine-if-a-string-is-a-number-with-c bool is_number(const std::string& s) { std::string::const_iterator it = s.begin(); while (it != s.end() && std::isdigit(*it)) ++it; return !s.empty() && it == s.end(); } /// Check if any of the weightnames is not a number bool AnalysisHandler::haveNamedWeights() { bool dec=false; for (unsigned int i=0;i<_weightNames.size();++i) { string s = _weightNames[i]; if (!is_number(s)) { dec=true; break; } } return dec; } void AnalysisHandler::init(const GenEvent& ge) { if (_initialised) throw UserError("AnalysisHandler::init has already been called: cannot re-initialize!"); // TODO TODO // should the rivet analysis objects know about weight names? setRunBeams(Rivet::beams(ge)); MSG_DEBUG("Initialising the analysis handler"); _eventNumber = ge.event_number(); setWeightNames(ge); if (haveNamedWeights()) MSG_INFO("Using named weights"); else MSG_INFO("NOT using named weights. Using first weight as nominal weight"); _eventCounter = CounterPtr(weightNames(), Counter("_EVTCOUNT")); // set the cross section based on what is reported by this event. // if no cross section if (ge.cross_section()) { MSG_TRACE("getting cross section."); double xs = ge.cross_section()->cross_section(); double xserr = ge.cross_section()->cross_section_error(); setCrossSection(xs, xserr); } // Check that analyses are beam-compatible, and remove those that aren't const size_t num_anas_requested = analysisNames().size(); vector anamestodelete; for (const AnaHandle a : _analyses) { if (!_ignoreBeams && !a->isCompatible(beams())) { //MSG_DEBUG(a->name() << " requires beams " << a->requiredBeams() << " @ " << a->requiredEnergies() << " GeV"); anamestodelete.push_back(a->name()); } } for (const string& aname : anamestodelete) { MSG_WARNING("Analysis '" << aname << "' is incompatible with the provided beams: removing"); removeAnalysis(aname); } if (num_anas_requested > 0 && analysisNames().empty()) { cerr << "All analyses were incompatible with the first event's beams\n" << "Exiting, since this probably wasn't intentional!" << '\n'; exit(1); } // Warn if any analysis' status is not unblemished for (const AnaHandle a : analyses()) { if (toUpper(a->status()) == "PRELIMINARY") { MSG_WARNING("Analysis '" << a->name() << "' is preliminary: be careful, it may change and/or be renamed!"); } else if (toUpper(a->status()) == "OBSOLETE") { MSG_WARNING("Analysis '" << a->name() << "' is obsolete: please update!"); } else if (toUpper(a->status()).find("UNVALIDATED") != string::npos) { MSG_WARNING("Analysis '" << a->name() << "' is unvalidated: be careful, it may be broken!"); } } // Initialize the remaining analyses _stage = Stage::INIT; for (AnaHandle a : _analyses) { MSG_DEBUG("Initialising analysis: " << a->name()); try { // Allow projection registration in the init phase onwards a->_allowProjReg = true; a->init(); //MSG_DEBUG("Checking consistency of analysis: " << a->name()); //a->checkConsistency(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::init method: " << err.what() << '\n'; exit(1); } MSG_DEBUG("Done initialising analysis: " << a->name()); } _stage = Stage::OTHER; _initialised = true; MSG_DEBUG("Analysis handler initialised"); } void AnalysisHandler::setWeightNames(const GenEvent& ge) { /// reroute the print output to a std::stringstream and process /// The iteration is done over a map in hepmc2 so this is safe std::ostringstream stream; ge.weights().print(stream); // Super lame, I know string str = stream.str(); std::regex re("(([^()]+))"); // Regex for stuff enclosed by parentheses () size_t idx = 0; for(std::sregex_iterator i = std::sregex_iterator(str.begin(), str.end(), re); i != std::sregex_iterator(); ++i ) { std::smatch m = *i; vector temp = ::split(m.str(), "[,]"); if (temp.size() ==2) { MSG_DEBUG("Name of weight #" << _weightNames.size() << ": " << temp[0]); // store the default weight based on weight names if (temp[0] == "Weight" || temp[0] == "0" || temp[0] == "Default") { MSG_DEBUG(_weightNames.size() << " is being used as the nominal."); _weightNames.push_back(""); _defaultWeightIdx = idx; } else _weightNames.push_back(temp[0]); idx++; } } } void AnalysisHandler::analyze(const GenEvent& ge) { // Call init with event as template if not already initialised if (!_initialised) init(ge); assert(_initialised); // Ensure that beam details match those from the first event (if we're checking beams) if ( !_ignoreBeams ) { const PdgIdPair beams = Rivet::beamIds(ge); const double sqrts = Rivet::sqrtS(ge); if (!compatible(beams, _beams) || !fuzzyEquals(sqrts, sqrtS())) { cerr << "Event beams mismatch: " << PID::toBeamsString(beams) << " @ " << sqrts/GeV << " GeV" << " vs. first beams " << this->beams() << " @ " << this->sqrtS()/GeV << " GeV" << '\n'; exit(1); } } // Create the Rivet event wrapper /// @todo Filter/normalize the event here Event event(ge); // set the cross section based on what is reported by this event. // if no cross section MSG_TRACE("getting cross section."); if (ge.cross_section()) { MSG_TRACE("getting cross section from GenEvent."); double xs = ge.cross_section()->cross_section(); double xserr = ge.cross_section()->cross_section_error(); setCrossSection(xs, xserr); } // won't happen for first event because _eventNumber is set in // init() if (_eventNumber != ge.event_number()) { /// @todo /// can we get away with not passing a matrix? MSG_TRACE("AnalysisHandler::analyze(): Pushing _eventCounter to persistent."); _eventCounter.get()->pushToPersistent(_subEventWeights); // if this is indeed a new event, push the temporary // histograms and reset for (const AnaHandle& a : _analyses) { for (auto ao : a->analysisObjects()) { MSG_TRACE("AnalysisHandler::analyze(): Pushing " << a->name() << "'s " << ao->name() << " to persistent."); ao.get()->pushToPersistent(_subEventWeights); } MSG_TRACE("AnalysisHandler::analyze(): finished pushing " << a->name() << "'s objects to persistent."); } _eventNumber = ge.event_number(); MSG_DEBUG("nominal event # " << _eventCounter.get()->_persistent[0]->numEntries()); MSG_DEBUG("nominal sum of weights: " << _eventCounter.get()->_persistent[0]->sumW()); MSG_DEBUG("Event has " << _subEventWeights.size() << " sub events."); _subEventWeights.clear(); } MSG_TRACE("starting new sub event"); _eventCounter.get()->newSubEvent(); for (const AnaHandle& a : _analyses) { for (auto ao : a->analysisObjects()) { ao.get()->newSubEvent(); } } _subEventWeights.push_back(event.weights()); MSG_DEBUG("Analyzing subevent #" << _subEventWeights.size() - 1 << "."); _eventCounter->fill(); // Run the analyses for (AnaHandle a : _analyses) { MSG_TRACE("About to run analysis " << a->name()); try { a->analyze(event); } catch (const Error& err) { cerr << "Error in " << a->name() << "::analyze method: " << err.what() << '\n'; exit(1); } MSG_TRACE("Finished running analysis " << a->name()); } } void AnalysisHandler::analyze(const GenEvent* ge) { if (ge == NULL) { MSG_ERROR("AnalysisHandler received null pointer to GenEvent"); //throw Error("AnalysisHandler received null pointer to GenEvent"); } analyze(*ge); } void AnalysisHandler::finalize() { if (!_initialised) return; MSG_INFO("Finalising analyses"); MSG_TRACE("AnalysisHandler::finalize(): Pushing analysis objects to persistent."); _eventCounter.get()->pushToPersistent(_subEventWeights); for (const AnaHandle& a : _analyses) { for (auto ao : a->analysisObjects()) ao.get()->pushToPersistent(_subEventWeights); } for (const AnaHandle& a : _analyses) { for (size_t iW = 0; iW < numWeights(); iW++) { _eventCounter.get()->setActiveWeightIdx(iW); _xs.get()->setActiveWeightIdx(iW); for (auto ao : a->analysisObjects()) ao.get()->setActiveWeightIdx(iW); MSG_TRACE("running " << a->name() << "::finalize() for weight " << iW << "."); try { a->finalize(); } catch (const Error& err) { cerr << "Error in " << a->name() << "::finalize method: " << err.what() << '\n'; exit(1); } } } // Print out number of events processed MSG_INFO("Processed " << numEvents() << " event" << (numEvents() == 1 ? "" : "s")); // Print out MCnet boilerplate cout << '\n'; cout << "The MCnet usage guidelines apply to Rivet: see http://www.montecarlonet.org/GUIDELINES" << '\n'; cout << "Please acknowledge plots made with Rivet analyses, and cite arXiv:1003.0694 (http://arxiv.org/abs/1003.0694)" << '\n'; } AnalysisHandler& AnalysisHandler::addAnalysis(const string& analysisname) { // Check for a duplicate analysis /// @todo Might we want to be able to run an analysis twice, with different params? /// Requires avoiding histo tree clashes, i.e. storing the histos on the analysis objects. for (const AnaHandle& a : _analyses) { if (a->name() == analysisname) { MSG_WARNING("Analysis '" << analysisname << "' already registered: skipping duplicate"); return *this; } } AnaHandle analysis( AnalysisLoader::getAnalysis(analysisname) ); if (analysis.get() != 0) { // < Check for null analysis. MSG_DEBUG("Adding analysis '" << analysisname << "'"); analysis->_analysishandler = this; _analyses.insert(analysis); } else { MSG_WARNING("Analysis '" << analysisname << "' not found."); } // MSG_WARNING(_analyses.size()); // for (const AnaHandle& a : _analyses) MSG_WARNING(a->name()); return *this; } AnalysisHandler& AnalysisHandler::removeAnalysis(const string& analysisname) { std::shared_ptr toremove; for (const AnaHandle a : _analyses) { if (a->name() == analysisname) { toremove = a; break; } } if (toremove.get() != 0) { MSG_DEBUG("Removing analysis '" << analysisname << "'"); _analyses.erase(toremove); } return *this; } void AnalysisHandler::addData(const std::vector& aos) { for (const YODA::AnalysisObjectPtr ao : aos) { const string path = ao->path(); if (path.size() > 1) { // path > "/" try { const string ananame = ::split(path, "/")[0]; AnaHandle a = analysis(ananame); //MultiweightAOPtr mao = ????; /// @todo generate right Multiweight object from ao //a->addAnalysisObject(mao); /// @todo Need to statistically merge... } catch (const Error& e) { MSG_WARNING(e.what()); } } } } void AnalysisHandler::readData(const string& filename) { vector aos; try { /// @todo Use new YODA SFINAE to fill the smart ptr vector directly vector aos_raw; YODA::ReaderYODA::read(filename, aos_raw); for (YODA::AnalysisObject* aor : aos_raw) aos.push_back(YODA::AnalysisObjectPtr(aor)); } catch (const YODA::ReadError & e) { throw UserError("Unexpected error in reading file: " + filename); } if (!aos.empty()) addData(aos); } vector AnalysisHandler::getRivetAOs() const { vector rtn; for (AnaHandle a : _analyses) { for (const auto & ao : a->analysisObjects()) { rtn.push_back(ao); } } rtn.push_back(_eventCounter); rtn.push_back(_xs); return rtn; } vector AnalysisHandler::getYodaAOs() const { vector rtn; for (auto rao : getRivetAOs()) { // need to set the index // before we can search the PATH rao.get()->setActiveWeightIdx(_defaultWeightIdx); if (rao->path().find("/TMP/") != string::npos) continue; for (size_t iW = 0; iW < numWeights(); iW++) { rao.get()->setActiveWeightIdx(iW); rtn.push_back(rao.get()->activeYODAPtr()); } } sort(rtn.begin(), rtn.end(), [](YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) { return a->path() < b->path(); } ); return rtn; } vector AnalysisHandler::getData() const { return getYodaAOs(); } void AnalysisHandler::writeData(const string& filename) const { const vector aos = getData(); try { YODA::WriterYODA::write(filename, aos.begin(), aos.end()); } catch ( YODA::WriteError ) { throw UserError("Unexpected error in writing file: " + filename); } } string AnalysisHandler::runName() const { return _runname; } size_t AnalysisHandler::numEvents() const { return _eventCounter->numEntries(); } std::vector AnalysisHandler::analysisNames() const { std::vector rtn; for (AnaHandle a : _analyses) { rtn.push_back(a->name()); } return rtn; } const AnaHandle AnalysisHandler::analysis(const std::string& analysisname) const { for (const AnaHandle a : analyses()) if (a->name() == analysisname) return a; throw Error("No analysis named '" + analysisname + "' registered in AnalysisHandler"); } AnalysisHandler& AnalysisHandler::addAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { //MSG_DEBUG("Adding analysis '" << aname << "'"); addAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::removeAnalyses(const std::vector& analysisnames) { for (const string& aname : analysisnames) { removeAnalysis(aname); } return *this; } AnalysisHandler& AnalysisHandler::setCrossSection(double xs, double xserr) { _xs = Scatter1DPtr(weightNames(), Scatter1D("_XSEC")); _eventCounter.get()->setActiveWeightIdx(_defaultWeightIdx); double nomwgt = sumOfWeights(); - // the cross section of each weight variation is the nominal cross section + // The cross section of each weight variation is the nominal cross section // times the sumOfWeights(variation) / sumOfWeights(nominal). - // this way the cross section will work correctly - for (int iW = 0; iW < numWeights(); iW++) { - _eventCounter.get()->setActiveWeightIdx(iW); - double s = sumOfWeights() / nomwgt; - _xs.get()->setActiveWeightIdx(iW); - _xs->addPoint(xs*s, xserr*s); + // This way the cross section will work correctly + for (size_t iW = 0; iW < numWeights(); iW++) { + _eventCounter.get()->setActiveWeightIdx(iW); + double s = sumOfWeights() / nomwgt; + _xs.get()->setActiveWeightIdx(iW); + _xs->addPoint(xs*s, xserr*s); } _eventCounter.get()->unsetActiveWeight(); _xs.get()->unsetActiveWeight(); return *this; } + AnalysisHandler& AnalysisHandler::addAnalysis(Analysis* analysis) { analysis->_analysishandler = this; _analyses.insert(AnaHandle(analysis)); return *this; } PdgIdPair AnalysisHandler::beamIds() const { return Rivet::beamIds(beams()); } double AnalysisHandler::sqrtS() const { return Rivet::sqrtS(beams()); } void AnalysisHandler::setIgnoreBeams(bool ignore) { _ignoreBeams=ignore; } } diff --git a/src/Projections/DressedLeptons.cc b/src/Projections/DressedLeptons.cc --- a/src/Projections/DressedLeptons.cc +++ b/src/Projections/DressedLeptons.cc @@ -1,137 +1,137 @@ // -*- C++ -*- #include "Rivet/Projections/DressedLeptons.hh" namespace Rivet { // On DressedLepton helper class //{ DressedLepton::DressedLepton(const Particle& dlepton) : Particle(dlepton) - { + { setConstituents({{dlepton}}); //< bare lepton is first constituent } DressedLepton::DressedLepton(const Particle& lepton, const Particles& photons, bool momsum) : Particle(lepton.pid(), lepton.momentum()) { setConstituents({{lepton}}); //< bare lepton is first constituent addConstituents(photons, momsum); } void DressedLepton::addPhoton(const Particle& p, bool momsum) { if (p.pid() != PID::PHOTON) throw Error("Clustering a non-photon on to a DressedLepton:"+to_string(p.pid())); addConstituent(p, momsum); } const Particle& DressedLepton::bareLepton() const { const Particle& l = constituents().front(); if (!l.isChargedLepton()) throw Error("First constituent of a DressedLepton is not a bare lepton: oops"); return l; } //} // Separate-FS version DressedLeptons::DressedLeptons(const FinalState& photons, const FinalState& bareleptons, double dRmax, const Cut& cut, bool useDecayPhotons) : FinalState(cut), _dRmax(dRmax), _fromDecay(useDecayPhotons) { setName("DressedLeptons"); IdentifiedFinalState photonfs(photons, PID::PHOTON); declare(photonfs, "Photons"); IdentifiedFinalState leptonfs(bareleptons); leptonfs.acceptIdPairs({PID::ELECTRON, PID::MUON, PID::TAU}); declare(leptonfs, "Leptons"); } // Single-FS version DressedLeptons::DressedLeptons(const FinalState& barefs, double dRmax, const Cut& cut, bool useDecayPhotons) : DressedLeptons(barefs, barefs, dRmax, cut, useDecayPhotons) { } CmpState DressedLeptons::compare(const Projection& p) const { // Compare the two as final states (for pT and eta cuts) const DressedLeptons& other = dynamic_cast(p); CmpState fscmp = FinalState::compare(other); if (fscmp != CmpState::EQ) return fscmp; const PCmp phcmp = mkNamedPCmp(p, "Photons"); if (phcmp != CmpState::EQ) return phcmp; const PCmp sigcmp = mkNamedPCmp(p, "Leptons"); if (sigcmp != CmpState::EQ) return sigcmp; return (cmp(_dRmax, other._dRmax) || cmp(_fromDecay, other._fromDecay)); } void DressedLeptons::project(const Event& e) { _theParticles.clear(); // Get bare leptons const FinalState& signal = applyProjection(e, "Leptons"); Particles bareleptons = signal.particles(); if (bareleptons.empty()) return; // Initialise DL collection with bare leptons vector allClusteredLeptons; allClusteredLeptons.reserve(bareleptons.size()); for (const Particle& bl : bareleptons) { Particle dl(bl.pid(), bl.momentum()); dl.setConstituents({bl}); allClusteredLeptons += dl; } // If the radius is 0 or negative, don't even attempt to cluster if (_dRmax > 0) { // Match each photon to its closest charged lepton within the dR cone const FinalState& photons = applyProjection(e, "Photons"); for (const Particle& photon : photons.particles()) { // Ignore photon if it's from a hadron/tau decay and we're avoiding those - if (!_fromDecay && photon.fromDecay()) continue; + if (!_fromDecay && !photon.isPrompt()) continue; const FourMomentum& p_P = photon.momentum(); double dRmin = _dRmax; int idx = -1; for (size_t i = 0; i < bareleptons.size(); ++i) { const Particle& bl = bareleptons[i]; // Only cluster photons around *charged* signal particles if (bl.charge3() == 0) continue; // Find the closest lepton double dR = deltaR(bl, p_P); if (dR < dRmin) { dRmin = dR; idx = i; } } if (idx > -1) allClusteredLeptons[idx].addConstituent(photon, true); } } // Fill the canonical particles collection with the composite DL Particles for (const Particle& lepton : allClusteredLeptons) { const bool acc = accept(lepton); - MSG_TRACE("Clustered lepton " << lepton - << " with constituents = " << lepton.constituents() + MSG_TRACE("Clustered lepton " << lepton + << " with constituents = " << lepton.constituents() << ", cut-pass = " << std::boolalpha << acc); if (acc) _theParticles.push_back(lepton); } - MSG_DEBUG("#dressed leptons = " << allClusteredLeptons.size() + MSG_DEBUG("#dressed leptons = " << allClusteredLeptons.size() << " -> " << _theParticles.size() << " after cuts"); } } diff --git a/src/Projections/FastJets.cc b/src/Projections/FastJets.cc --- a/src/Projections/FastJets.cc +++ b/src/Projections/FastJets.cc @@ -1,215 +1,215 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/HeavyHadrons.hh" #include "Rivet/Projections/TauFinder.hh" namespace Rivet { void FastJets::_initBase() { setName("FastJets"); declare(HeavyHadrons(), "HFHadrons"); declare(TauFinder(TauFinder::DecayMode::HADRONIC), "Taus"); } void FastJets::_initJdef(Algo alg, double rparameter, double seed_threshold) { MSG_DEBUG("JetAlg = " << static_cast(alg)); MSG_DEBUG("R parameter = " << rparameter); MSG_DEBUG("Seed threshold = " << seed_threshold); if (alg == KT) { _jdef = fastjet::JetDefinition(fastjet::kt_algorithm, rparameter, fastjet::E_scheme); } else if (alg == CAM) { _jdef = fastjet::JetDefinition(fastjet::cambridge_algorithm, rparameter, fastjet::E_scheme); } else if (alg == ANTIKT) { _jdef = fastjet::JetDefinition(fastjet::antikt_algorithm, rparameter, fastjet::E_scheme); } else if (alg == DURHAM) { _jdef = fastjet::JetDefinition(fastjet::ee_kt_algorithm, fastjet::E_scheme); } else if (alg == GENKTEE) { _jdef = fastjet::JetDefinition(fastjet::ee_genkt_algorithm, rparameter, -1); } else { // Plugins: if (alg == SISCONE) { const double OVERLAP_THRESHOLD = 0.75; _plugin.reset(new fastjet::SISConePlugin(rparameter, OVERLAP_THRESHOLD)); // } else if (alg == PXCONE) { // string msg = "PxCone currently not supported, since FastJet doesn't install it by default. "; // msg += "Please notify the Rivet authors if this behaviour should be changed."; // throw Error(msg); // _plugin.reset(new fastjet::PxConePlugin(rparameter)); } else if (alg == ATLASCONE) { const double OVERLAP_THRESHOLD = 0.5; _plugin.reset(new fastjet::ATLASConePlugin(rparameter, seed_threshold, OVERLAP_THRESHOLD)); } else if (alg == CMSCONE) { _plugin.reset(new fastjet::CMSIterativeConePlugin(rparameter, seed_threshold)); } else if (alg == CDFJETCLU) { const double OVERLAP_THRESHOLD = 0.75; _plugin.reset(new fastjet::CDFJetCluPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == CDFMIDPOINT) { const double OVERLAP_THRESHOLD = 0.5; _plugin.reset(new fastjet::CDFMidPointPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == D0ILCONE) { const double min_jet_Et = 6.0; _plugin.reset(new fastjet::D0RunIIConePlugin(rparameter, min_jet_Et)); } else if (alg == JADE) { _plugin.reset(new fastjet::JadePlugin()); } else if (alg == TRACKJET) { _plugin.reset(new fastjet::TrackJetPlugin(rparameter)); } _jdef = fastjet::JetDefinition(_plugin.get()); } } CmpState FastJets::compare(const Projection& p) const { const FastJets& other = dynamic_cast(p); return \ cmp(_useMuons, other._useMuons) || cmp(_useInvisibles, other._useInvisibles) || mkNamedPCmp(other, "FS") || cmp(_jdef.jet_algorithm(), other._jdef.jet_algorithm()) || cmp(_jdef.recombination_scheme(), other._jdef.recombination_scheme()) || cmp(_jdef.plugin(), other._jdef.plugin()) || cmp(_jdef.R(), other._jdef.R()) || cmp(_adef, other._adef); } // STATIC PseudoJets FastJets::mkClusterInputs(const Particles& fsparticles, const Particles& tagparticles) { PseudoJets pjs; /// @todo Use FastJet3's UserInfo system to store Particle pointers directly? // Store 4 vector data about each particle into FastJet's PseudoJets for (size_t i = 0; i < fsparticles.size(); ++i) { fastjet::PseudoJet pj = fsparticles[i]; pj.set_user_index(i+1); pjs.push_back(pj); } // And the same for ghost tagging particles (with negative user indices) for (size_t i = 0; i < tagparticles.size(); ++i) { fastjet::PseudoJet pj = tagparticles[i]; pj *= 1e-20; ///< Ghostify the momentum pj.set_user_index(-i-1); pjs.push_back(pj); } return pjs; } // STATIC Jet FastJets::mkJet(const PseudoJet& pj, const Particles& fsparticles, const Particles& tagparticles) { const PseudoJets pjconstituents = pj.constituents(); Particles constituents, tags; constituents.reserve(pjconstituents.size()); for (const fastjet::PseudoJet& pjc : pjconstituents) { // Pure ghosts don't have corresponding particles if (pjc.has_area() && pjc.is_pure_ghost()) continue; // Default user index = 0 doesn't give valid particle lookup if (pjc.user_index() == 0) continue; // Split by index sign into constituent & tag lookup if (pjc.user_index() > 0) { // Find constituents if index > 0 const size_t i = pjc.user_index() - 1; if (i >= fsparticles.size()) throw RangeError("FS particle lookup failed in jet construction"); constituents.push_back(fsparticles.at(i)); } else if (!tagparticles.empty()) { // Find tags if index < 0 const size_t i = abs(pjc.user_index()) - 1; if (i >= tagparticles.size()) throw RangeError("Tag particle lookup failed in jet construction"); tags.push_back(tagparticles.at(i)); } } return Jet(pj, constituents, tags); } // STATIC Jets FastJets::mkJets(const PseudoJets& pjs, const Particles& fsparticles, const Particles& tagparticles) { Jets rtn; rtn.reserve(pjs.size()); for (const PseudoJet pj : pjs) { rtn.push_back(FastJets::mkJet(pj, fsparticles, tagparticles)); } return rtn; } void FastJets::project(const Event& e) { // Assemble final state particles const string fskey = (_useInvisibles == JetAlg::Invisibles::NONE) ? "VFS" : "FS"; Particles fsparticles = applyProjection(e, fskey).particles(); // Remove prompt invisibles if needed (already done by VFS if using NO_INVISIBLES) if (_useInvisibles == JetAlg::Invisibles::DECAY) { - ifilter_discard(fsparticles, [](const Particle& p) { return !(p.isVisible() || p.fromDecay()); }); + ifilter_discard(fsparticles, [](const Particle& p) { return !p.isVisible() && p.isPrompt(); }); } // Remove prompt/all muons if needed if (_useMuons == JetAlg::Muons::DECAY) { - ifilter_discard(fsparticles, [](const Particle& p) { return isMuon(p) && !p.fromDecay(); }); + ifilter_discard(fsparticles, [](const Particle& p) { return isMuon(p) && p.isPrompt(); }); } else if (_useMuons == JetAlg::Muons::NONE) { ifilter_discard(fsparticles, isMuon); } // Tagging particles const Particles chadrons = applyProjection(e, "HFHadrons").cHadrons(); const Particles bhadrons = applyProjection(e, "HFHadrons").bHadrons(); const Particles taus = applyProjection(e, "Taus").particles(); calc(fsparticles, chadrons+bhadrons+taus); } void FastJets::calc(const Particles& fsparticles, const Particles& tagparticles) { MSG_DEBUG("Finding jets from " << fsparticles.size() << " input particles + " << tagparticles.size() << " tagging particles"); _fsparticles = fsparticles; _tagparticles = tagparticles; // Make pseudojets, with mapping info to Rivet FS and tag particles PseudoJets pjs = mkClusterInputs(_fsparticles, _tagparticles); // Run either basic or area-calculating cluster sequence as reqd. if (_adef) { _cseq.reset(new fastjet::ClusterSequenceArea(pjs, _jdef, *_adef)); } else { _cseq.reset(new fastjet::ClusterSequence(pjs, _jdef)); } MSG_DEBUG("ClusterSequence constructed; Njets_tot = " << _cseq->inclusive_jets().size() << ", Njets(pT > 10 GeV) = " << _cseq->inclusive_jets(10*GeV).size()); } void FastJets::reset() { _yscales.clear(); _fsparticles.clear(); _tagparticles.clear(); /// @todo _cseq = fastjet::ClusterSequence(); } Jets FastJets::_jets() const { /// @todo Cache? return mkJets(pseudojets(), _fsparticles, _tagparticles); } /// @todo "Automate" trimming as part of project() with pre-registered Filters Jet FastJets::trimJet(const Jet& input, const fastjet::Filter& trimmer) const { if (input.pseudojet().associated_cluster_sequence() != clusterSeq().get()) throw Error("To trim a Rivet::Jet, its associated PseudoJet must have come from this FastJets' ClusterSequence"); PseudoJet pj = trimmer(input); return mkJet(pj, _fsparticles, _tagparticles); } PseudoJets FastJets::pseudoJets(double ptmin) const { return clusterSeq() ? clusterSeq()->inclusive_jets(ptmin) : PseudoJets(); } } diff --git a/src/Projections/FinalPartons.cc b/src/Projections/FinalPartons.cc --- a/src/Projections/FinalPartons.cc +++ b/src/Projections/FinalPartons.cc @@ -1,42 +1,42 @@ // -*- C++ -*- #include "Rivet/Projections/FinalPartons.hh" #include "Rivet/Tools/RivetHepMC.hh" namespace Rivet { bool FinalPartons::accept(const Particle& p) const { // Reject if *not* a parton if (!isParton(p)) return false; // Accept partons if they end on a standard hadronization vertex if (p.genParticle()->end_vertex() != nullptr && p.genParticle()->end_vertex()->id() == 5) return true; // Reject if p has a parton child. for (const Particle& c : p.children()) if (isParton(c)) return false; // Reject if from a hadron or tau decay. - if (p.fromDecay()) + if (p.fromHadron() || p.fromTau()) return false; return _cuts->accept(p); } void FinalPartons::project(const Event& e) { _theParticles.clear(); for (const GenParticle* gp : Rivet::particles(e.genEvent())) { if (!gp) continue; const Particle p(gp); if (accept(p)) _theParticles.push_back(p); } } }