diff --git a/analyses/pluginALICE/ALICE_2014_I1300380.cc b/analyses/pluginALICE/ALICE_2014_I1300380.cc --- a/analyses/pluginALICE/ALICE_2014_I1300380.cc +++ b/analyses/pluginALICE/ALICE_2014_I1300380.cc @@ -1,119 +1,119 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/UnstableFinalState.hh" namespace Rivet { class ALICE_2014_I1300380 : public Analysis { public: ALICE_2014_I1300380() : Analysis("ALICE_2014_I1300380") {} public: void init() { const UnstableFinalState cfs(Cuts::absrap<0.5); declare(cfs, "CFS"); // Plots from the paper book(_histPtSigmaStarPlus ,"d01-x01-y01"); // Sigma*+ book(_histPtSigmaStarMinus ,"d01-x01-y02"); // Sigma*- book(_histPtSigmaStarPlusAnti ,"d01-x01-y03"); // anti Sigma*- book(_histPtSigmaStarMinusAnti ,"d01-x01-y04"); // anti Sigma*+ book(_histPtXiStar ,"d02-x01-y01"); // 0.5 * (xi star + anti xi star) book(_histAveragePt ,"d03-x01-y01"); // profile } void analyze(const Event& event) { const UnstableFinalState& cfs = apply(event, "CFS"); for (const Particle& p : cfs.particles()) { // protections against mc generators decaying long-lived particles if ( !(p.hasAncestor(310) || p.hasAncestor(-310) || // K0s p.hasAncestor(130) || p.hasAncestor(-130) || // K0l p.hasAncestor(3322) || p.hasAncestor(-3322) || // Xi0 p.hasAncestor(3122) || p.hasAncestor(-3122) || // Lambda p.hasAncestor(3222) || p.hasAncestor(-3222) || // Sigma+/- p.hasAncestor(3312) || p.hasAncestor(-3312) || // Xi-/+ p.hasAncestor(3334) || p.hasAncestor(-3334) )) // Omega-/+ { - int aid = abs(p.pdgId()); + int aid = abs(p.pid()); if (aid == 211 || // pi+ aid == 321 || // K+ aid == 313 || // K*(892)0 aid == 2212 || // proton aid == 333 ) { // phi(1020) _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); } } // end if "rejection of long-lived particles" - switch (p.pdgId()) { + switch (p.pid()) { case 3224: _histPtSigmaStarPlus->fill(p.pT()/GeV); _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case -3224: _histPtSigmaStarPlusAnti->fill(p.pT()/GeV); _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case 3114: _histPtSigmaStarMinus->fill(p.pT()/GeV); _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case -3114: _histPtSigmaStarMinusAnti->fill(p.pT()/GeV); _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case 3324: _histPtXiStar->fill(p.pT()/GeV); _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case -3324: _histPtXiStar->fill(p.pT()/GeV); _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case 3312: _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case -3312: _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case 3334: _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; case -3334: _histAveragePt->fill(p.mass()/GeV, p.pT()/GeV); break; } } } void finalize() { scale(_histPtSigmaStarPlus, 1./sumOfWeights()); scale(_histPtSigmaStarPlusAnti, 1./sumOfWeights()); scale(_histPtSigmaStarMinus, 1./sumOfWeights()); scale(_histPtSigmaStarMinusAnti, 1./sumOfWeights()); scale(_histPtXiStar, 1./sumOfWeights()/ 2.); } private: // plots from the paper Histo1DPtr _histPtSigmaStarPlus; Histo1DPtr _histPtSigmaStarPlusAnti; Histo1DPtr _histPtSigmaStarMinus; Histo1DPtr _histPtSigmaStarMinusAnti; Histo1DPtr _histPtXiStar; Profile1DPtr _histAveragePt; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2014_I1300380); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1298811.cc b/analyses/pluginATLAS/ATLAS_2014_I1298811.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1298811.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1298811.cc @@ -1,196 +1,196 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2014_I1298811 : public Analysis { public: ATLAS_2014_I1298811() : Analysis("ATLAS_2014_I1298811") { } void init() { // Configure projections const FinalState fs(-4.8, 4.8, 0*MeV); declare(fs, "FS"); const FastJets jets(fs, FastJets::ANTIKT, 0.4); declare(jets, "Jets"); // Book histograms for (size_t itopo = 0; itopo < 2; ++itopo) { // Profiles for (size_t iregion = 0; iregion < 3; ++iregion) { book(_p_ptsumch_vs_ptlead[itopo][iregion] ,1+iregion, 1, itopo+1); book(_p_nch_vs_ptlead[itopo][iregion] ,4+iregion, 1, itopo+1); } book(_p_etsum25_vs_ptlead_trans[itopo] ,7, 1, itopo+1); book(_p_etsum48_vs_ptlead_trans[itopo] ,8, 1, itopo+1); book(_p_chratio_vs_ptlead_trans[itopo] ,9, 1, itopo+1); book(_p_ptmeanch_vs_ptlead_trans[itopo] ,10, 1, itopo+1); // 1D histos for (size_t iregion = 0; iregion < 3; ++iregion) { for (size_t ipt = 0; ipt < 4; ++ipt) { book(_h_ptsumch[ipt][itopo][iregion] ,13+3*ipt+iregion, 1, itopo+1); book(_h_nch[ipt][itopo][iregion] ,25+3*ipt+iregion, 1, itopo+1); } } } book(_p_ptmeanch_vs_nch_trans[0], 11, 1, 1); book(_p_ptmeanch_vs_nch_trans[1], 12, 1, 1); } void analyze(const Event& event) { // Find the jets with pT > 20 GeV and *rapidity* within 2.8 /// @todo Use Cuts instead rather than an eta cut in the proj and a y cut after const Jets alljets = apply(event, "Jets").jetsByPt(20*GeV); Jets jets; for (const Jet& j : alljets) if (j.absrap() < 2.8) jets.push_back(j); // Require at least one jet in the event if (jets.empty()) vetoEvent; // Identify the leading jet and its phi and pT const FourMomentum plead = jets[0].momentum(); const double philead = plead.phi(); const double etalead = plead.eta(); const double ptlead = plead.pT(); MSG_DEBUG("Leading object: pT = " << ptlead << ", eta = " << etalead << ", phi = " << philead); // Sum particle properties in the transverse regions int tmpnch[2] = {0,0}; double tmpptsum[2] = {0,0}; double tmpetsum48[2] = {0,0}; double tmpetsum25[2] = {0,0}; const Particles particles = apply(event, "FS").particles(); for (const Particle& p : particles) { // Only consider the transverse region(s), not toward or away if (!inRange(deltaPhi(p.phi(), philead), PI/3.0, TWOPI/3.0)) continue; // Work out which transverse side this particle is on const size_t iside = (mapAngleMPiToPi(p.phi() - philead) > 0) ? 0 : 1; MSG_TRACE(p.phi() << " vs. " << philead << ": " << iside); // Charged or neutral particle? - const bool charged = PID::charge3(p.pdgId()) != 0; + const bool charged = PID::charge3(p.pid()) != 0; // Track observables if (charged && fabs(p.eta()) < 2.5 && p.pT() > 500*MeV) { tmpnch[iside] += 1; tmpptsum[iside] += p.pT(); } // Cluster observables if ((charged && p.p3().mod() > 200*MeV) || (!charged && p.p3().mod() > 500*MeV)) { tmpetsum48[iside] += p.pT(); if (fabs(p.eta()) < 2.5) tmpetsum25[iside] += p.pT(); } } // Construct tot/max/min counts (for trans/max/min, indexed by iregion) const int nch[3] = { tmpnch[0] + tmpnch[1], std::max(tmpnch[0], tmpnch[1]), std::min(tmpnch[0], tmpnch[1]) }; const double ptsum[3] = { tmpptsum[0] + tmpptsum[1], std::max(tmpptsum[0], tmpptsum[1]), std::min(tmpptsum[0], tmpptsum[1]) }; const double etsum48[3] = { tmpetsum48[0] + tmpetsum48[1], std::max(tmpetsum48[0], tmpetsum48[1]), std::min(tmpetsum48[0], tmpetsum48[1]) }; const double etsum25[3] = { tmpetsum25[0] + tmpetsum25[1], std::max(tmpetsum25[0], tmpetsum25[1]), std::min(tmpetsum25[0], tmpetsum25[1]) }; ////////////////////////////////////////////////////////// // Now fill the histograms with the computed quantities // phi sizes of each trans/max/min region (for indexing by iregion) const double dphi[3] = { 2*PI/3.0, PI/3.0, PI/3.0 }; // Loop over inclusive jet and exclusive dijet configurations for (size_t itopo = 0; itopo < 2; ++itopo) { // Exit early if in the exclusive dijet iteration and the exclusive dijet cuts are not met if (itopo == 1) { if (jets.size() != 2) continue; const FourMomentum psublead = jets[1].momentum(); // Delta(phi) cut const double phisublead = psublead.phi(); if (deltaPhi(philead, phisublead) < 2.5) continue; // pT fraction cut const double ptsublead = psublead.pT(); if (ptsublead < 0.5*ptlead) continue; MSG_DEBUG("Exclusive dijet event"); } // Plot profiles and distributions which have no max/min region definition _p_etsum25_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum25[0]/5.0/dphi[0]/GeV); _p_etsum48_vs_ptlead_trans[itopo]->fill(ptlead/GeV, etsum48[0]/9.6/dphi[0]/GeV); if (etsum25[0] > 0) { _p_chratio_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptsum[0]/etsum25[0]); } const double ptmean = safediv(ptsum[0], nch[0], -1); ///< Return -1 if div by zero if (ptmean >= 0) { _p_ptmeanch_vs_ptlead_trans[itopo]->fill(ptlead/GeV, ptmean/GeV); _p_ptmeanch_vs_nch_trans[itopo]->fill(nch[0], ptmean/GeV); } // Plot remaining profile and 1D observables, which are defined in all 3 tot/max/min regions for (size_t iregion = 0; iregion < 3; ++iregion) { _p_ptsumch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, ptsum[iregion]/5.0/dphi[iregion]/GeV); _p_nch_vs_ptlead[itopo][iregion]->fill(ptlead/GeV, nch[iregion]/5.0/dphi[iregion]); for (size_t ipt = 0; ipt < 4; ++ipt) { if (ipt == 1 && !inRange(ptlead/GeV, 20, 60)) continue; if (ipt == 2 && !inRange(ptlead/GeV, 60, 210)) continue; if (ipt == 3 && ptlead/GeV < 210) continue; _h_ptsumch[ipt][itopo][iregion]->fill(ptsum[iregion]/5.0/dphi[iregion]/GeV); _h_nch[ipt][itopo][iregion]->fill(nch[iregion]/5.0/dphi[iregion]); } } } } void finalize() { for (size_t iregion = 0; iregion < 3; ++iregion) { for (size_t itopo = 0; itopo < 2; ++itopo) { for (size_t ipt = 0; ipt < 4; ++ipt) { normalize(_h_ptsumch[ipt][itopo][iregion], 1.0); normalize(_h_nch[ipt][itopo][iregion], 1.0); } } } } private: /// @name Histogram arrays //@{ Profile1DPtr _p_ptsumch_vs_ptlead[2][3]; Profile1DPtr _p_nch_vs_ptlead[2][3]; Profile1DPtr _p_ptmeanch_vs_ptlead_trans[2]; Profile1DPtr _p_etsum25_vs_ptlead_trans[2]; Profile1DPtr _p_etsum48_vs_ptlead_trans[2]; Profile1DPtr _p_chratio_vs_ptlead_trans[2]; Profile1DPtr _p_ptmeanch_vs_nch_trans[2]; Histo1DPtr _h_ptsumch[4][2][3]; Histo1DPtr _h_nch[4][2][3]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1298811); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1625109.cc b/analyses/pluginATLAS/ATLAS_2017_I1625109.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1625109.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1625109.cc @@ -1,307 +1,307 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2017_I1625109 : public Analysis { public: /// Constructor /// @brief measurement of on-shell ZZ at 13 TeV DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1625109); /// @name Analysis methods //@{ struct Dilepton { Dilepton() {}; Dilepton(const ParticlePair & _leptons) : leptons(_leptons) {} FourMomentum momentum() const { return leptons.first.mom() + leptons.second.mom(); } ParticlePair leptons; }; struct Quadruplet { vector getLeptonsSortedByPt() const { vector out = { leadingDilepton.leptons.first, leadingDilepton.leptons.second, subleadingDilepton.leptons.first, subleadingDilepton.leptons.second }; std::sort(out.begin(), out.end(), cmpMomByPt); return out; } Quadruplet(const Dilepton& dilepton1, const Dilepton& dilepton2) { if (dilepton1.momentum().pt() > dilepton2.momentum().pt()) { leadingDilepton = dilepton1; subleadingDilepton = dilepton2; } else { leadingDilepton = dilepton2; subleadingDilepton = dilepton1; } leptonsSortedByPt = getLeptonsSortedByPt(); } FourMomentum momentum() const { return leadingDilepton.momentum() + subleadingDilepton.momentum(); } double distanceFromZMass() const { return abs(leadingDilepton.momentum().mass() - Z_mass) + abs(subleadingDilepton.momentum().mass() - Z_mass); } Dilepton leadingDilepton; Dilepton subleadingDilepton; vector leptonsSortedByPt; }; typedef vector Quadruplets; typedef std::pair IndexPair; vector getOppositeChargePairsIndices(const vector& leptons) { vector indices = {}; if (leptons.size() < 2) return indices; for (size_t i = 0; i < leptons.size(); ++i) { for (size_t k = i+1; k < leptons.size(); ++k) { const auto charge_i = leptons.at(i).charge(); const auto charge_k = leptons.at(k).charge(); if (charge_i == -charge_k) { indices.push_back(std::make_pair(i, k)); } } } return indices; } bool indicesOverlap(IndexPair a, IndexPair b) { return (a.first == b.first || a.first == b.second || a.second == b.first || a.second == b.second); } bool passesHierarchicalPtRequirements(const Quadruplet& quadruplet) { const auto& sorted_leptons = quadruplet.leptonsSortedByPt; if (sorted_leptons.at(0).pt() < 20*GeV) return false; if (sorted_leptons.at(1).pt() < 15*GeV) return false; if (sorted_leptons.at(2).pt() < 10*GeV) return false; return true; } bool passesDileptonMinimumMassRequirement(const Quadruplet& quadruplet) { const auto& leptons = quadruplet.leptonsSortedByPt; for (const auto& l1 : leptons) { for (const auto& l2 : leptons) { if (l1 == l2) continue; - if ((l1.pdgId() + l2.pdgId() == 0) && ((l1.mom() + l2.mom()).mass() < 5.0*GeV)) return false; + if ((l1.pid() + l2.pid() == 0) && ((l1.mom() + l2.mom()).mass() < 5.0*GeV)) return false; } } return true; } bool passesLeptonDeltaRRequirements(const Quadruplet& quadruplet) { const auto& leptons = quadruplet.leptonsSortedByPt; for (const auto& l1 : leptons) { for (const auto& l2 : leptons) { if (l1 == l2) continue; // Any lepton flavour: if (deltaR(l1.mom(), l2.mom()) < 0.1) return false; // Different lepton flavour: if ((l1.abspid() - l2.abspid() != 0) && (deltaR(l1.mom(), l2.mom()) < 0.2)) return false; } } return true; } Quadruplets formQuadrupletsByChannel(const vector& same_flavour_leptons, vector indices) { Quadruplets quadruplets = {}; for (size_t i = 0; i < indices.size(); ++i) { for (size_t k = i+1; k < indices.size(); ++k) { const auto& pair_i = indices.at(i); const auto& pair_k = indices.at(k); if (indicesOverlap(pair_i, pair_k)) continue; const auto d1 = Dilepton({same_flavour_leptons.at(pair_i.first), same_flavour_leptons.at(pair_i.second)}); const auto d2 = Dilepton({same_flavour_leptons.at(pair_k.first), same_flavour_leptons.at(pair_k.second)}); const auto quadruplet = Quadruplet(d1, d2); if (passesHierarchicalPtRequirements(quadruplet)) quadruplets.push_back(quadruplet); } } return quadruplets; } Quadruplets formQuadrupletsByChannel(const vector& electrons, vector e_indices, const vector& muons, vector m_indices) { Quadruplets quadruplets = {}; for (const auto& pair_e : e_indices) { for (const auto& pair_m : m_indices) { const auto d1 = Dilepton({electrons.at(pair_e.first), electrons.at(pair_e.second)}); const auto d2 = Dilepton({muons.at(pair_m.first), muons.at(pair_m.second)}); const auto quadruplet = Quadruplet(d1, d2); if (passesHierarchicalPtRequirements(quadruplet)) quadruplets.push_back(quadruplet); } } return quadruplets; } Quadruplets getQuadruplets(const vector& electrons, const vector& muons) { const auto oc_electrons_indices = getOppositeChargePairsIndices(electrons); const auto oc_muons_indices = getOppositeChargePairsIndices(muons); const auto electron_quadruplets = formQuadrupletsByChannel(electrons, oc_electrons_indices); const auto muon_quadruplets = formQuadrupletsByChannel(muons, oc_muons_indices); const auto mixed_quadruplets = formQuadrupletsByChannel(electrons, oc_electrons_indices, muons, oc_muons_indices); auto quadruplets = electron_quadruplets; quadruplets.insert(quadruplets.end(), muon_quadruplets.begin(), muon_quadruplets.end()); quadruplets.insert(quadruplets.end(), mixed_quadruplets.begin(), mixed_quadruplets.end()); return quadruplets; } Quadruplet selectQuadruplet(const Quadruplets& quadruplets) { if (quadruplets.empty()) throw std::logic_error("Expect at least one quadruplet! The user should veto events without quadruplets."); Quadruplets sortedQuadruplets = quadruplets; std::sort(sortedQuadruplets.begin(), sortedQuadruplets.end(), [](const Quadruplet& a, const Quadruplet& b) { return a.distanceFromZMass() < b.distanceFromZMass(); }); return sortedQuadruplets.at(0); } /// Book histograms and initialise projections before the run void init() { const Cut presel = Cuts::abseta < 5 && Cuts::pT > 100*MeV; const FinalState fs(presel); // Prompt leptons, photons, neutrinos // Excluding those from tau decay const PromptFinalState photons(presel && Cuts::abspid == PID::PHOTON, false); const PromptFinalState bare_elecs(presel && Cuts::abspid == PID::ELECTRON, false); const PromptFinalState bare_muons(presel && Cuts::abspid == PID::MUON, false); // Baseline lepton and jet declaration const Cut lepton_baseline_cuts = Cuts::abseta < 2.7 && Cuts::pT > 5*GeV; const DressedLeptons elecs = DressedLeptons(photons, bare_elecs, 0.1, lepton_baseline_cuts); const DressedLeptons muons = DressedLeptons(photons, bare_muons, 0.1, lepton_baseline_cuts); declare(elecs, "electrons"); declare(muons, "muons"); VetoedFinalState jet_input(fs); jet_input.addVetoOnThisFinalState(elecs); jet_input.addVetoOnThisFinalState(muons); declare(FastJets(jet_input, FastJets::ANTIKT, 0.4), "jets"); // // Book histograms book(_h["pT_4l"], 2, 1, 1); book(_h["pT_leading_dilepton"], 8, 1, 1); book(_h["pT_subleading_dilepton"], 14, 1, 1); book(_h["pT_lepton1"], 20, 1, 1); book(_h["pT_lepton2"], 26, 1, 1); book(_h["pT_lepton3"], 32, 1, 1); book(_h["pT_lepton4"], 38, 1, 1); book(_h["absy_4l"], 44, 1, 1); book(_h["deltay_dileptons"], 50, 1, 1); book(_h["deltaphi_dileptons"], 56, 1, 1); book(_h["N_jets"], 62, 1, 1); book(_h["N_central_jets"], 68, 1, 1); book(_h["N_jets60"], 74, 1, 1); book(_h["mass_dijet"], 80, 1, 1); book(_h["deltay_dijet"], 86, 1, 1); book(_h["scalarpTsum_jets"], 92, 1, 1); book(_h["abseta_jet1"], 98, 1, 1); book(_h["abseta_jet2"], 104, 1, 1); book(_h["pT_jet1"], 110, 1, 1); book(_h["pT_jet2"], 116, 1, 1); } /// Perform the per-event analysis void analyze(Event const & event) { const auto& baseline_electrons = apply(event, "electrons").dressedLeptons(); const auto& baseline_muons = apply(event, "muons").dressedLeptons(); // Form all possible quadruplets passing hierarchical lepton pT cuts const auto quadruplets = getQuadruplets(baseline_electrons, baseline_muons); if (quadruplets.empty()) vetoEvent; // Select the best quadruplet, the one minimising the total distance from the Z pole mass auto const quadruplet = selectQuadruplet(quadruplets); // Event selection on the best quadruplet if (!passesDileptonMinimumMassRequirement(quadruplet)) vetoEvent; if (!passesLeptonDeltaRRequirements(quadruplet)) vetoEvent; if (!inRange(quadruplet.leadingDilepton.momentum().mass(), 66*GeV, 116*GeV)) vetoEvent; if (!inRange(quadruplet.subleadingDilepton.momentum().mass(), 66*GeV, 116*GeV)) vetoEvent; // Select jets Jets alljets = apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV); for (const DressedLepton& lep : quadruplet.leptonsSortedByPt) ifilter_discard(alljets, deltaRLess(lep, 0.4)); const Jets jets = alljets; const Jets centralJets = filterBy(jets, Cuts::abseta < 2.4); const Jets pt60Jets = filterBy(jets, Cuts::pT > 60*GeV); const auto& leadingDilepton = quadruplet.leadingDilepton.momentum(); const auto& subleadingDilepton = quadruplet.subleadingDilepton.momentum(); _h["pT_4l"]->fill((leadingDilepton + subleadingDilepton).pt()/GeV); _h["pT_leading_dilepton"]->fill(leadingDilepton.pt()/GeV); _h["pT_subleading_dilepton"]->fill(subleadingDilepton.pt()/GeV); _h["pT_lepton1"]->fill(quadruplet.leptonsSortedByPt.at(0).pt()/GeV); _h["pT_lepton2"]->fill(quadruplet.leptonsSortedByPt.at(1).pt()/GeV); _h["pT_lepton3"]->fill(quadruplet.leptonsSortedByPt.at(2).pt()/GeV); _h["pT_lepton4"]->fill(quadruplet.leptonsSortedByPt.at(3).pt()/GeV); _h["absy_4l"]->fill((leadingDilepton + subleadingDilepton).absrapidity()); _h["deltay_dileptons"]->fill(fabs(leadingDilepton.rapidity() - subleadingDilepton.rapidity())); _h["deltaphi_dileptons"]->fill(deltaPhi(leadingDilepton, subleadingDilepton)/pi); _h["N_jets"]->fill(jets.size()); _h["N_central_jets"]->fill(centralJets.size()); _h["N_jets60"]->fill(pt60Jets.size()); // If at least one jet present if (jets.empty()) vetoEvent; _h["scalarpTsum_jets"]->fill(sum(jets, pT, 0.)/GeV); _h["abseta_jet1"]->fill(jets.front().abseta()); _h["pT_jet1"]->fill(jets.front().pt()/GeV); // If at least two jets present if (jets.size() < 2) vetoEvent; _h["mass_dijet"]->fill((jets.at(0).mom() + jets.at(1).mom()).mass()/GeV); _h["deltay_dijet"]->fill(fabs(jets.at(0).rapidity() - jets.at(1).rapidity())); _h["abseta_jet2"]->fill(jets.at(1).abseta()); _h["pT_jet2"]->fill(jets.at(1).pt()/GeV); } /// Normalise histograms etc., after the run void finalize() { // Normalise histograms to cross section const double sf = crossSectionPerEvent() / femtobarn; for (map::iterator it = _h.begin(); it != _h.end(); ++it) { scale(it->second, sf); } } //@} private: /// @name Histograms //@{ map _h; static constexpr double Z_mass = 91.1876; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1625109); } diff --git a/analyses/pluginCMS/CMS_2012_I1111014.cc b/analyses/pluginCMS/CMS_2012_I1111014.cc --- a/analyses/pluginCMS/CMS_2012_I1111014.cc +++ b/analyses/pluginCMS/CMS_2012_I1111014.cc @@ -1,185 +1,185 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/JetShape.hh" namespace Rivet { /// @brief CMS jet shape analysis /// @author Andreas Hinzmann class CMS_2012_I1111014 : public Analysis { public: /// Constructor CMS_2012_I1111014() : Analysis("CMS_2012_I1111014") { } /// @name Analysis methods //@{ void init() { // Set up projections const FinalState fs(-5.0, 5.0); declare(fs, "FS"); FastJets fj5(fs, FastJets::ANTIKT, 0.5); declare(fj5, "Jets5"); FastJets fj7(fs, FastJets::ANTIKT, 0.7); declare(fj7, "Jets7"); // Specify pT bins _ptedges = {{ 20.,25.,30.,40.,50.,60.,70.,80.,90.,100.,110.,125.,140.,160.,180.,200.,225.,250.,300.,400.,500.,600.,1000. }}; _yedges = {{ 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0 }}; // Register a jet shape projection and histogram for each pT bin unsigned int histo_counter=1; for (size_t j = 0; j < 6; ++j) { for (size_t i = 0; i < 22; ++i) { if (i > 20 && j == 3) continue; if (i > 18 && j >= 4) continue; // Set up projections for each (pT,y) bin _jsnames_pT[i][j] = "JetShape" + to_str(i) + "_" + to_str(j); const JetShape jsp(fj7, 0.0, 0.7, 7, _ptedges[i], _ptedges[i+1], _yedges[j], _yedges[j+1], RAPIDITY); declare(jsp, _jsnames_pT[i][j]); // Book profile histograms for each (pT,y) bin book(_profhistRho_pT[i][j], histo_counter, 1, 1); histo_counter+=1; } } book(_profhistNch[0], 126, 1, 1); book(_profhistNch[1], 126, 1, 2); book(_profhistDr[0], 127, 1, 1); book(_profhistDr[1], 127, 1, 2); book(_profhistDeta, "TMP/Deta", refData(127,1,1)); book(_profhistDphi, "TMP/Dphi", refData(127,1,1)); book(_profhistAsym, "d128-x01-y01", true); } /// Do the analysis void analyze(const Event& evt) { // Get jets and require at least one to pass pT and y cuts Jets jets7 = apply(evt, "Jets7") .jetsByPt(Cuts::ptIn(_ptedges.front()*GeV, _ptedges.back()*GeV) && Cuts::absrap < 3.0); if(jets7.size()>2) jets7.resize(2); // Use only the first two jets MSG_DEBUG("Jet (R=0.7) multiplicity before cuts = " << jets7.size()); if (jets7.size() == 0) { MSG_DEBUG("No jets (R=0.7) found in required pT and rapidity range"); vetoEvent; } // Calculate and histogram jet shapes for (size_t jy = 0; jy < 6; ++jy) { for (size_t ipt = 0; ipt < 22; ++ipt) { if (ipt > 20 && jy == 3) continue; if (ipt > 18 && jy >= 4) continue; JetShape jsipt = apply(evt, _jsnames_pT[ipt][jy]); jsipt.calc(jets7); for (size_t ijet = 0; ijet < jsipt.numJets(); ++ijet) { for (size_t rbin = 0; rbin < jsipt.numBins(); ++rbin) { const double r_rho = jsipt.rBinMid(rbin); _profhistRho_pT[ipt][jy]->fill(r_rho, (1./0.1)*jsipt.diffJetShape(ijet, rbin)); } } } } // Get jets and require at least one to pass pT and y cuts Jets jets5 = apply(evt, "Jets5") .jetsByPt(Cuts::ptIn(50*GeV, 1000*GeV) && Cuts::absrap < 2.0); // Calculate and histogram charged jet shapes for (const Jet& jet : jets5) { double ncharge = 0; double eta=0; double phi=0; double sumpt=0; for (const Particle& p : jet.particles()) { - if ((p.pT() < 0.5) || (p.charge3()==0) || (abs(p.pdgId())==11) || (abs(p.pdgId())==13)) continue; + if ((p.pT() < 0.5) || (p.charge3()==0) || (abs(p.pid())==11) || (abs(p.pid())==13)) continue; ncharge++; sumpt+=p.pT(); eta+=p.pT()*p.eta(); phi+=p.pT()*mapAngleMPiToPi(p.phi()-jet.phi()); } if (jet.absrap()<1.0) { _profhistNch[0]->fill(jet.pT(), ncharge ); } else if (jet.absrap()<2.0) { _profhistNch[1]->fill(jet.pT(), ncharge ); } if (sumpt==0) continue; eta/=sumpt; phi/=sumpt; double deta=0; double dphi=0; for (const Particle& p : jet.particles()) { - if ((p.pT() < 0.5) || (p.charge3()==0) || (abs(p.pdgId())==11) || (abs(p.pdgId())==13)) continue; + if ((p.pT() < 0.5) || (p.charge3()==0) || (abs(p.pid())==11) || (abs(p.pid())==13)) continue; deta+=p.pT()*pow(p.eta()-eta,2); dphi+=p.pT()*pow(mapAngleMPiToPi(p.phi()-phi-jet.phi()),2); } deta/=sumpt; dphi/=sumpt; if ((dphi==0)||(deta==0)) continue; if (jet.absrap()<1.0) { _profhistDr[0]->fill(jet.pT(), deta+dphi ); _profhistDeta->fill(jet.pT(), deta ); _profhistDphi->fill(jet.pT(), dphi ); } else if (jet.absrap()<2.0) { _profhistDr[1]->fill(jet.pT(), deta+dphi ); } } } // Finalize void finalize() { for (unsigned int i = 0; i < _profhistAsym->numPoints(); ++i) { if((_profhistDeta->bin(i).numEntries()<2)||(_profhistDphi->bin(i).numEntries()<2)) continue; if((_profhistDeta->bin(i).mean()==0)||(_profhistDphi->bin(i).mean()==0)) continue; double mean_ratio=_profhistDeta->bin(i).mean() / _profhistDphi->bin(i).mean(); double mean_error=mean_ratio*sqrt(pow(_profhistDeta->bin(i).stdErr()/_profhistDeta->bin(i).mean(),2)+pow(_profhistDphi->bin(i).stdErr()/_profhistDphi->bin(i).mean(),2)); _profhistAsym->point(i).setY(mean_ratio,mean_error); } } //@} private: /// @name Analysis data //@{ /// Jet \f$ p_\perp\f$ bins. vector _ptedges; // This can't be a raw array if we want to initialise it non-painfully vector _yedges; /// JetShape projection name for each \f$p_\perp\f$ bin. string _jsnames_pT[22][6]; //@} /// @name Histograms //@{ Profile1DPtr _profhistRho_pT[22][6]; Profile1DPtr _profhistNch[2]; Profile1DPtr _profhistDr[2]; Profile1DPtr _profhistDeta; Profile1DPtr _profhistDphi; Scatter2DPtr _profhistAsym; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2012_I1111014); } diff --git a/analyses/pluginCMS/CMS_2015_I1370682.cc b/analyses/pluginCMS/CMS_2015_I1370682.cc --- a/analyses/pluginCMS/CMS_2015_I1370682.cc +++ b/analyses/pluginCMS/CMS_2015_I1370682.cc @@ -1,604 +1,604 @@ #include "Rivet/Analysis.hh" #include "Rivet/Math/LorentzTrans.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { namespace { //< only visible in this compilation unit /// @brief Pseudo top finder /// /// Find top quark in the particle level. /// The definition is based on the agreement at the LHC working group. class PseudoTop : public FinalState { public: /// @name Standard constructors and destructors. //@{ /// The default constructor. May specify the minimum and maximum /// pseudorapidity \f$ \eta \f$ and the min \f$ p_T \f$ (in GeV). PseudoTop(double lepR = 0.1, double lepMinPt = 20, double lepMaxEta = 2.4, double jetR = 0.4, double jetMinPt = 30, double jetMaxEta = 4.7) : FinalState(-DBL_MAX, DBL_MAX, 0*GeV), _lepR(lepR), _lepMinPt(lepMinPt), _lepMaxEta(lepMaxEta), _jetR(jetR), _jetMinPt(jetMinPt), _jetMaxEta(jetMaxEta) { setName("PseudoTop"); } enum TTbarMode {CH_NONE=-1, CH_FULLHADRON = 0, CH_SEMILEPTON, CH_FULLLEPTON}; enum DecayMode {CH_HADRON = 0, CH_MUON, CH_ELECTRON}; TTbarMode mode() const { if (!_isValid) return CH_NONE; if (_mode1 == CH_HADRON && _mode2 == CH_HADRON) return CH_FULLHADRON; else if ( _mode1 != CH_HADRON && _mode2 != CH_HADRON) return CH_FULLLEPTON; else return CH_SEMILEPTON; } DecayMode mode1() const {return _mode1;} DecayMode mode2() const {return _mode2;} /// Clone on the heap. virtual unique_ptr clone() const { return unique_ptr(new PseudoTop(*this)); } //@} public: Particle t1() const {return _t1;} Particle t2() const {return _t2;} Particle b1() const {return _b1;} Particle b2() const {return _b2;} ParticleVector wDecays1() const {return _wDecays1;} ParticleVector wDecays2() const {return _wDecays2;} Jets jets() const {return _jets;} Jets bjets() const {return _bjets;} Jets ljets() const {return _ljets;} protected: // Apply the projection to the event void project(const Event& e); // override; ///< @todo Re-enable when C++11 allowed void cleanup(std::map >& v, const bool doCrossCleanup=false) const; private: const double _lepR, _lepMinPt, _lepMaxEta; const double _jetR, _jetMinPt, _jetMaxEta; //constexpr ///< @todo Re-enable when C++11 allowed static double _tMass; // = 172.5*GeV; ///< @todo Re-enable when C++11 allowed //constexpr ///< @todo Re-enable when C++11 allowed static double _wMass; // = 80.4*GeV; ///< @todo Re-enable when C++11 allowed private: bool _isValid; DecayMode _mode1, _mode2; Particle _t1, _t2; Particle _b1, _b2; ParticleVector _wDecays1, _wDecays2; Jets _jets, _bjets, _ljets; }; // More implementation below the analysis code } /// Pseudo-top analysis from CMS class CMS_2015_I1370682 : public Analysis { public: CMS_2015_I1370682() : Analysis("CMS_2015_I1370682"), _applyCorrection(true), _doShapeOnly(true) { } void init() { declare(PseudoTop(0.1, 20, 2.4, 0.5, 30, 2.4), "ttbar"); // Lepton + Jet channel book(_hSL_topPt ,"d15-x01-y01"); // 1/sigma dsigma/dpt(top) book(_hSL_topPtTtbarSys ,"d16-x01-y01"); // 1/sigma dsigma/dpt*(top) book(_hSL_topY ,"d17-x01-y01"); // 1/sigma dsigma/dy(top) book(_hSL_ttbarDelPhi ,"d18-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar) book(_hSL_topPtLead ,"d19-x01-y01"); // 1/sigma dsigma/dpt(t1) book(_hSL_topPtSubLead ,"d20-x01-y01"); // 1/sigma dsigma/dpt(t2) book(_hSL_ttbarPt ,"d21-x01-y01"); // 1/sigma dsigma/dpt(ttbar) book(_hSL_ttbarY ,"d22-x01-y01"); // 1/sigma dsigma/dy(ttbar) book(_hSL_ttbarMass ,"d23-x01-y01"); // 1/sigma dsigma/dm(ttbar) // Dilepton channel book(_hDL_topPt ,"d24-x01-y01"); // 1/sigma dsigma/dpt(top) book(_hDL_topPtTtbarSys ,"d25-x01-y01"); // 1/sigma dsigma/dpt*(top) book(_hDL_topY ,"d26-x01-y01"); // 1/sigma dsigma/dy(top) book(_hDL_ttbarDelPhi ,"d27-x01-y01"); // 1/sigma dsigma/ddeltaphi(t,tbar) book(_hDL_topPtLead ,"d28-x01-y01"); // 1/sigma dsigma/dpt(t1) book(_hDL_topPtSubLead ,"d29-x01-y01"); // 1/sigma dsigma/dpt(t2) book(_hDL_ttbarPt ,"d30-x01-y01"); // 1/sigma dsigma/dpt(ttbar) book(_hDL_ttbarY ,"d31-x01-y01"); // 1/sigma dsigma/dy(ttbar) book(_hDL_ttbarMass ,"d32-x01-y01"); // 1/sigma dsigma/dm(ttbar) } void analyze(const Event& event) { // Get the ttbar candidate const PseudoTop& ttbar = apply(event, "ttbar"); if ( ttbar.mode() == PseudoTop::CH_NONE ) vetoEvent; const FourMomentum& t1P4 = ttbar.t1().momentum(); const FourMomentum& t2P4 = ttbar.t2().momentum(); const double pt1 = std::max(t1P4.pT(), t2P4.pT()); const double pt2 = std::min(t1P4.pT(), t2P4.pT()); const double dPhi = deltaPhi(t1P4, t2P4); const FourMomentum ttP4 = t1P4 + t2P4; const FourMomentum t1P4AtCM = LorentzTransform::mkFrameTransformFromBeta(ttP4.betaVec()).transform(t1P4); const double weight = 1.0; if ( ttbar.mode() == PseudoTop::CH_SEMILEPTON ) { const Particle lCand1 = ttbar.wDecays1()[0]; // w1 dau0 is the lepton in the PseudoTop if (lCand1.pT() < 33*GeV || lCand1.abseta() > 2.1) vetoEvent; _hSL_topPt->fill(t1P4.pT(), weight); _hSL_topPt->fill(t2P4.pT(), weight); _hSL_topPtTtbarSys->fill(t1P4AtCM.pT(), weight); _hSL_topY->fill(t1P4.rapidity(), weight); _hSL_topY->fill(t2P4.rapidity(), weight); _hSL_ttbarDelPhi->fill(dPhi, weight); _hSL_topPtLead->fill(pt1, weight); _hSL_topPtSubLead->fill(pt2, weight); _hSL_ttbarPt->fill(ttP4.pT(), weight); _hSL_ttbarY->fill(ttP4.rapidity(), weight); _hSL_ttbarMass->fill(ttP4.mass(), weight); } else if ( ttbar.mode() == PseudoTop::CH_FULLLEPTON ) { const Particle lCand1 = ttbar.wDecays1()[0]; // dau0 are the lepton in the PseudoTop const Particle lCand2 = ttbar.wDecays2()[0]; // dau0 are the lepton in the PseudoTop if (lCand1.pT() < 20*GeV || lCand1.abseta() > 2.4) vetoEvent; if (lCand2.pT() < 20*GeV || lCand2.abseta() > 2.4) vetoEvent; _hDL_topPt->fill(t1P4.pT(), weight); _hDL_topPt->fill(t2P4.pT(), weight); _hDL_topPtTtbarSys->fill(t1P4AtCM.pT(), weight); _hDL_topY->fill(t1P4.rapidity(), weight); _hDL_topY->fill(t2P4.rapidity(), weight); _hDL_ttbarDelPhi->fill(dPhi, weight); _hDL_topPtLead->fill(pt1, weight); _hDL_topPtSubLead->fill(pt2, weight); _hDL_ttbarPt->fill(ttP4.pT(), weight); _hDL_ttbarY->fill(ttP4.rapidity(), weight); _hDL_ttbarMass->fill(ttP4.mass(), weight); } } void finalize() { if ( _applyCorrection ) { // Correction functions for TOP-12-028 paper, (parton bin height)/(pseudotop bin height) const double ch15[] = { 5.473609, 4.941048, 4.173346, 3.391191, 2.785644, 2.371346, 2.194161, 2.197167, }; const double ch16[] = { 5.470905, 4.948201, 4.081982, 3.225532, 2.617519, 2.239217, 2.127878, 2.185918, }; const double ch17[] = { 10.003667, 4.546519, 3.828115, 3.601018, 3.522194, 3.524694, 3.600951, 3.808553, 4.531891, 9.995370, }; const double ch18[] = { 4.406683, 4.054041, 3.885393, 4.213646, }; const double ch19[] = { 6.182537, 5.257703, 4.422280, 3.568402, 2.889408, 2.415878, 2.189974, 2.173210, }; const double ch20[] = { 5.199874, 4.693318, 3.902882, 3.143785, 2.607877, 2.280189, 2.204124, 2.260829, }; const double ch21[] = { 6.053523, 3.777506, 3.562251, 3.601356, 3.569347, 3.410472, }; const double ch22[] = { 11.932351, 4.803773, 3.782709, 3.390775, 3.226806, 3.218982, 3.382678, 3.773653, 4.788191, 11.905338, }; const double ch23[] = { 7.145255, 5.637595, 4.049882, 3.025917, 2.326430, 1.773824, 1.235329, }; const double ch24[] = { 2.268193, 2.372063, 2.323975, 2.034655, 1.736793, }; const double ch25[] = { 2.231852, 2.383086, 2.341894, 2.031318, 1.729672, 1.486993, }; const double ch26[] = { 3.993526, 2.308249, 2.075136, 2.038297, 2.036302, 2.078270, 2.295817, 4.017713, }; const double ch27[] = { 2.205978, 2.175010, 2.215376, 2.473144, }; const double ch28[] = { 2.321077, 2.371895, 2.338871, 2.057821, 1.755382, }; const double ch29[] = { 2.222707, 2.372591, 2.301688, 1.991162, 1.695343, }; const double ch30[] = { 2.599677, 2.026855, 2.138620, 2.229553, }; const double ch31[] = { 5.791779, 2.636219, 2.103642, 1.967198, 1.962168, 2.096514, 2.641189, 5.780828, }; const double ch32[] = { 2.006685, 2.545525, 2.477745, 2.335747, 2.194226, 2.076500, }; applyCorrection(_hSL_topPt, ch15); applyCorrection(_hSL_topPtTtbarSys, ch16); applyCorrection(_hSL_topY, ch17); applyCorrection(_hSL_ttbarDelPhi, ch18); applyCorrection(_hSL_topPtLead, ch19); applyCorrection(_hSL_topPtSubLead, ch20); applyCorrection(_hSL_ttbarPt, ch21); applyCorrection(_hSL_ttbarY, ch22); applyCorrection(_hSL_ttbarMass, ch23); applyCorrection(_hDL_topPt, ch24); applyCorrection(_hDL_topPtTtbarSys, ch25); applyCorrection(_hDL_topY, ch26); applyCorrection(_hDL_ttbarDelPhi, ch27); applyCorrection(_hDL_topPtLead, ch28); applyCorrection(_hDL_topPtSubLead, ch29); applyCorrection(_hDL_ttbarPt, ch30); applyCorrection(_hDL_ttbarY, ch31); applyCorrection(_hDL_ttbarMass, ch32); } if ( _doShapeOnly ) { normalize(_hSL_topPt ); normalize(_hSL_topPtTtbarSys); normalize(_hSL_topY ); normalize(_hSL_ttbarDelPhi ); normalize(_hSL_topPtLead ); normalize(_hSL_topPtSubLead ); normalize(_hSL_ttbarPt ); normalize(_hSL_ttbarY ); normalize(_hSL_ttbarMass ); normalize(_hDL_topPt ); normalize(_hDL_topPtTtbarSys); normalize(_hDL_topY ); normalize(_hDL_ttbarDelPhi ); normalize(_hDL_topPtLead ); normalize(_hDL_topPtSubLead ); normalize(_hDL_ttbarPt ); normalize(_hDL_ttbarY ); normalize(_hDL_ttbarMass ); } else { const double s = 1./sumOfWeights(); scale(_hSL_topPt , s); scale(_hSL_topPtTtbarSys, s); scale(_hSL_topY , s); scale(_hSL_ttbarDelPhi , s); scale(_hSL_topPtLead , s); scale(_hSL_topPtSubLead , s); scale(_hSL_ttbarPt , s); scale(_hSL_ttbarY , s); scale(_hSL_ttbarMass , s); scale(_hDL_topPt , s); scale(_hDL_topPtTtbarSys, s); scale(_hDL_topY , s); scale(_hDL_ttbarDelPhi , s); scale(_hDL_topPtLead , s); scale(_hDL_topPtSubLead , s); scale(_hDL_ttbarPt , s); scale(_hDL_ttbarY , s); scale(_hDL_ttbarMass , s); } } void applyCorrection(Histo1DPtr h, const double* cf) { vector& bins = h->bins(); for (size_t i=0, n=bins.size(); i >& v, const bool doCrossCleanup) const { vector >::iterator> toErase; set usedLeg1, usedLeg2; if ( !doCrossCleanup ) { /// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) { for (map >::iterator key = v.begin(); key != v.end(); ++key) { const size_t leg1 = key->second.first; const size_t leg2 = key->second.second; if (usedLeg1.find(leg1) == usedLeg1.end() and usedLeg2.find(leg2) == usedLeg2.end()) { usedLeg1.insert(leg1); usedLeg2.insert(leg2); } else { toErase.push_back(key); } } } else { /// @todo Reinstate when C++11 allowed: for (auto key = v.begin(); key != v.end(); ++key) { for (map >::iterator key = v.begin(); key != v.end(); ++key) { const size_t leg1 = key->second.first; const size_t leg2 = key->second.second; if (usedLeg1.find(leg1) == usedLeg1.end() and usedLeg1.find(leg2) == usedLeg1.end()) { usedLeg1.insert(leg1); usedLeg1.insert(leg2); } else { toErase.push_back(key); } } } /// @todo Reinstate when C++11 allowed: for (auto& key : toErase) v.erase(key); for (size_t i = 0; i < toErase.size(); ++i) v.erase(toErase[i]); } void PseudoTop::project(const Event& e) { // Leptons : do the lepton clustering anti-kt R=0.1 using stable photons and leptons not from hadron decay // Neutrinos : neutrinos not from hadron decay // MET : vector sum of all invisible particles in x-y plane // Jets : anti-kt R=0.4 using all particles excluding neutrinos and particles used in lepton clustering // add ghost B hadrons during the jet clustering to identify B jets. // W->lv : dressed lepton and neutrino pairs // W->jj : light flavored dijet // W candidate : select lv or jj pairs which minimise |mW1-80.4|+|mW2-80.4| // lepton-neutrino pair will be selected with higher priority // t->Wb : W candidate + b jet // t candidate : select Wb pairs which minimise |mtop1-172.5|+|mtop2-172.5| _isValid = false; _theParticles.clear(); _wDecays1.clear(); _wDecays2.clear(); _jets.clear(); _bjets.clear(); _ljets.clear(); _mode1 = _mode2 = CH_HADRON; // Collect final state particles Particles pForLep, pForJet; Particles neutrinos; // Prompt neutrinos /// @todo Avoid this unsafe jump into HepMC -- all this can be done properly via VisibleFS and HeavyHadrons projections for (const GenParticle* p : Rivet::particles(e.genEvent())) { const int status = p->status(); - const int pdgId = p->pdg_id(); + const int pid = p->pdg_id(); if (status == 1) { Particle rp = *p; - if (!PID::isHadron(pdgId) && !rp.fromHadron()) { + if (!PID::isHadron(pid) && !rp.fromHadron()) { // Collect particles not from hadron decay if (rp.isNeutrino()) { // Prompt neutrinos are kept in separate collection neutrinos.push_back(rp); - } else if (pdgId == 22 || rp.isLepton()) { + } else if (pid == 22 || rp.isLepton()) { // Leptons and photons for the dressing pForLep.push_back(rp); } } else if (!rp.isNeutrino()) { // Use all particles from hadron decay pForJet.push_back(rp); } - } else if (PID::isHadron(pdgId) && PID::hasBottom(pdgId)) { + } else if (PID::isHadron(pid) && PID::hasBottom(pid)) { // NOTE: Consider B hadrons with pT > 5GeV - not in CMS proposal //if ( p->momentum().perp() < 5 ) continue; // Do unstable particles, to be used in the ghost B clustering // Use last B hadrons only bool isLast = true; for (const GenParticlePtr pp : Rivet::particles(p->end_vertex(), HepMC::children)) { if (PID::hasBottom(pp->pdg_id())) { isLast = false; break; } } if (!isLast) continue; // Rescale momentum by 10^-20 - Particle ghost(pdgId, FourMomentum(p->momentum())*1e-20/p->momentum().rho()); + Particle ghost(pid, FourMomentum(p->momentum())*1e-20/p->momentum().rho()); pForJet.push_back(ghost); } } // Start object building from trivial thing - prompt neutrinos sortByPt(neutrinos); // Proceed to lepton dressing FastJets fjLep(FinalState(), FastJets::ANTIKT, _lepR); fjLep.calc(pForLep); Jets leptons; vector leptonsId; set dressedIdxs; for (const Jet& lep : fjLep.jetsByPt(_lepMinPt)) { if (lep.abseta() > _lepMaxEta) continue; double leadingPt = -1; int leptonId = 0; for (const Particle& p : lep.particles()) { /// @warning Barcodes aren't future-proof in HepMC dressedIdxs.insert(p.genParticle()->barcode()); if (p.isLepton() && p.pT() > leadingPt) { leadingPt = p.pT(); leptonId = p.pid(); } } if (leptonId == 0) continue; leptons.push_back(lep); leptonsId.push_back(leptonId); } // Re-use particles not used in lepton dressing for (const Particle& rp : pForLep) { /// @warning Barcodes aren't future-proof in HepMC const int barcode = rp.genParticle()->barcode(); // Skip if the particle is used in dressing if (dressedIdxs.find(barcode) != dressedIdxs.end()) continue; // Put back to be used in jet clustering pForJet.push_back(rp); } // Then do the jet clustering FastJets fjJet(FinalState(), FastJets::ANTIKT, _jetR); //fjJet.useInvisibles(); // NOTE: CMS proposal to remove neutrinos (AB: wouldn't work anyway, since they were excluded from clustering inputs) fjJet.calc(pForJet); for (const Jet& jet : fjJet.jetsByPt(_jetMinPt)) { if (jet.abseta() > _jetMaxEta) continue; _jets.push_back(jet); bool isBJet = false; for (const Particle& rp : jet.particles()) { - if (PID::hasBottom(rp.pdgId())) { + if (PID::hasBottom(rp.pid())) { isBJet = true; break; } } if ( isBJet ) _bjets.push_back(jet); else _ljets.push_back(jet); } // Every building blocks are ready. Continue to pseudo-W and pseudo-top combination if (_bjets.size() < 2) return; // Ignore single top for now map > wLepCandIdxs; map > wHadCandIdxs; // Collect leptonic-decaying W's for (size_t iLep = 0, nLep = leptons.size(); iLep < nLep; ++iLep) { const Jet& lep = leptons.at(iLep); for (size_t iNu = 0, nNu = neutrinos.size(); iNu < nNu; ++iNu) { const Particle& nu = neutrinos.at(iNu); const double m = (lep.momentum()+nu.momentum()).mass(); const double dm = std::abs(m-_wMass); wLepCandIdxs[dm] = make_pair(iLep, iNu); } } // Continue to hadronic decaying W's for (size_t i = 0, nLjet = _ljets.size(); i < nLjet; ++i) { const Jet& ljet1 = _ljets[i]; for (size_t j = i+1; j < nLjet; ++j) { const Jet& ljet2 = _ljets[j]; const double m = (ljet1.momentum()+ljet2.momentum()).mass(); const double dm = std::abs(m-_wMass); wHadCandIdxs[dm] = make_pair(i, j); } } // Cleanup W candidate, choose pairs with minimum dm if they share decay products cleanup(wLepCandIdxs); cleanup(wHadCandIdxs, true); const size_t nWLepCand = wLepCandIdxs.size(); const size_t nWHadCand = wHadCandIdxs.size(); if (nWLepCand + nWHadCand < 2) return; // We skip single top int w1Q = 1, w2Q = -1; int w1dau1Id = 1, w2dau1Id = -1; FourMomentum w1dau1LVec, w1dau2LVec; FourMomentum w2dau1LVec, w2dau2LVec; if (nWLepCand == 0) { // Full hadronic case const pair& idPair1 = wHadCandIdxs.begin()->second; const pair& idPair2 = (++wHadCandIdxs.begin())->second; ///< @todo Reinstate std::next const Jet& w1dau1 = _ljets[idPair1.first]; const Jet& w1dau2 = _ljets[idPair1.second]; const Jet& w2dau1 = _ljets[idPair2.first]; const Jet& w2dau2 = _ljets[idPair2.second]; w1dau1LVec = w1dau1.momentum(); w1dau2LVec = w1dau2.momentum(); w2dau1LVec = w2dau1.momentum(); w2dau2LVec = w2dau2.momentum(); } else if (nWLepCand == 1) { // Semi-leptonic case const pair& idPair1 = wLepCandIdxs.begin()->second; const pair& idPair2 = wHadCandIdxs.begin()->second; const Jet& w1dau1 = leptons[idPair1.first]; const Particle& w1dau2 = neutrinos[idPair1.second]; const Jet& w2dau1 = _ljets[idPair2.first]; const Jet& w2dau2 = _ljets[idPair2.second]; w1dau1LVec = w1dau1.momentum(); w1dau2LVec = w1dau2.momentum(); w2dau1LVec = w2dau1.momentum(); w2dau2LVec = w2dau2.momentum(); w1dau1Id = leptonsId[idPair1.first]; w1Q = w1dau1Id > 0 ? -1 : 1; w2Q = -w1Q; switch (w1dau1Id) { case 13: case -13: _mode1 = CH_MUON; break; case 11: case -11: _mode1 = CH_ELECTRON; break; } } else { // Full leptonic case const pair& idPair1 = wLepCandIdxs.begin()->second; const pair& idPair2 = (++wLepCandIdxs.begin())->second; ///< @todo Reinstate std::next const Jet& w1dau1 = leptons[idPair1.first]; const Particle& w1dau2 = neutrinos[idPair1.second]; const Jet& w2dau1 = leptons[idPair2.first]; const Particle& w2dau2 = neutrinos[idPair2.second]; w1dau1LVec = w1dau1.momentum(); w1dau2LVec = w1dau2.momentum(); w2dau1LVec = w2dau1.momentum(); w2dau2LVec = w2dau2.momentum(); w1dau1Id = leptonsId[idPair1.first]; w2dau1Id = leptonsId[idPair2.first]; w1Q = w1dau1Id > 0 ? -1 : 1; w2Q = w2dau1Id > 0 ? -1 : 1; switch (w1dau1Id) { case 13: case -13: _mode1 = CH_MUON; break; case 11: case -11: _mode1 = CH_ELECTRON; break; } switch (w2dau1Id) { case 13: case -13: _mode2 = CH_MUON; break; case 11: case -11: _mode2 = CH_ELECTRON; break; } } const FourMomentum w1LVec = w1dau1LVec+w1dau2LVec; const FourMomentum w2LVec = w2dau1LVec+w2dau2LVec; // Combine b jets double sumDm = 1e9; FourMomentum b1LVec, b2LVec; for (size_t i = 0, n = _bjets.size(); i < n; ++i) { const Jet& bjet1 = _bjets[i]; const double mtop1 = (w1LVec+bjet1.momentum()).mass(); const double dmtop1 = std::abs(mtop1-_tMass); for (size_t j=0; j= 1e9) return; // Failed to make top, but this should not happen. const FourMomentum t1LVec = w1LVec + b1LVec; const FourMomentum t2LVec = w2LVec + b2LVec; // Put all of them into candidate collection _t1 = Particle(w1Q*6, t1LVec); _b1 = Particle(w1Q*5, b1LVec); _wDecays1.push_back(Particle(w1dau1Id, w1dau1LVec)); _wDecays1.push_back(Particle(-w1dau1Id+w1Q, w1dau2LVec)); _t2 = Particle(w2Q*6, t2LVec); _b2 = Particle(w2Q*5, b2LVec); _wDecays2.push_back(Particle(w2dau1Id, w2dau1LVec)); _wDecays2.push_back(Particle(-w2dau1Id+w2Q, w2dau2LVec)); _isValid = true; } } } diff --git a/analyses/pluginCMS/CMS_2016_I1413748.cc b/analyses/pluginCMS/CMS_2016_I1413748.cc --- a/analyses/pluginCMS/CMS_2016_I1413748.cc +++ b/analyses/pluginCMS/CMS_2016_I1413748.cc @@ -1,328 +1,328 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/PartonicTops.hh" namespace Rivet { /// CMS 8 TeV dilepton channel ttbar spin correlations and polarisation analysis class CMS_2016_I1413748 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1413748); /// Book histograms and initialise projections void init() { // Complete final state FinalState fs(-DBL_MAX, DBL_MAX, 0*GeV); // Projection for dressed electrons and muons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); declare(electrons, "Electrons"); DressedLeptons dressed_electrons(photons, electrons, 0.1); declare(dressed_electrons, "DressedElectrons"); IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); declare(muons, "Muons"); DressedLeptons dressed_muons(photons, muons, 0.1); declare(dressed_muons, "DressedMuons"); // Parton-level top quarks declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicPartonTops"); // Booking of histograms // This histogram is independent of the parton-level information, and is an addition to the original analysis. // It is compared to the same data as the parton-level delta_phi histogram d02-x01-y01. book(_h_dphidressedleptons, "d00-x01-y01", _bins_dphi); // The remaining histos use parton-level information book(_h_dphi, "d02-x01-y01", _bins_dphi); book(_h_cos_opening_angle, "d05-x01-y01", _bins_cos_opening_angle); book(_h_c1c2, "d08-x01-y01", _bins_c1c2); book(_h_lep_costheta, "d11-x01-y01", _bins_lep_costheta); book(_h_lep_costheta_CPV, "d14-x01-y01", _bins_lep_costheta_CPV); // 2D histos book(_h_dphi_var[0], "d20-x01-y01", _bins_dphi, _bins_tt_mass); book(_h_cos_opening_angle_var[0], "d26-x01-y01", _bins_cos_opening_angle, _bins_tt_mass); book(_h_c1c2_var[0], "d32-x01-y01", _bins_c1c2, _bins_tt_mass); book(_h_lep_costheta_var[0], "d38-x01-y01", _bins_lep_costheta, _bins_tt_mass); book(_h_lep_costheta_CPV_var[0], "d44-x01-y01", _bins_lep_costheta_CPV, _bins_tt_mass); book(_h_dphi_var[1], "d50-x01-y01", _bins_dphi, _bins_tt_pT); book(_h_cos_opening_angle_var[1], "d56-x01-y01", _bins_cos_opening_angle, _bins_tt_pT); book(_h_c1c2_var[1], "d62-x01-y01", _bins_c1c2, _bins_tt_pT); book(_h_lep_costheta_var[1], "d68-x01-y01", _bins_lep_costheta, _bins_tt_pT); book(_h_lep_costheta_CPV_var[1], "d74-x01-y01", _bins_lep_costheta_CPV, _bins_tt_pT); book(_h_dphi_var[2], "d80-x01-y01", _bins_dphi, _bins_tt_absrapidity); book(_h_cos_opening_angle_var[2], "d86-x01-y01", _bins_cos_opening_angle, _bins_tt_absrapidity); book(_h_c1c2_var[2], "d92-x01-y01", _bins_c1c2, _bins_tt_absrapidity); book(_h_lep_costheta_var[2], "d98-x01-y01", _bins_lep_costheta, _bins_tt_absrapidity); book(_h_lep_costheta_CPV_var[2], "d104-x01-y01", _bins_lep_costheta_CPV, _bins_tt_absrapidity); // Profile histos for asymmetries book(_h_dphi_profile[0], "d17-x01-y01", _bins_tt_mass); book(_h_cos_opening_angle_profile[0], "d23-x01-y01", _bins_tt_mass); book(_h_c1c2_profile[0], "d29-x01-y01", _bins_tt_mass); book(_h_lep_costheta_profile[0], "d35-x01-y01", _bins_tt_mass); book(_h_lep_costheta_CPV_profile[0], "d41-x01-y01", _bins_tt_mass); book(_h_dphi_profile[1], "d47-x01-y01", _bins_tt_pT); book(_h_cos_opening_angle_profile[1], "d53-x01-y01", _bins_tt_pT); book(_h_c1c2_profile[1], "d59-x01-y01", _bins_tt_pT); book(_h_lep_costheta_profile[1], "d65-x01-y01", _bins_tt_pT); book(_h_lep_costheta_CPV_profile[1], "d71-x01-y01", _bins_tt_pT); book(_h_dphi_profile[2], "d77-x01-y01", _bins_tt_absrapidity); book(_h_cos_opening_angle_profile[2], "d83-x01-y01", _bins_tt_absrapidity); book(_h_c1c2_profile[2], "d89-x01-y01", _bins_tt_absrapidity); book(_h_lep_costheta_profile[2], "d95-x01-y01", _bins_tt_absrapidity); book(_h_lep_costheta_CPV_profile[2], "d101-x01-y01", _bins_tt_absrapidity); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = 1.0; // Use particle-level leptons for the first histogram const DressedLeptons& dressed_electrons = applyProjection(event, "DressedElectrons"); const DressedLeptons& dressed_muons = applyProjection(event, "DressedMuons"); const vector dressedels = dressed_electrons.dressedLeptons(); const vector dressedmus = dressed_muons.dressedLeptons(); const size_t ndressedel = dressedels.size(); const size_t ndressedmu = dressedmus.size(); // For the particle-level histogram, require exactly one electron and exactly one muon, to select // the ttbar->emu channel. Note this means ttbar->emu events with additional PromptFinalState // dilepton pairs from the shower are vetoed - for PYTHIA8, this affects ~0.5% of events, so the // effect is well below the level of sensitivity of the measured distribution. if ( ndressedel == 1 && ndressedmu == 1 ) { const int electrontouse = 0, muontouse = 0; // Opposite-charge leptons only if ( sameSign(dressedels[electrontouse],dressedmus[muontouse]) ) { MSG_INFO("Error, e and mu have same charge, skipping event"); } else { //Get the four-momenta of the positively- and negatively-charged leptons FourMomentum lepPlus = dressedels[electrontouse].charge() > 0 ? dressedels[electrontouse] : dressedmus[muontouse]; FourMomentum lepMinus = dressedels[electrontouse].charge() > 0 ? dressedmus[muontouse] : dressedels[electrontouse]; // Now calculate the variable double dphi_temp = deltaPhi(lepPlus,lepMinus); fillWithUFOF( _h_dphidressedleptons, dphi_temp, weight ); } } // The remaining variables use parton-level information. // Get the leptonically decaying tops const Particles& leptonicpartontops = apply(event, "LeptonicPartonTops").particlesByPt(); Particles chargedleptons; unsigned int ntrueleptonictops = 0; bool oppositesign = false; if ( leptonicpartontops.size() == 2 ) { for (size_t k = 0; k < leptonicpartontops.size(); ++k) { // Get the lepton const Particle lepTop = leptonicpartontops[k]; const auto isPromptChargedLepton = [](const Particle& p){return (isChargedLepton(p) && isPrompt(p, false, false));}; Particles lepton_candidates = lepTop.allDescendants(firstParticleWith(isPromptChargedLepton), false); if ( lepton_candidates.size() < 1 ) MSG_WARNING("error, PartonicTops::DecayMode::E_MU top quark had no daughter lepton candidate, skipping event."); // In some cases there is no lepton from the W decay but only leptons from the decay of a radiated gamma. // These hadronic PartonicTops are currently being mistakenly selected by PartonicTops::DecayMode::E_MU (as of April 2017), and need to be rejected. // PartonicTops::DecayMode::E_MU is being fixed in Rivet, and when it is the veto below should do nothing. /// @todo Should no longer be necessary -- remove bool istrueleptonictop = false; for (size_t i = 0; i < lepton_candidates.size(); ++i) { const Particle& lepton_candidate = lepton_candidates[i]; if ( lepton_candidate.hasParent(PID::PHOTON) ) { MSG_DEBUG("Found gamma parent, top: " << k+1 << " of " << leptonicpartontops.size() << " , lepton: " << i+1 << " of " << lepton_candidates.size()); continue; } if ( !istrueleptonictop && sameSign(lepTop,lepton_candidate) ) { chargedleptons.push_back(lepton_candidate); istrueleptonictop = true; } else MSG_WARNING("Found extra prompt charged lepton from top decay (and without gamma parent), ignoring it."); } if ( istrueleptonictop ) ++ntrueleptonictops; } } if ( ntrueleptonictops == 2 ) { oppositesign = !( sameSign(chargedleptons[0],chargedleptons[1]) ); if ( !oppositesign ) MSG_WARNING("error, same charge tops, skipping event."); } if ( ntrueleptonictops == 2 && oppositesign ) { // Get the four-momenta of the positively- and negatively-charged leptons FourMomentum lepPlus = chargedleptons[0].charge() > 0 ? chargedleptons[0] : chargedleptons[1]; FourMomentum lepMinus = chargedleptons[0].charge() > 0 ? chargedleptons[1] : chargedleptons[0]; const double dphi_temp = deltaPhi(lepPlus,lepMinus); // Get the four-momenta of the positively- and negatively-charged tops - FourMomentum topPlus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[0] : leptonicpartontops[1]; - FourMomentum topMinus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[1] : leptonicpartontops[0]; + FourMomentum topPlus_p4 = leptonicpartontops[0].pid() > 0 ? leptonicpartontops[0] : leptonicpartontops[1]; + FourMomentum topMinus_p4 = leptonicpartontops[0].pid() > 0 ? leptonicpartontops[1] : leptonicpartontops[0]; const FourMomentum ttbar_p4 = topPlus_p4 + topMinus_p4; const double tt_mass_temp = ttbar_p4.mass(); const double tt_absrapidity_temp = ttbar_p4.absrapidity(); const double tt_pT_temp = ttbar_p4.pT(); // Lorentz transformations to calculate the spin observables in the helicity basis // Transform everything to the ttbar CM frame LorentzTransform ttCM; ttCM.setBetaVec(-ttbar_p4.betaVec()); topPlus_p4 = ttCM.transform(topPlus_p4); topMinus_p4 = ttCM.transform(topMinus_p4); lepPlus = ttCM.transform(lepPlus); lepMinus = ttCM.transform(lepMinus); // Now boost the leptons to their parent top CM frames LorentzTransform topPlus, topMinus; topPlus.setBetaVec(-topPlus_p4.betaVec()); topMinus.setBetaVec(-topMinus_p4.betaVec()); lepPlus = topPlus.transform(lepPlus); lepMinus = topMinus.transform(lepMinus); const double lepPlus_costheta_temp = lepPlus.vector3().dot(topPlus_p4.vector3()) / (lepPlus.vector3().mod() * topPlus_p4.vector3().mod()); const double lepMinus_costheta_temp = lepMinus.vector3().dot(topMinus_p4.vector3()) / (lepMinus.vector3().mod() * topMinus_p4.vector3().mod()); const double c1c2_temp = lepPlus_costheta_temp * lepMinus_costheta_temp; const double cos_opening_angle_temp = lepPlus.vector3().dot(lepMinus.vector3()) / (lepPlus.vector3().mod() * lepMinus.vector3().mod()); // Fill parton-level histos fillWithUFOF( _h_dphi, dphi_temp, weight ); fillWithUFOF( _h_cos_opening_angle, cos_opening_angle_temp, weight ); fillWithUFOF( _h_c1c2, c1c2_temp, weight ); fillWithUFOF( _h_lep_costheta, lepPlus_costheta_temp, weight ); fillWithUFOF( _h_lep_costheta, lepMinus_costheta_temp, weight ); fillWithUFOF( _h_lep_costheta_CPV, lepPlus_costheta_temp, weight ); fillWithUFOF( _h_lep_costheta_CPV, -lepMinus_costheta_temp, weight ); // Now fill the same variables in the 2D and profile histos vs ttbar invariant mass, pT, and absolute rapidity for (int i_var = 0; i_var < 3; ++i_var) { double var; if ( i_var == 0 ) { var = tt_mass_temp; } else if ( i_var == 1 ) { var = tt_pT_temp; } else { var = tt_absrapidity_temp; } fillWithUFOF( _h_dphi_var[i_var], dphi_temp, var, weight ); fillWithUFOF( _h_cos_opening_angle_var[i_var], cos_opening_angle_temp, var, weight ); fillWithUFOF( _h_c1c2_var[i_var], c1c2_temp, var, weight ); fillWithUFOF( _h_lep_costheta_var[i_var], lepPlus_costheta_temp, var, weight ); fillWithUFOF( _h_lep_costheta_var[i_var], lepMinus_costheta_temp, var, weight ); fillWithUFOF( _h_lep_costheta_CPV_var[i_var], lepPlus_costheta_temp, var, weight ); fillWithUFOF( _h_lep_costheta_CPV_var[i_var], -lepMinus_costheta_temp, var, weight ); fillWithUFOF( _h_dphi_profile[i_var], dphi_temp, var, weight, (_h_dphi->xMax() + _h_dphi->xMin())/2. ); fillWithUFOF( _h_cos_opening_angle_profile[i_var], cos_opening_angle_temp, var, weight, (_h_cos_opening_angle->xMax() + _h_cos_opening_angle->xMin())/2. ); fillWithUFOF( _h_c1c2_profile[i_var], c1c2_temp, var, weight, (_h_c1c2->xMax() + _h_c1c2->xMin())/2. ); fillWithUFOF( _h_lep_costheta_profile[i_var], lepPlus_costheta_temp, var, weight, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. ); fillWithUFOF( _h_lep_costheta_profile[i_var], lepMinus_costheta_temp, var, weight, (_h_lep_costheta->xMax() + _h_lep_costheta->xMin())/2. ); fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], lepPlus_costheta_temp, var, weight, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. ); fillWithUFOF( _h_lep_costheta_CPV_profile[i_var], -lepMinus_costheta_temp, var, weight, (_h_lep_costheta_CPV->xMax() + _h_lep_costheta_CPV->xMin())/2. ); } } } /// Normalise histograms to unit area void finalize() { normalize(_h_dphidressedleptons); normalize(_h_dphi); normalize(_h_cos_opening_angle); normalize(_h_c1c2); normalize(_h_lep_costheta); normalize(_h_lep_costheta_CPV); for (int i_var = 0; i_var < 3; ++i_var) { normalize(_h_dphi_var[i_var]); normalize(_h_cos_opening_angle_var[i_var]); normalize(_h_c1c2_var[i_var]); normalize(_h_lep_costheta_var[i_var]); normalize(_h_lep_costheta_CPV_var[i_var]); } } private: Histo1DPtr _h_dphidressedleptons, _h_dphi, _h_lep_costheta, _h_lep_costheta_CPV, _h_c1c2, _h_cos_opening_angle; Histo2DPtr _h_dphi_var[3], _h_lep_costheta_var[3], _h_lep_costheta_CPV_var[3], _h_c1c2_var[3], _h_cos_opening_angle_var[3]; Profile1DPtr _h_dphi_profile[3], _h_lep_costheta_profile[3], _h_lep_costheta_CPV_profile[3], _h_c1c2_profile[3], _h_cos_opening_angle_profile[3]; const vector _bins_tt_mass = {300., 430., 530., 1200.}; const vector _bins_tt_pT = {0., 41., 92., 300.}; const vector _bins_tt_absrapidity = {0., 0.34, 0.75, 1.5}; const vector _bins_dphi = {0., 5.*M_PI/60., 10.*M_PI/60., 15.*M_PI/60., 20.*M_PI/60., 25.*M_PI/60., 30.*M_PI/60., 35.*M_PI/60., 40.*M_PI/60., 45.*M_PI/60., 50.*M_PI/60., 55.*M_PI/60., M_PI}; const vector _bins_lep_costheta = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.}; const vector _bins_lep_costheta_CPV = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.}; const vector _bins_c1c2 = {-1., -0.4, -10./60., 0., 10./60., 0.4, 1.}; const vector _bins_cos_opening_angle = {-1., -2./3., -1./3., 0., 1./3., 2./3., 1.}; void fillWithUFOF(Histo1DPtr h, double x, double w) { h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), w); } void fillWithUFOF(Histo2DPtr h, double x, double y, double w) { h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), std::max(std::min(y, h->yMax()-1e-9),h->yMin()+1e-9), w); } void fillWithUFOF(Profile1DPtr h, double x, double y, double w, double c) { h->fill(std::max(std::min(y, h->xMax()-1e-9),h->xMin()+1e-9), float(x > c) - float(x < c), w); } }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_I1413748); } diff --git a/analyses/pluginCMS/CMS_2016_I1430892.cc b/analyses/pluginCMS/CMS_2016_I1430892.cc --- a/analyses/pluginCMS/CMS_2016_I1430892.cc +++ b/analyses/pluginCMS/CMS_2016_I1430892.cc @@ -1,259 +1,259 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/PartonicTops.hh" namespace Rivet { /// CMS 8 TeV dilepton channel ttbar charge asymmetry analysis class CMS_2016_I1430892 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_I1430892); /// Book histograms and initialise projections void init() { // Complete final state FinalState fs(-DBL_MAX, DBL_MAX, 0*GeV); // Projection for dressed electrons and muons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); declare(electrons, "Electrons"); DressedLeptons dressed_electrons(photons, electrons, 0.1); declare(dressed_electrons, "DressedElectrons"); IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); declare(muons, "Muons"); DressedLeptons dressed_muons(photons, muons, 0.1); declare(dressed_muons, "DressedMuons"); // Parton-level top quarks declare(PartonicTops(PartonicTops::DecayMode::E_MU, false), "LeptonicPartonTops"); // Booking of histograms // This histogram is independent of the parton-level information, and is an // addition to the original analysis. It is compared to the same data as // the parton-level delta_abseta histogram d05-x01-y01. book(_h_dabsetadressedleptons, "d00-x01-y01", _bins_dabseta); // The remaining histos use parton-level information book(_h_dabseta, "d05-x01-y01", _bins_dabseta); book(_h_dabsrapidity, "d02-x01-y01", _bins_dabsrapidity); // 2D histos book(_h_dabsrapidity_var[0], "d11-x01-y01", _bins_dabsrapidity, _bins_tt_mass); book(_h_dabseta_var[0], "d17-x01-y01", _bins_dabseta, _bins_tt_mass); book(_h_dabsrapidity_var[1], "d23-x01-y01", _bins_dabsrapidity, _bins_tt_pT); book(_h_dabseta_var[1], "d29-x01-y01", _bins_dabseta, _bins_tt_pT); book(_h_dabsrapidity_var[2], "d35-x01-y01", _bins_dabsrapidity, _bins_tt_absrapidity); book(_h_dabseta_var[2], "d41-x01-y01", _bins_dabseta, _bins_tt_absrapidity); // Profile histos for asymmetries book(_h_dabsrapidity_profile[0], "d08-x01-y01", _bins_tt_mass); book(_h_dabseta_profile[0], "d14-x01-y01", _bins_tt_mass); book(_h_dabsrapidity_profile[1], "d20-x01-y01", _bins_tt_pT); book(_h_dabseta_profile[1], "d26-x01-y01", _bins_tt_pT); book(_h_dabsrapidity_profile[2], "d32-x01-y01", _bins_tt_absrapidity); book(_h_dabseta_profile[2], "d38-x01-y01", _bins_tt_absrapidity); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = 1.0; // Use particle-level leptons for the first histogram const DressedLeptons& dressed_electrons = applyProjection(event, "DressedElectrons"); const DressedLeptons& dressed_muons = applyProjection(event, "DressedMuons"); const vector dressedels = dressed_electrons.dressedLeptons(); const vector dressedmus = dressed_muons.dressedLeptons(); const size_t ndressedel = dressedels.size(); const size_t ndressedmu = dressedmus.size(); // For the particle-level histogram, require exactly one electron and exactly one muon, to select the ttbar->emu channel. // Note this means ttbar->emu events with additional PromptFinalState dilepton pairs from the shower are vetoed - for PYTHIA8, // this affects ~0.5% of events, so the effect is well below the level of sensitivity of the measured distribution. if ( ndressedel == 1 && ndressedmu == 1 ) { const int electrontouse = 0, muontouse = 0; // Opposite-charge leptons only if ( sameSign(dressedels[electrontouse], dressedmus[muontouse]) ) { MSG_INFO("Error, e and mu have same charge, skipping event"); } else { // Get the four-momenta of the positively- and negatively-charged leptons FourMomentum lepPlus = dressedels[electrontouse].charge() > 0 ? dressedels[electrontouse] : dressedmus[muontouse]; FourMomentum lepMinus = dressedels[electrontouse].charge() > 0 ? dressedmus[muontouse] : dressedels[electrontouse]; // Now calculate the variable double dabseta_temp = lepPlus.abseta() - lepMinus.abseta(); fillWithUFOF( _h_dabsetadressedleptons, dabseta_temp, weight ); } } // The remaining variables use parton-level information. // Get the leptonically decaying tops const Particles leptonicpartontops = apply(event, "LeptonicPartonTops").particlesByPt(); Particles chargedleptons; unsigned int ntrueleptonictops = 0; bool oppositesign = false; if ( leptonicpartontops.size() == 2 ) { for (size_t k = 0; k < leptonicpartontops.size(); ++k) { // Get the lepton const Particle lepTop = leptonicpartontops[k]; const auto isPromptChargedLepton = [](const Particle& p){return (isChargedLepton(p) && isPrompt(p, false, false));}; Particles lepton_candidates = lepTop.allDescendants(firstParticleWith(isPromptChargedLepton), false); if ( lepton_candidates.size() < 1 ) MSG_WARNING("error, PartonicTops::DecayMode::E_MU top quark had no daughter lepton candidate, skipping event."); // In some cases there is no lepton from the W decay but only leptons from the decay of a radiated gamma. // These hadronic PartonicTops are currently being mistakenly selected by PartonicTops::DecayMode::E_MU (as of April 2017), and need to be rejected. // PartonicTops::DecayMode::E_MU is being fixed in Rivet, and when it is the veto below should do nothing. /// @todo Should no longer be necessary -- remove bool istrueleptonictop = false; for (size_t i = 0; i < lepton_candidates.size(); ++i) { const Particle& lepton_candidate = lepton_candidates[i]; if ( lepton_candidate.hasParent(PID::PHOTON) ) { MSG_DEBUG("Found gamma parent, top: " << k+1 << " of " << leptonicpartontops.size() << " , lepton: " << i+1 << " of " << lepton_candidates.size()); continue; } if ( !istrueleptonictop && sameSign(lepTop,lepton_candidate) ) { chargedleptons.push_back(lepton_candidate); istrueleptonictop = true; } else MSG_WARNING("Error, found extra prompt charged lepton from top decay (and without gamma parent), ignoring it."); } if ( istrueleptonictop ) ++ntrueleptonictops; } } if ( ntrueleptonictops == 2 ) { oppositesign = !( sameSign(chargedleptons[0],chargedleptons[1]) ); if ( !oppositesign ) MSG_WARNING("error, same charge tops, skipping event."); } if ( ntrueleptonictops == 2 && oppositesign ) { // Get the four-momenta of the positively- and negatively-charged leptons const FourMomentum lepPlus = chargedleptons[0].charge() > 0 ? chargedleptons[0] : chargedleptons[1]; const FourMomentum lepMinus = chargedleptons[0].charge() > 0 ? chargedleptons[1] : chargedleptons[0]; const double dabseta_temp = lepPlus.abseta() - lepMinus.abseta(); // Get the four-momenta of the positively- and negatively-charged tops - const FourMomentum topPlus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[0] : leptonicpartontops[1]; - const FourMomentum topMinus_p4 = leptonicpartontops[0].pdgId() > 0 ? leptonicpartontops[1] : leptonicpartontops[0]; + const FourMomentum topPlus_p4 = leptonicpartontops[0].pid() > 0 ? leptonicpartontops[0] : leptonicpartontops[1]; + const FourMomentum topMinus_p4 = leptonicpartontops[0].pid() > 0 ? leptonicpartontops[1] : leptonicpartontops[0]; const FourMomentum ttbar_p4 = topPlus_p4 + topMinus_p4; const double tt_mass_temp = ttbar_p4.mass(); const double tt_absrapidity_temp = ttbar_p4.absrapidity(); const double tt_pT_temp = ttbar_p4.pT(); const double dabsrapidity_temp = topPlus_p4.absrapidity() - topMinus_p4.absrapidity(); // Fill parton-level histos fillWithUFOF( _h_dabseta, dabseta_temp, weight ); fillWithUFOF( _h_dabsrapidity, dabsrapidity_temp, weight ); // Now fill the same variables in the 2D and profile histos vs ttbar invariant mass, pT, and absolute rapidity for (int i_var = 0; i_var < 3; ++i_var) { double var; if ( i_var == 0 ) { var = tt_mass_temp; } else if ( i_var == 1 ) { var = tt_pT_temp; } else { var = tt_absrapidity_temp; } fillWithUFOF( _h_dabsrapidity_var[i_var], dabsrapidity_temp, var, weight ); fillWithUFOF( _h_dabseta_var[i_var], dabseta_temp, var, weight ); fillWithUFOF( _h_dabsrapidity_profile[i_var], dabsrapidity_temp, var, weight, (_h_dabsrapidity->xMax() + _h_dabsrapidity->xMin())/2. ); fillWithUFOF( _h_dabseta_profile[i_var], dabseta_temp, var, weight, (_h_dabseta->xMax() + _h_dabseta->xMin())/2. ); } } } /// Normalise histograms to unit area void finalize() { normalize(_h_dabsetadressedleptons); normalize(_h_dabseta); normalize(_h_dabsrapidity); for (int i_var = 0; i_var < 3; ++i_var) { normalize(_h_dabsrapidity_var[i_var]); normalize(_h_dabseta_var[i_var]); } } private: Histo1DPtr _h_dabsetadressedleptons, _h_dabseta, _h_dabsrapidity; Histo2DPtr _h_dabseta_var[3], _h_dabsrapidity_var[3]; Profile1DPtr _h_dabseta_profile[3], _h_dabsrapidity_profile[3]; const vector _bins_tt_mass = {300., 430., 530., 1200.}; const vector _bins_tt_pT = {0., 41., 92., 300.}; const vector _bins_tt_absrapidity = {0., 0.34, 0.75, 1.5}; const vector _bins_dabseta = { -2., -68./60., -48./60., -32./60., -20./60., -8./60., 0., 8./60., 20./60., 32./60., 48./60., 68./60., 2.}; const vector _bins_dabsrapidity = {-2., -44./60., -20./60., 0., 20./60., 44./60., 2.}; void fillWithUFOF(Histo1DPtr h, double x, double w) { h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), w); } void fillWithUFOF(Histo2DPtr h, double x, double y, double w) { h->fill(std::max(std::min(x, h->xMax()-1e-9),h->xMin()+1e-9), std::max(std::min(y, h->yMax()-1e-9),h->yMin()+1e-9), w); } void fillWithUFOF(Profile1DPtr h, double x, double y, double w, double c) { h->fill(std::max(std::min(y, h->xMax()-1e-9),h->xMin()+1e-9), float(x > c) - float(x < c), w); } }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_I1430892); } diff --git a/analyses/pluginCMS/CMS_2016_I1491950.cc b/analyses/pluginCMS/CMS_2016_I1491950.cc --- a/analyses/pluginCMS/CMS_2016_I1491950.cc +++ b/analyses/pluginCMS/CMS_2016_I1491950.cc @@ -1,500 +1,500 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/ParticleIdUtils.hh" namespace Rivet { namespace { //< only visible in this compilation unit /// @brief Special dressed lepton finder /// /// Find dressed leptons by clustering all leptons and photons class SpecialDressedLeptons : public FinalState { public: /// The default constructor. May specify cuts SpecialDressedLeptons(const FinalState& fs, const Cut& cut) : FinalState(cut) { setName("SpecialDressedLeptons"); IdentifiedFinalState ifs(fs); ifs.acceptIdPair(PID::PHOTON); ifs.acceptIdPair(PID::ELECTRON); ifs.acceptIdPair(PID::MUON); declare(ifs, "IFS"); declare(FastJets(ifs, FastJets::ANTIKT, 0.1), "LeptonJets"); } /// Clone on the heap. virtual unique_ptr clone() const { return unique_ptr(new SpecialDressedLeptons(*this)); } /// Retrieve the dressed leptons const vector& dressedLeptons() const { return _clusteredLeptons; } private: /// Container which stores the clustered lepton objects vector _clusteredLeptons; public: void project(const Event& e) { _theParticles.clear(); _clusteredLeptons.clear(); vector allClusteredLeptons; const Jets jets = applyProjection(e, "LeptonJets").jetsByPt(5.*GeV); for (const Jet& jet : jets) { Particle lepCand; for (const Particle& cand : jet.particles()) { - const int absPdgId = abs(cand.pdgId()); + const int absPdgId = abs(cand.pid()); if (absPdgId == PID::ELECTRON || absPdgId == PID::MUON) { if (cand.pt() > lepCand.pt()) lepCand = cand; } } //Central lepton must be the major component - if ((lepCand.pt() < jet.pt()/2.) || (lepCand.pdgId() == 0)) continue; + if ((lepCand.pt() < jet.pt()/2.) || (lepCand.pid() == 0)) continue; DressedLepton lepton = DressedLepton(lepCand); for (const Particle& cand : jet.particles()) { if (cand == lepCand) continue; if (cand.pid() != PID::PHOTON) continue; lepton.addPhoton(cand, true); } allClusteredLeptons.push_back(lepton); } for (const DressedLepton& lepton : allClusteredLeptons) { if (accept(lepton)) { _clusteredLeptons.push_back(lepton); _theParticles.push_back(lepton.constituentLepton()); _theParticles += lepton.constituentPhotons(); } } } }; } class CMS_2016_I1491950 : public Analysis { public: /// Constructor CMS_2016_I1491950() : Analysis("CMS_2016_I1491950") { } /// Book histograms and initialise projections before the run void init() { FinalState fs(Cuts::pT > 0. && Cuts::abseta < 6.); PromptFinalState prompt_fs(fs); prompt_fs.acceptMuonDecays(true); prompt_fs.acceptTauDecays(true); // Projection for dressed electrons and muons Cut leptonCuts = Cuts::abseta < 2.5 and Cuts::pt > 30.*GeV; SpecialDressedLeptons dressedleptons(prompt_fs, leptonCuts); declare(dressedleptons, "DressedLeptons"); // Neutrinos IdentifiedFinalState neutrinos(prompt_fs); neutrinos.acceptNeutrinos(); declare(neutrinos, "Neutrinos"); // Projection for jets VetoedFinalState fsForJets(fs); fsForJets.addVetoOnThisFinalState(dressedleptons); fsForJets.addVetoOnThisFinalState(neutrinos); declare(FastJets(fsForJets, FastJets::ANTIKT, 0.4, JetAlg::Muons::DECAY, JetAlg::Invisibles::DECAY), "Jets"); //book hists book(_hist_thadpt, "d01-x02-y01"); book(_hist_thady, "d03-x02-y01"); book(_hist_tleppt, "d05-x02-y01"); book(_hist_tlepy, "d07-x02-y01"); book(_hist_ttpt, "d09-x02-y01"); book(_hist_tty, "d13-x02-y01"); book(_hist_ttm, "d11-x02-y01"); book(_hist_njet, "d15-x02-y01"); book(_hist_njets_thadpt_1, "d17-x02-y01"); book(_hist_njets_thadpt_2, "d18-x02-y01"); book(_hist_njets_thadpt_3, "d19-x02-y01"); book(_hist_njets_thadpt_4, "d20-x02-y01"); book(_hist_njets_ttpt_1, "d22-x02-y01"); book(_hist_njets_ttpt_2, "d23-x02-y01"); book(_hist_njets_ttpt_3, "d24-x02-y01"); book(_hist_njets_ttpt_4, "d25-x02-y01"); book(_hist_thady_thadpt_1, "d27-x02-y01"); book(_hist_thady_thadpt_2, "d28-x02-y01"); book(_hist_thady_thadpt_3, "d29-x02-y01"); book(_hist_thady_thadpt_4, "d30-x02-y01"); book(_hist_ttm_tty_1, "d32-x02-y01"); book(_hist_ttm_tty_2, "d33-x02-y01"); book(_hist_ttm_tty_3, "d34-x02-y01"); book(_hist_ttm_tty_4, "d35-x02-y01"); book(_hist_ttpt_ttm_1, "d37-x02-y01"); book(_hist_ttpt_ttm_2, "d38-x02-y01"); book(_hist_ttpt_ttm_3, "d39-x02-y01"); book(_hist_ttpt_ttm_4, "d40-x02-y01"); book(_histnorm_thadpt, "d42-x02-y01"); book(_histnorm_thady, "d44-x02-y01"); book(_histnorm_tleppt, "d46-x02-y01"); book(_histnorm_tlepy, "d48-x02-y01"); book(_histnorm_ttpt, "d50-x02-y01"); book(_histnorm_tty, "d54-x02-y01"); book(_histnorm_ttm, "d52-x02-y01"); book(_histnorm_njet, "d56-x02-y01"); book(_histnorm_njets_thadpt_1, "d58-x02-y01"); book(_histnorm_njets_thadpt_2, "d59-x02-y01"); book(_histnorm_njets_thadpt_3, "d60-x02-y01"); book(_histnorm_njets_thadpt_4, "d61-x02-y01"); book(_histnorm_njets_ttpt_1, "d63-x02-y01"); book(_histnorm_njets_ttpt_2, "d64-x02-y01"); book(_histnorm_njets_ttpt_3, "d65-x02-y01"); book(_histnorm_njets_ttpt_4, "d66-x02-y01"); book(_histnorm_thady_thadpt_1, "d68-x02-y01"); book(_histnorm_thady_thadpt_2, "d69-x02-y01"); book(_histnorm_thady_thadpt_3, "d70-x02-y01"); book(_histnorm_thady_thadpt_4, "d71-x02-y01"); book(_histnorm_ttm_tty_1, "d73-x02-y01"); book(_histnorm_ttm_tty_2, "d74-x02-y01"); book(_histnorm_ttm_tty_3, "d75-x02-y01"); book(_histnorm_ttm_tty_4, "d76-x02-y01"); book(_histnorm_ttpt_ttm_1, "d78-x02-y01"); book(_histnorm_ttpt_ttm_2, "d79-x02-y01"); book(_histnorm_ttpt_ttm_3, "d80-x02-y01"); book(_histnorm_ttpt_ttm_4, "d81-x02-y01"); } /// Perform the per-event analysis void analyze(const Event& event) { // leptons const SpecialDressedLeptons& dressedleptons_proj = applyProjection(event, "DressedLeptons"); std::vector dressedLeptons = dressedleptons_proj.dressedLeptons(); if(dressedLeptons.size() != 1) return; // neutrinos const Particles neutrinos = applyProjection(event, "Neutrinos").particlesByPt(); _nusum = FourMomentum(0., 0., 0., 0.); for(const Particle& neutrino : neutrinos) { _nusum += neutrino.momentum(); } _wl = _nusum + dressedLeptons[0].momentum(); // jets Cut jet_cut = (Cuts::abseta < 2.5) and (Cuts::pT > 25.*GeV); const Jets jets = applyProjection(event, "Jets").jetsByPt(jet_cut); Jets allJets; for (const Jet& jet : jets) { allJets.push_back(jet); } Jets bJets; for (const Jet& jet : allJets) { if (jet.bTagged()) bJets.push_back(jet); } if(bJets.size() < 2 || allJets.size() < 4) return; //construct top quark proxies double Kmin = numeric_limits::max(); for(const Jet& itaj : allJets) { for(const Jet& itbj : allJets) { if (itaj.momentum() == itbj.momentum()) continue; FourMomentum wh(itaj.momentum() + itbj.momentum()); for(const Jet& ithbj : bJets) { if(itaj.momentum() == ithbj.momentum() || itbj.momentum() == ithbj.momentum()) continue; FourMomentum th(wh + ithbj.momentum()); for(const Jet& itlbj : bJets) { if(itaj.momentum() == itlbj.momentum() || itbj.momentum() == itlbj.momentum() || ithbj.momentum() == itlbj.momentum()) continue; FourMomentum tl(_wl + itlbj.momentum()); double K = pow(wh.mass() - 80.4, 2) + pow(th.mass() - 172.5, 2) + pow(tl.mass() - 172.5, 2); if(K < Kmin) { Kmin = K; _tl = tl; _th = th; _wh = wh; } } } } } _hist_thadpt->fill(_th.pt()); _hist_thady->fill(abs(_th.rapidity()) ); _hist_tleppt->fill(_tl.pt() ); _hist_tlepy->fill(abs(_tl.rapidity()) ); _histnorm_thadpt->fill(_th.pt()); _histnorm_thady->fill(abs(_th.rapidity()) ); _histnorm_tleppt->fill(_tl.pt() ); _histnorm_tlepy->fill(abs(_tl.rapidity()) ); FourMomentum tt(_tl+_th); _hist_ttpt->fill(tt.pt() ); _hist_tty->fill(abs(tt.rapidity()) ); _hist_ttm->fill(tt.mass() ); _hist_njet->fill(min(allJets.size()-4., 4.)); _histnorm_ttpt->fill(tt.pt() ); _histnorm_tty->fill(abs(tt.rapidity()) ); _histnorm_ttm->fill(tt.mass() ); _histnorm_njet->fill(min(allJets.size()-4., 4.)); if(allJets.size() == 4) { _hist_njets_thadpt_1->fill(_th.pt()); _hist_njets_ttpt_1->fill(tt.pt()); _histnorm_njets_thadpt_1->fill(_th.pt()); _histnorm_njets_ttpt_1->fill(tt.pt()); } else if(allJets.size() == 5) { _hist_njets_thadpt_2->fill(_th.pt()); _hist_njets_ttpt_2->fill(tt.pt()); _histnorm_njets_thadpt_2->fill(_th.pt()); _histnorm_njets_ttpt_2->fill(tt.pt()); } else if(allJets.size() == 6) { _hist_njets_thadpt_3->fill(_th.pt()); _hist_njets_ttpt_3->fill(tt.pt()); _histnorm_njets_thadpt_3->fill(_th.pt()); _histnorm_njets_ttpt_3->fill(tt.pt()); } else //>= 4 jets { _hist_njets_thadpt_4->fill(_th.pt()); _hist_njets_ttpt_4->fill(tt.pt()); _histnorm_njets_thadpt_4->fill(_th.pt()); _histnorm_njets_ttpt_4->fill(tt.pt()); } if(abs(_th.rapidity()) < 0.5) { _hist_thady_thadpt_1->fill(_th.pt()); _histnorm_thady_thadpt_1->fill(_th.pt()); } else if(abs(_th.rapidity()) < 1.0) { _hist_thady_thadpt_2->fill(_th.pt()); _histnorm_thady_thadpt_2->fill(_th.pt()); } else if(abs(_th.rapidity()) < 1.5) { _hist_thady_thadpt_3->fill(_th.pt()); _histnorm_thady_thadpt_3->fill(_th.pt()); } else if(abs(_th.rapidity()) < 2.5) { _hist_thady_thadpt_4->fill(_th.pt()); _histnorm_thady_thadpt_4->fill(_th.pt()); } if(tt.mass() >= 300. && tt.mass() < 450.) { _hist_ttm_tty_1->fill(abs(tt.rapidity())); _histnorm_ttm_tty_1->fill(abs(tt.rapidity())); } else if(tt.mass() >= 450. && tt.mass() < 625.) { _hist_ttm_tty_2->fill(abs(tt.rapidity())); _histnorm_ttm_tty_2->fill(abs(tt.rapidity())); } else if(tt.mass() >= 625. && tt.mass() < 850.) { _hist_ttm_tty_3->fill(abs(tt.rapidity())); _histnorm_ttm_tty_3->fill(abs(tt.rapidity())); } else if(tt.mass() >= 850. && tt.mass() < 2000.) { _hist_ttm_tty_4->fill(abs(tt.rapidity())); _histnorm_ttm_tty_4->fill(abs(tt.rapidity())); } if(tt.pt() < 35.) { _hist_ttpt_ttm_1->fill(tt.mass()); _histnorm_ttpt_ttm_1->fill(tt.mass()); } else if(tt.pt() < 80.) { _hist_ttpt_ttm_2->fill(tt.mass()); _histnorm_ttpt_ttm_2->fill(tt.mass()); } else if(tt.pt() < 140.) { _hist_ttpt_ttm_3->fill(tt.mass()); _histnorm_ttpt_ttm_3->fill(tt.mass()); } else if(tt.pt() < 500.) { _hist_ttpt_ttm_4->fill(tt.mass()); _histnorm_ttpt_ttm_4->fill(tt.mass()); } } /// Normalise histograms etc., after the run void finalize() { scale(_hist_thadpt, crossSection()/sumOfWeights()); scale(_hist_thady, crossSection()/sumOfWeights()); scale(_hist_tleppt, crossSection()/sumOfWeights()); scale(_hist_tlepy, crossSection()/sumOfWeights()); scale(_hist_ttpt, crossSection()/sumOfWeights()); scale(_hist_tty, crossSection()/sumOfWeights()); scale(_hist_ttm, crossSection()/sumOfWeights()); scale(_hist_njet, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_1, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_2, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_3, crossSection()/sumOfWeights()); scale(_hist_njets_thadpt_4, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_1, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_2, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_3, crossSection()/sumOfWeights()); scale(_hist_njets_ttpt_4, crossSection()/sumOfWeights()); scale(_hist_thady_thadpt_1, crossSection()/sumOfWeights()/0.5); scale(_hist_thady_thadpt_2, crossSection()/sumOfWeights()/0.5); scale(_hist_thady_thadpt_3, crossSection()/sumOfWeights()/0.5); scale(_hist_thady_thadpt_4, crossSection()/sumOfWeights()/1.0); scale(_hist_ttm_tty_1, crossSection()/sumOfWeights()/150.); scale(_hist_ttm_tty_2, crossSection()/sumOfWeights()/175.); scale(_hist_ttm_tty_3, crossSection()/sumOfWeights()/225.); scale(_hist_ttm_tty_4, crossSection()/sumOfWeights()/1150.); scale(_hist_ttpt_ttm_1, crossSection()/sumOfWeights()/35.); scale(_hist_ttpt_ttm_2, crossSection()/sumOfWeights()/45.); scale(_hist_ttpt_ttm_3, crossSection()/sumOfWeights()/60.); scale(_hist_ttpt_ttm_4, crossSection()/sumOfWeights()/360.); scale(_histnorm_thadpt, 1./_histnorm_thadpt->sumW(false)); scale(_histnorm_thady, 1./_histnorm_thady->sumW(false)); scale(_histnorm_tleppt, 1./_histnorm_tleppt->sumW(false)); scale(_histnorm_tlepy, 1./_histnorm_tlepy->sumW(false)); scale(_histnorm_ttpt, 1./_histnorm_ttpt->sumW(false)); scale(_histnorm_tty, 1./_histnorm_tty->sumW(false)); scale(_histnorm_ttm, 1./_histnorm_ttm->sumW(false)); scale(_histnorm_njet, 1./_histnorm_njet->sumW(false)); double sum_njets_thadpt = _histnorm_njets_thadpt_1->sumW(false) + _histnorm_njets_thadpt_2->sumW(false) + _histnorm_njets_thadpt_3->sumW(false) + _histnorm_njets_thadpt_4->sumW(false); scale(_histnorm_njets_thadpt_1, 1./sum_njets_thadpt); scale(_histnorm_njets_thadpt_2, 1./sum_njets_thadpt); scale(_histnorm_njets_thadpt_3, 1./sum_njets_thadpt); scale(_histnorm_njets_thadpt_4, 1./sum_njets_thadpt); double sum_njets_ttpt = _histnorm_njets_ttpt_1->sumW(false) + _histnorm_njets_ttpt_2->sumW(false) + _histnorm_njets_ttpt_3->sumW(false) + _histnorm_njets_ttpt_4->sumW(false); scale(_histnorm_njets_ttpt_1, 1./sum_njets_ttpt); scale(_histnorm_njets_ttpt_2, 1./sum_njets_ttpt); scale(_histnorm_njets_ttpt_3, 1./sum_njets_ttpt); scale(_histnorm_njets_ttpt_4, 1./sum_njets_ttpt); double sum_thady_thadpt = _histnorm_thady_thadpt_1->sumW(false) + _histnorm_thady_thadpt_2->sumW(false) + _histnorm_thady_thadpt_3->sumW(false) + _histnorm_thady_thadpt_4->sumW(false); scale(_histnorm_thady_thadpt_1, 1./sum_thady_thadpt/0.5); scale(_histnorm_thady_thadpt_2, 1./sum_thady_thadpt/0.5); scale(_histnorm_thady_thadpt_3, 1./sum_thady_thadpt/0.5); scale(_histnorm_thady_thadpt_4, 1./sum_thady_thadpt/1.0); double sum_ttm_tty = _histnorm_ttm_tty_1->sumW(false) + _histnorm_ttm_tty_2->sumW(false) + _histnorm_ttm_tty_3->sumW(false) + _histnorm_ttm_tty_4->sumW(false); scale(_histnorm_ttm_tty_1, 1./sum_ttm_tty/150.); scale(_histnorm_ttm_tty_2, 1./sum_ttm_tty/175.); scale(_histnorm_ttm_tty_3, 1./sum_ttm_tty/225.); scale(_histnorm_ttm_tty_4, 1./sum_ttm_tty/1150.); double sum_ttpt_ttm = _histnorm_ttpt_ttm_1->sumW(false) + _histnorm_ttpt_ttm_2->sumW(false) + _histnorm_ttpt_ttm_3->sumW(false) + _histnorm_ttpt_ttm_4->sumW(false); scale(_histnorm_ttpt_ttm_1, 1./sum_ttpt_ttm/35.); scale(_histnorm_ttpt_ttm_2, 1./sum_ttpt_ttm/45.); scale(_histnorm_ttpt_ttm_3, 1./sum_ttpt_ttm/60.); scale(_histnorm_ttpt_ttm_4, 1./sum_ttpt_ttm/360.); } private: FourMomentum _tl; FourMomentum _th; FourMomentum _wl; FourMomentum _wh; FourMomentum _nusum; Histo1DPtr _hist_thadpt; Histo1DPtr _hist_thady; Histo1DPtr _hist_tleppt; Histo1DPtr _hist_tlepy; Histo1DPtr _hist_ttpt; Histo1DPtr _hist_tty; Histo1DPtr _hist_ttm; Histo1DPtr _hist_njet; Histo1DPtr _hist_njets_thadpt_1; Histo1DPtr _hist_njets_thadpt_2; Histo1DPtr _hist_njets_thadpt_3; Histo1DPtr _hist_njets_thadpt_4; Histo1DPtr _hist_njets_ttpt_1; Histo1DPtr _hist_njets_ttpt_2; Histo1DPtr _hist_njets_ttpt_3; Histo1DPtr _hist_njets_ttpt_4; Histo1DPtr _hist_thady_thadpt_1; Histo1DPtr _hist_thady_thadpt_2; Histo1DPtr _hist_thady_thadpt_3; Histo1DPtr _hist_thady_thadpt_4; Histo1DPtr _hist_ttm_tty_1; Histo1DPtr _hist_ttm_tty_2; Histo1DPtr _hist_ttm_tty_3; Histo1DPtr _hist_ttm_tty_4; Histo1DPtr _hist_ttpt_ttm_1; Histo1DPtr _hist_ttpt_ttm_2; Histo1DPtr _hist_ttpt_ttm_3; Histo1DPtr _hist_ttpt_ttm_4; Histo1DPtr _histnorm_thadpt; Histo1DPtr _histnorm_thady; Histo1DPtr _histnorm_tleppt; Histo1DPtr _histnorm_tlepy; Histo1DPtr _histnorm_ttpt; Histo1DPtr _histnorm_tty; Histo1DPtr _histnorm_ttm; Histo1DPtr _histnorm_njet; Histo1DPtr _histnorm_njets_thadpt_1; Histo1DPtr _histnorm_njets_thadpt_2; Histo1DPtr _histnorm_njets_thadpt_3; Histo1DPtr _histnorm_njets_thadpt_4; Histo1DPtr _histnorm_njets_ttpt_1; Histo1DPtr _histnorm_njets_ttpt_2; Histo1DPtr _histnorm_njets_ttpt_3; Histo1DPtr _histnorm_njets_ttpt_4; Histo1DPtr _histnorm_thady_thadpt_1; Histo1DPtr _histnorm_thady_thadpt_2; Histo1DPtr _histnorm_thady_thadpt_3; Histo1DPtr _histnorm_thady_thadpt_4; Histo1DPtr _histnorm_ttm_tty_1; Histo1DPtr _histnorm_ttm_tty_2; Histo1DPtr _histnorm_ttm_tty_3; Histo1DPtr _histnorm_ttm_tty_4; Histo1DPtr _histnorm_ttpt_ttm_1; Histo1DPtr _histnorm_ttpt_ttm_2; Histo1DPtr _histnorm_ttpt_ttm_3; Histo1DPtr _histnorm_ttpt_ttm_4; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_I1491950); } diff --git a/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc b/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc --- a/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc +++ b/analyses/pluginCMS/CMS_2016_PAS_TOP_15_006.cc @@ -1,183 +1,183 @@ #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Tools/ParticleName.hh" #include "Rivet/Tools/ParticleIdUtils.hh" namespace Rivet { namespace { //< only visible in this compilation unit /// @brief Special dressed lepton finder /// /// Find dressed leptons by clustering all leptons and photons class SpecialDressedLeptons : public FinalState { public: /// Constructor SpecialDressedLeptons(const FinalState& fs, const Cut& cut) : FinalState(cut) { setName("SpecialDressedLeptons"); IdentifiedFinalState ifs(fs); ifs.acceptIdPair(PID::PHOTON); ifs.acceptIdPair(PID::ELECTRON); ifs.acceptIdPair(PID::MUON); declare(ifs, "IFS"); declare(FastJets(ifs, FastJets::ANTIKT, 0.1), "LeptonJets"); } /// Clone on the heap virtual unique_ptr clone() const { return unique_ptr(new SpecialDressedLeptons(*this)); } /// Retrieve the dressed leptons const vector& dressedLeptons() const { return _clusteredLeptons; } /// Perform the calculation void project(const Event& e) { _theParticles.clear(); _clusteredLeptons.clear(); vector allClusteredLeptons; const Jets jets = applyProjection(e, "LeptonJets").jetsByPt(5*GeV); for (const Jet& jet : jets) { Particle lepCand; for (const Particle& cand : jet.particles()) { const int absPdgId = cand.abspid(); if (absPdgId == PID::ELECTRON || absPdgId == PID::MUON) { if (cand.pt() > lepCand.pt()) lepCand = cand; } } // Central lepton must be the major component - if ((lepCand.pt() < jet.pt()/2.) || (lepCand.pdgId() == 0)) continue; + if ((lepCand.pt() < jet.pt()/2.) || (lepCand.pid() == 0)) continue; DressedLepton lepton = DressedLepton(lepCand); for (const Particle& cand : jet.particles()) { if (cand == lepCand) continue; lepton.addPhoton(cand, true); } allClusteredLeptons.push_back(lepton); } for (const DressedLepton& lepton : allClusteredLeptons) { if (accept(lepton)) { _clusteredLeptons.push_back(lepton); _theParticles.push_back(lepton.constituentLepton()); _theParticles += lepton.constituentPhotons(); } } } private: /// Container which stores the clustered lepton objects vector _clusteredLeptons; }; } /// Jet multiplicity in lepton+jets ttbar at 8 TeV class CMS_2016_PAS_TOP_15_006 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2016_PAS_TOP_15_006); /// @name Analysis methods //@{ /// Set up projections and book histograms void init() { // Complete final state FinalState fs; Cut superLooseLeptonCuts = Cuts::pt > 5*GeV; SpecialDressedLeptons dressedleptons(fs, superLooseLeptonCuts); declare(dressedleptons, "DressedLeptons"); // Projection for jets VetoedFinalState fsForJets(fs); fsForJets.addVetoOnThisFinalState(dressedleptons); declare(FastJets(fsForJets, FastJets::ANTIKT, 0.5), "Jets"); // Booking of histograms book(_normedElectronMuonHisto, "normedElectronMuonHisto", 7, 3.5, 10.5, "Normalized differential cross section in lepton+jets channel", "Jet multiplicity", "Normed units"); book(_absXSElectronMuonHisto , "absXSElectronMuonHisto", 7, 3.5, 10.5, "Differential cross section in lepton+jets channel", "Jet multiplicity", "pb"); } /// Per-event analysis void analyze(const Event& event) { // Select ttbar -> lepton+jets const SpecialDressedLeptons& dressedleptons = applyProjection(event, "DressedLeptons"); vector selleptons; for (const DressedLepton& dressedlepton : dressedleptons.dressedLeptons()) { // Select good leptons if (dressedlepton.pT() > 30*GeV && dressedlepton.abseta() < 2.4) selleptons += dressedlepton.mom(); // Veto loose leptons else if (dressedlepton.pT() > 15*GeV && dressedlepton.abseta() < 2.5) vetoEvent; } if (selleptons.size() != 1) vetoEvent; // Identify hardest tight lepton const FourMomentum lepton = selleptons[0]; // Jets const FastJets& jets = applyProjection(event, "Jets"); const Jets jets30 = jets.jetsByPt(30*GeV); int nJets = 0, nBJets = 0; for (const Jet& jet : jets30) { if (jet.abseta() > 2.5) continue; if (deltaR(jet.momentum(), lepton) < 0.5) continue; nJets += 1; if (jet.bTagged(Cuts::pT > 5*GeV)) nBJets += 1; } // Require >= 4 resolved jets, of which two must be b-tagged if (nJets < 4 || nBJets < 2) vetoEvent; // Fill histograms _normedElectronMuonHisto->fill(min(nJets, 10)); _absXSElectronMuonHisto ->fill(min(nJets, 10)); } 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"); const double xsPerWeight = ttbarXS/picobarn / sumOfWeights(); scale(_absXSElectronMuonHisto, xsPerWeight); normalize(_normedElectronMuonHisto); } //@} private: /// Histograms Histo1DPtr _normedElectronMuonHisto, _absXSElectronMuonHisto; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2016_PAS_TOP_15_006); } diff --git a/analyses/pluginCMS/CMS_2017_I1605749.cc b/analyses/pluginCMS/CMS_2017_I1605749.cc --- a/analyses/pluginCMS/CMS_2017_I1605749.cc +++ b/analyses/pluginCMS/CMS_2017_I1605749.cc @@ -1,145 +1,145 @@ // -*- C++ -*- // Rivet framework #include "Rivet/Analysis.hh" // Projections #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { using namespace Cuts; class CMS_2017_I1605749 : public Analysis { public: // Constructor CMS_2017_I1605749() : Analysis("CMS_2017_I1605749") { } // Book histograms and initialise projections before the run void init() { // Projections const FinalState fs(-5.0, 5.0, 0.0*GeV); declare(FastJets(fs, FastJets::ANTIKT, 0.5), "Jets"); // Jet Charge Histos for (int i = 1; i <= 18; i++) { book(_h_Charge[i - 1], i, 1, 1); } } // Perform the per-event analysis void analyze(const Event& event) { const Jets& jets = applyProjection(event, "Jets").jetsByPt(10.0*GeV); if (jets.size() < 2) vetoEvent; double leadingpt = jets[0].pt()/GeV; double subleadingpt = jets[1].pt()/GeV; if (jets.size() < 2 || jets[0].abseta() >= 1.5 || jets[1].abseta() >= 1.5 || leadingpt < 400.0 || subleadingpt < 100.0) { vetoEvent; } vector constituents1 = jets[0].constituents(); std::vector numerator(9, 0), denominator(9, 0); double t_jetcharge1, t_jetcharge1k6, t_jetcharge1k3; double t_jetchargeL1, t_jetchargeL1k6, t_jetchargeL1k3; double t_jetchargeT1, t_jetchargeT1k6, t_jetchargeT1k3; denominator[0] = leadingpt; denominator[1] = std::pow(leadingpt, 0.6); denominator[2] = std::pow(leadingpt, 0.3); if (constituents1.size() > 0) { for (unsigned j = 0; j < constituents1.size(); j++) { - if (std::abs(constituents1[j].pdgId()) > 9 && - std::abs(constituents1[j].pdgId())!= 21) { + if (std::abs(constituents1[j].pid()) > 9 && + std::abs(constituents1[j].pid())!= 21) { if (constituents1[j].pt() > 1*GeV) { double charge = constituents1[j].charge(); double mom = constituents1[j].pt(); double dotproduct = constituents1[j].p3().dot(jets[0].p3()) / jets[0].p(); double crossproduct = constituents1[j].p3().cross(jets[0].p3()).mod() / jets[0].p(); numerator[0] += (mom * charge); numerator[1] += ((std::pow(mom, 0.6)) * charge); numerator[2] += ((std::pow(mom, 0.3)) * charge); numerator[3] += (dotproduct * charge); numerator[4] += ((std::pow(dotproduct, 0.6)) * charge); numerator[5] += ((std::pow(dotproduct, 0.3)) * charge); denominator[3] += dotproduct; denominator[4] += (std::pow(dotproduct, 0.6)); denominator[5] += (std::pow(dotproduct, 0.3)); numerator[6] += (crossproduct * charge); numerator[7] += ((std::pow(crossproduct, 0.6)) * charge); numerator[8] += ((std::pow(crossproduct, 0.3)) * charge); denominator[6] += crossproduct; denominator[7] += (std::pow(crossproduct, 0.6)); denominator[8] += (std::pow(crossproduct, 0.3)); } } } } t_jetcharge1 = (denominator[0] > 0) ? numerator[0] / denominator[0] : 0; t_jetcharge1k6 = (denominator[1] > 0) ? numerator[1] / denominator[1] : 0; t_jetcharge1k3 = (denominator[2] > 0) ? numerator[2] / denominator[2] : 0; t_jetchargeL1 = (denominator[3] > 0) ? numerator[3] / denominator[3] : 0; t_jetchargeL1k6 = (denominator[4] > 0) ? numerator[4] / denominator[4] : 0; t_jetchargeL1k3 = (denominator[5] > 0) ? numerator[5] / denominator[5] : 0; t_jetchargeT1 = (denominator[6] > 0) ? numerator[6] / denominator[6] : 0; t_jetchargeT1k6 = (denominator[7] > 0) ? numerator[7] / denominator[7] : 0; t_jetchargeT1k3 = (denominator[8] > 0) ? numerator[8] / denominator[8] : 0; _h_Charge[0]->fill(t_jetcharge1); _h_Charge[1]->fill(t_jetcharge1k6); _h_Charge[2]->fill(t_jetcharge1k3); _h_Charge[3]->fill(t_jetchargeL1); _h_Charge[4]->fill(t_jetchargeL1k6); _h_Charge[5]->fill(t_jetchargeL1k3); _h_Charge[6]->fill(t_jetchargeT1); _h_Charge[7]->fill(t_jetchargeT1k6); _h_Charge[8]->fill(t_jetchargeT1k3); if (leadingpt > 400 && leadingpt < 700) { _h_Charge[9]->fill(t_jetcharge1k6); _h_Charge[12]->fill(t_jetchargeL1k6); _h_Charge[15]->fill(t_jetchargeT1k6); } else if (leadingpt > 700 && leadingpt < 1000) { _h_Charge[10]->fill(t_jetcharge1k6); _h_Charge[13]->fill(t_jetchargeL1k6); _h_Charge[16]->fill(t_jetchargeT1k6); } else if (leadingpt > 1000 && leadingpt < 1800) { _h_Charge[11]->fill(t_jetcharge1k6); _h_Charge[14]->fill(t_jetchargeL1k6); _h_Charge[17]->fill(t_jetchargeT1k6); } } // Normalise histograms etc., after the run void finalize() { for (int j = 0; j < 18; j++) { normalize(_h_Charge[j]); for (size_t i = 0; i < _h_Charge[j]-> numBins(); i++) { _h_Charge[j]->bin(i).scaleW(1.0 / _h_Charge[j]->bin(i).width()); } } } private: Histo1DPtr _h_Charge[18]; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2017_I1605749); } diff --git a/analyses/pluginCMS/CMS_2018_I1663958.cc b/analyses/pluginCMS/CMS_2018_I1663958.cc --- a/analyses/pluginCMS/CMS_2018_I1663958.cc +++ b/analyses/pluginCMS/CMS_2018_I1663958.cc @@ -1,452 +1,452 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" namespace Rivet { class CMS_2018_I1663958 : public Analysis { public: struct Histo1DGroup { Histo1DGroup() {} Histo1DGroup(CMS_2018_I1663958* an, const vector& hnames, const vector& xbinranges); void fill(double x, double y, double w = 1.); void scale(double s, bool diff = true); void norm(bool diff = true); void gapfractionfromjetpt(Scatter2DPtr hgap, int njet); CMS_2018_I1663958* m_an; vector m_xbins; vector m_histos; }; CMS_2018_I1663958() : Analysis("CMS_2018_I1663958"), m_thad_decay(3), m_tlep_decay(3), m_tt_jets(4) {} virtual void init() { const FinalState fs(Cuts::abseta < 6.); declare(fs, "FS"); const VisibleFinalState vfs(Cuts::abseta < 6.); declare(vfs, "vFS"); VetoedFinalState invisibles(fs); invisibles.addVetoOnThisFinalState(vfs); declare(invisibles, "Invisibles"); IdentifiedFinalState all_photons(vfs); all_photons.acceptId(PID::PHOTON); IdentifiedFinalState leptons(vfs); leptons.acceptIds({PID::ELECTRON, -PID::ELECTRON, PID::MUON,-PID::MUON}); DressedLeptons dressed_leptons(all_photons, leptons, m_lepdressdr, Cuts::abseta < m_lepetamax && Cuts::pT > m_vetolepptmin*GeV, true); declare(dressed_leptons, "MyLeptons"); VetoedFinalState photons(all_photons); photons.addVetoOnThisFinalState(dressed_leptons); declare(photons, "MyPhotons"); VetoedFinalState isolationparticles(vfs); isolationparticles.addVetoOnThisFinalState(dressed_leptons); declare(isolationparticles, "IsoParticles"); declare(FastJets(vfs, FastJets::ANTIKT, m_jetdr), "Jets"); book(m_hist_thadpt , "d01-x01-y01"); book(m_hist_thady , "d03-x01-y01"); book(m_hist_tleppt , "d05-x01-y01"); book(m_hist_tlepy , "d07-x01-y01"); book(m_hist_ttm , "d09-x01-y01"); book(m_hist_ttpt , "d11-x01-y01"); book(m_hist_tty , "d13-x01-y01"); book(m_hist_njet , "d15-x01-y01"); /// @todo Memory leak m_hist_njet_ttm = Histo1DGroup(this, {"d17-x01-y01", "d18-x01-y01", "d19-x01-y01", "d20-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5}); m_hist_njet_thadpt = Histo1DGroup(this, {"d22-x01-y01", "d23-x01-y01", "d24-x01-y01", "d25-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5}); m_hist_njet_ttpt = Histo1DGroup(this, {"d27-x01-y01", "d28-x01-y01", "d29-x01-y01", "d30-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5}); m_hist_thady_thadpt = Histo1DGroup(this, {"d32-x01-y01", "d33-x01-y01", "d34-x01-y01", "d35-x01-y01"}, {0.0,0.5, 1.0, 1.5, 2.5}); m_hist_ttm_tty = Histo1DGroup(this, {"d37-x01-y01", "d38-x01-y01", "d39-x01-y01", "d40-x01-y01"}, {300., 450., 625., 850., 2000.}); m_hist_thadpt_ttm = Histo1DGroup(this, {"d42-x01-y01", "d43-x01-y01", "d44-x01-y01", "d45-x01-y01"}, {0., 90., 180., 270., 800.}); m_hist_jetspt = Histo1DGroup(this, {"d47-x01-y01", "d48-x01-y01", "d49-x01-y01", "d50-x01-y01", "d51-x01-y01", "d52-x01-y01", "d53-x01-y01", "d54-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_hist_jetseta = Histo1DGroup(this, {"d56-x01-y01", "d57-x01-y01", "d58-x01-y01", "d59-x01-y01", "d60-x01-y01", "d61-x01-y01", "d62-x01-y01", "d63-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_hist_jetsdrjets = Histo1DGroup(this, {"d65-x01-y01", "d66-x01-y01", "d67-x01-y01", "d68-x01-y01", "d69-x01-y01", "d70-x01-y01", "d71-x01-y01", "d72-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_hist_jetsdrtops = Histo1DGroup(this, {"d74-x01-y01", "d75-x01-y01", "d76-x01-y01", "d77-x01-y01", "d78-x01-y01", "d79-x01-y01", "d80-x01-y01", "d81-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_hist_njetspt = Histo1DGroup(this, {"d169-x01-y01", "d170-x01-y01", "d171-x01-y01", "d172-x01-y01"}, {0., 40., 60., 80., 120.}); book(m_nhist_thadpt , "d83-x01-y01"); book(m_nhist_thady , "d85-x01-y01"); book(m_nhist_tleppt , "d87-x01-y01"); book(m_nhist_tlepy , "d89-x01-y01"); book(m_nhist_ttm , "d91-x01-y01"); book(m_nhist_ttpt , "d93-x01-y01"); book(m_nhist_tty , "d95-x01-y01"); book(m_nhist_njet , "d97-x01-y01"); /// @todo Memory leak m_nhist_njet_ttm = Histo1DGroup(this, {"d99-x01-y01", "d100-x01-y01", "d101-x01-y01", "d102-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5}); m_nhist_njet_thadpt = Histo1DGroup(this, {"d104-x01-y01", "d105-x01-y01", "d106-x01-y01", "d107-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5}); m_nhist_njet_ttpt = Histo1DGroup(this, {"d109-x01-y01", "d110-x01-y01", "d111-x01-y01", "d112-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5}); m_nhist_thady_thadpt = Histo1DGroup(this, {"d114-x01-y01", "d115-x01-y01", "d116-x01-y01", "d117-x01-y01"}, {0.0,0.5, 1.0, 1.5, 2.5}); m_nhist_ttm_tty = Histo1DGroup(this, {"d119-x01-y01", "d120-x01-y01", "d121-x01-y01", "d122-x01-y01"}, {300., 450., 625., 850., 2000.}); m_nhist_thadpt_ttm = Histo1DGroup(this, {"d124-x01-y01", "d125-x01-y01", "d126-x01-y01", "d127-x01-y01"}, {0., 90., 180., 270., 800.}); m_nhist_jetspt = Histo1DGroup(this, {"d129-x01-y01", "d130-x01-y01", "d131-x01-y01", "d132-x01-y01", "d133-x01-y01", "d134-x01-y01", "d135-x01-y01", "d136-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_nhist_jetseta = Histo1DGroup(this, {"d138-x01-y01", "d139-x01-y01", "d140-x01-y01", "d141-x01-y01", "d142-x01-y01", "d143-x01-y01", "d144-x01-y01", "d145-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_nhist_jetsdrjets = Histo1DGroup(this, {"d147-x01-y01", "d148-x01-y01", "d149-x01-y01", "d150-x01-y01", "d151-x01-y01", "d152-x01-y01", "d153-x01-y01", "d154-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); m_nhist_jetsdrtops = Histo1DGroup(this, {"d156-x01-y01", "d157-x01-y01", "d158-x01-y01", "d159-x01-y01", "d160-x01-y01", "d161-x01-y01", "d162-x01-y01", "d163-x01-y01"}, {-0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5}); book(m_hist_gap1, "d165-x01-y01"); book(m_hist_gap2, "d167-x01-y01"); } void analyze(const Event& event) { m_leptons.clear(); m_vetoleptons.clear(); m_neutrinos.clear(); m_bjets.clear(); m_ljets.clear(); m_additionalobjects.clear(); m_additionaljets.clear(); const Particles& isopars = applyProjection(event, "IsoParticles").particles(); const vector& dressedleptons = applyProjection(event, "MyLeptons").dressedLeptons(); for (const DressedLepton& lep : dressedleptons) { double isolation = accumulate(isopars.begin(), isopars.end(), 0., [&](double iso, const Particle& par) {if(deltaR(lep, par) < m_lepisodr){return iso + par.pt();} else {return iso;}}); isolation = isolation/lep.pt(); if (isolation > m_lepisomax) continue; if (lep.pt() > m_lepptmin && lep.abseta() < m_lepetamax) m_leptons.push_back(lep); else if(lep.pt() > m_vetolepptmin && lep.abseta() < m_vetolepetamax) m_vetoleptons.push_back(lep); } const Particles& photons = applyProjection(event, "MyPhotons").particles(Cuts::abseta < m_phetamax && Cuts::pT > m_phptmin*GeV); for (const Particle& ph : photons) { double isolation = accumulate(isopars.begin(), isopars.end(), 0., [&](double iso, const Particle& par) {if(deltaR(ph, par) < m_phisodr){return iso + par.pt();} else {return iso;}}); isolation = isolation/ph.pt() - 1.; if (isolation > m_phisomax) continue; m_additionalobjects.push_back(ph); } const Particles& invfspars = apply(event, "Invisibles").particles(); for (const Particle& par : invfspars) m_neutrinos.push_back(par); const Jets& allJets = apply(event, "Jets").jetsByPt(Cuts::abseta < m_jetetamax && Cuts::pT > m_jetptmin*GeV); for (const Jet& jet : allJets) { //clean jets from leptons if (find_if(m_leptons.begin(), m_leptons.end(), [&](const Particle& par){return deltaR(jet, par) < m_jetdr;}) != m_leptons.end()) continue; //clean jets from veto leptons if (find_if(m_vetoleptons.begin(), m_vetoleptons.end(), [&](const Particle& par){return deltaR(jet, par) < m_jetdr;}) != m_vetoleptons.end()) continue; //clean other objects (photons) if (find_if(m_additionalobjects.begin(), m_additionalobjects.end(), [&](const Particle& par){return deltaR(jet, par) < m_jetdr;}) != m_additionalobjects.end()) continue; if (jet.bTagged()) { m_bjets.push_back(jet); } else { m_ljets.push_back(jet); } } //Semi-leptonic reconstruction if (m_leptons.size() != 1 || m_vetoleptons.size() != 0 || m_bjets.size() < 2 || m_ljets.size() < 2) return; FourMomentum nusum = accumulate(m_neutrinos.begin(), m_neutrinos.end(), FourMomentum(0.,0.,0.,0.), [&](FourMomentum& invmom, const Particle& par) {return invmom += par.momentum();}); FourMomentum wl = nusum + m_leptons[0].momentum(); double Kmin = numeric_limits::max(); for (size_t a = 0 ; a < m_ljets.size() ; ++a){ const Jet& lja = m_ljets[a]; for (size_t b = 0 ; b < a ; ++b) { const Jet& ljb = m_ljets[b]; FourMomentum wh(lja.momentum() + ljb.momentum()); for (const Jet& bjh : m_bjets) { FourMomentum th(wh + bjh.momentum()); for (const Jet& bjl : m_bjets) { if (&bjh == &bjl) continue; FourMomentum tl(wl + bjl.momentum()); double K = pow(wh.mass() - 80.4, 2) + pow(th.mass() - 172.5, 2) + pow(tl.mass() - 172.5, 2); if (K < Kmin) { Kmin = K; m_thad = Particle(6, th); m_thad_decay[0] = Particle(5, bjh); m_thad_decay[1] = lja.pt() > ljb.pt() ? Particle(1, lja) : Particle(1, ljb); m_thad_decay[2] = lja.pt() <= ljb.pt() ? Particle(1, lja) : Particle(1, ljb); m_tlep = Particle(-6, tl); m_tlep_decay[0] = Particle(5, bjl); m_tlep_decay[1] = m_leptons[0]; - m_tlep_decay[2] = Particle(-1*(m_leptons[0].pdgId()+1), nusum); + m_tlep_decay[2] = Particle(-1*(m_leptons[0].pid()+1), nusum); } } } } } m_tt_jets[0] = m_tlep_decay[0]; m_tt_jets[1] = m_thad_decay[0]; m_tt_jets[2] = m_thad_decay[1]; m_tt_jets[3] = m_thad_decay[2]; const double eps = 1E-5; for (const Jet& jet : m_bjets) { if(jet.pt() < m_addjetptmin || jet.abseta() > m_addjetetamax) continue; if(find_if(m_tt_jets.begin(), m_tt_jets.end(), [&](const Particle& par){return deltaR(jet, par) < eps;}) != m_tt_jets.end()) continue; m_additionaljets.push_back(Particle(5, jet.momentum())); } for (const Jet& jet : m_ljets) { if(jet.pt() < m_addjetptmin || jet.abseta() > m_addjetetamax) continue; if(find_if(m_tt_jets.begin(), m_tt_jets.end(), [&](const Particle& par){return deltaR(jet, par) < eps;}) != m_tt_jets.end()) continue; if(jet.cTagged()) { m_additionaljets.push_back(Particle(4, jet.momentum())); } else { m_additionaljets.push_back(Particle(1, jet.momentum())); } } sort(m_additionaljets.begin(), m_additionaljets.end(), [](const Particle& ja, const Particle& jb) {return ja.pt() > jb.pt();}); FourMomentum tt(m_thad.momentum() + m_tlep.momentum()); m_hist_thadpt->fill(m_thad.pt()); m_nhist_thadpt->fill(m_thad.pt()); m_hist_thady->fill(abs(m_thad.rapidity())); m_nhist_thady->fill(abs(m_thad.rapidity())); m_hist_tleppt->fill(m_tlep.pt()); m_nhist_tleppt->fill(m_tlep.pt()); m_hist_tlepy->fill(abs(m_tlep.rapidity())); m_nhist_tlepy->fill(abs(m_tlep.rapidity())); m_hist_ttm->fill(tt.mass()); m_nhist_ttm->fill(tt.mass()); m_hist_ttpt->fill(tt.pt()); m_nhist_ttpt->fill(tt.pt()); m_hist_tty->fill(abs(tt.rapidity())); m_nhist_tty->fill(abs(tt.rapidity())); m_hist_njet->fill(min(m_additionaljets.size(), (size_t)5)); m_nhist_njet->fill(min(m_additionaljets.size(), (size_t)5)); int njet = min((size_t)3, m_additionaljets.size()); m_hist_njet_ttm.fill(njet, tt.mass()); m_nhist_njet_ttm.fill(njet, tt.mass()); m_hist_njet_thadpt.fill(njet, m_thad.pt()); m_nhist_njet_thadpt.fill(njet, m_thad.pt()); m_hist_njet_ttpt.fill(njet, tt.pt()); m_nhist_njet_ttpt.fill(njet, tt.pt()); m_hist_thady_thadpt.fill(abs(m_thad.rapidity()), m_thad.pt()); m_nhist_thady_thadpt.fill(abs(m_thad.rapidity()), m_thad.pt()); m_hist_ttm_tty.fill(tt.mass(), abs(tt.rapidity())); m_nhist_ttm_tty.fill(tt.mass(), abs(tt.rapidity())); m_hist_thadpt_ttm.fill(m_thad.pt(), tt.mass()); m_nhist_thadpt_ttm.fill(m_thad.pt(), tt.mass()); int jpos = -1; for (const Particles& jets : {m_tt_jets, m_additionaljets}) { for (const Particle& jet : jets) { jpos++; m_hist_jetspt.fill(jpos, jet.pt()); m_nhist_jetspt.fill(jpos, jet.pt()); m_hist_jetseta.fill(jpos, abs(jet.eta())); m_nhist_jetseta.fill(jpos, abs(jet.eta())); double drmin = 1E10; for (const Particle& tjet : m_tt_jets) { double dr = deltaR(jet, tjet); if(dr > eps && dr < drmin) { drmin = dr; } } m_hist_jetsdrjets.fill(jpos, drmin); m_nhist_jetsdrjets.fill(jpos, drmin); m_hist_jetsdrtops.fill(jpos, min(deltaR(jet, m_thad), deltaR(jet, m_tlep))); m_nhist_jetsdrtops.fill(jpos, min(deltaR(jet, m_thad), deltaR(jet, m_tlep))); } } for (double ptcut : {30, 50, 75, 100}) { m_hist_njetspt.fill(ptcut , count_if(m_additionaljets.begin(), m_additionaljets.end(), [&ptcut](const Particle& j) {return j.pt() > ptcut;})); } } virtual void finalize() { m_hist_jetspt.gapfractionfromjetpt(m_hist_gap1, 1); m_hist_jetspt.gapfractionfromjetpt(m_hist_gap2, 2); scale(m_hist_thadpt, crossSection()/sumOfWeights()); scale(m_hist_thady, crossSection()/sumOfWeights()); scale(m_hist_tleppt, crossSection()/sumOfWeights()); scale(m_hist_tlepy, crossSection()/sumOfWeights()); scale(m_hist_ttpt, crossSection()/sumOfWeights()); scale(m_hist_tty, crossSection()/sumOfWeights()); scale(m_hist_ttm, crossSection()/sumOfWeights()); scale(m_hist_njet, crossSection()/sumOfWeights()); m_hist_njet_ttm.scale(crossSection()/sumOfWeights(), false); m_hist_njet_thadpt.scale(crossSection()/sumOfWeights(), false); m_hist_njet_ttpt.scale(crossSection()/sumOfWeights(), false); m_hist_thady_thadpt.scale(crossSection()/sumOfWeights()); m_hist_ttm_tty.scale(crossSection()/sumOfWeights()); m_hist_thadpt_ttm.scale(crossSection()/sumOfWeights()); m_hist_jetspt.scale(crossSection()/sumOfWeights(), false); m_hist_jetseta.scale(crossSection()/sumOfWeights(), false); m_hist_jetsdrjets.scale(crossSection()/sumOfWeights(), false); m_hist_jetsdrtops.scale(crossSection()/sumOfWeights(), false); m_hist_njetspt.scale(crossSection()/sumOfWeights(), false); scale(m_nhist_thadpt, 1./m_nhist_thadpt->sumW(false)); scale(m_nhist_thady, 1./m_nhist_thady->sumW(false)); scale(m_nhist_tleppt, 1./m_nhist_tleppt->sumW(false)); scale(m_nhist_tlepy, 1./m_nhist_tlepy->sumW(false)); scale(m_nhist_ttpt, 1./m_nhist_ttpt->sumW(false)); scale(m_nhist_tty, 1./m_nhist_tty->sumW(false)); scale(m_nhist_ttm, 1./m_nhist_ttm->sumW(false)); scale(m_nhist_njet, 1./m_nhist_njet->sumW(false)); m_nhist_njet_ttm.norm(false); m_nhist_njet_thadpt.norm(false); m_nhist_njet_ttpt.norm(false); m_nhist_thady_thadpt.norm(true); m_nhist_ttm_tty.norm(true); m_nhist_thadpt_ttm.norm(true); m_nhist_jetspt.norm(false); m_nhist_jetseta.norm(false); m_nhist_jetsdrjets.norm(false); m_nhist_jetsdrtops.norm(false); } Histo1DPtr m_hist_thadpt; Histo1DPtr m_hist_thady; Histo1DPtr m_hist_tleppt; Histo1DPtr m_hist_tlepy; Histo1DPtr m_hist_ttpt; Histo1DPtr m_hist_tty; Histo1DPtr m_hist_ttm; Histo1DPtr m_hist_njet; Histo1DGroup m_hist_njet_ttm; Histo1DGroup m_hist_njet_thadpt; Histo1DGroup m_hist_njet_ttpt; Histo1DGroup m_hist_thady_thadpt; Histo1DGroup m_hist_ttm_tty; Histo1DGroup m_hist_thadpt_ttm; Histo1DGroup m_hist_jetspt; Histo1DGroup m_hist_jetseta; Histo1DGroup m_hist_jetsdrjets; Histo1DGroup m_hist_jetsdrtops; Histo1DGroup m_hist_njetspt; Histo1DPtr m_nhist_thadpt; Histo1DPtr m_nhist_thady; Histo1DPtr m_nhist_tleppt; Histo1DPtr m_nhist_tlepy; Histo1DPtr m_nhist_ttm; Histo1DPtr m_nhist_ttpt; Histo1DPtr m_nhist_tty; Histo1DPtr m_nhist_njet; Histo1DGroup m_nhist_njet_ttm; Histo1DGroup m_nhist_njet_thadpt; Histo1DGroup m_nhist_njet_ttpt; Histo1DGroup m_nhist_thady_thadpt; Histo1DGroup m_nhist_ttm_tty; Histo1DGroup m_nhist_thadpt_ttm; Histo1DGroup m_nhist_jetspt; Histo1DGroup m_nhist_jetseta; Histo1DGroup m_nhist_jetsdrjets; Histo1DGroup m_nhist_jetsdrtops; Scatter2DPtr m_hist_gap1; Scatter2DPtr m_hist_gap2; Particles m_neutrinos; Particles m_leptons; Particles m_vetoleptons; Jets m_bjets; Jets m_ljets; Particle m_thad; Particles m_thad_decay; Particle m_tlep; Particles m_tlep_decay; Particles m_tt_jets; Particles m_additionalobjects; Particles m_additionaljets; double m_jetptmin = 25.; double m_jetetamax = 2.4; double m_addjetptmin = 30.; double m_addjetetamax = 2.4; double m_jetdr = 0.4; double m_lepptmin = 30.; double m_lepetamax = 2.4; double m_vetolepptmin = 15.; double m_vetolepetamax = 2.4; double m_lepisomax = 0.35; double m_lepisodr = 0.4; double m_lepdressdr = 0.1; double m_phptmin = 15.; double m_phetamax = 2.4; double m_phisomax = 0.25; double m_phisodr = 0.4; }; DECLARE_RIVET_PLUGIN(CMS_2018_I1663958); CMS_2018_I1663958::Histo1DGroup::Histo1DGroup(CMS_2018_I1663958* an, const vector& hnames, const vector& xbinranges) : m_an(an), m_xbins(xbinranges) { for (const string& hname : hnames) { Histo1DPtr tmp; m_histos.push_back(an->book(tmp,hname)); } } void CMS_2018_I1663958::Histo1DGroup::fill(double x, double y, double w) { if (x < m_xbins[0]) {return;} if (x >= m_xbins[m_xbins.size()-1]) {return;} int xbin = upper_bound(m_xbins.begin(), m_xbins.end(), x) - m_xbins.begin()-1; m_histos[xbin]->fill(y, w); } void CMS_2018_I1663958::Histo1DGroup::scale(double s, bool diff) { for (size_t h = 0 ; h < m_histos.size() ; ++h) { double sc = s; if (diff) {sc /= m_xbins[h+1]-m_xbins[h];} m_an->scale(m_histos[h], sc); } } void CMS_2018_I1663958::Histo1DGroup::norm(bool diff) { double sum = 0.; for (Histo1DPtr hist : m_histos) {sum+=hist->sumW(false);} scale(1./sum, diff); } void CMS_2018_I1663958::Histo1DGroup::gapfractionfromjetpt(Scatter2DPtr hgap, int njet) { double total = m_histos[0]->integral(); int hn = njet+3; double totalj = m_histos[hn]->integral(); double acc = total-totalj; for (size_t nb = 0 ; nb < m_histos[hn]->numBins() ; ++nb) { double gf = acc/total; double bl = m_histos[njet+3]->bin(nb).xMin(); double bh = m_histos[njet+3]->bin(nb).xMax(); double bc = 0.5*(bh+bl); hgap->addPoint(bc, gf, bc-bl, bh-bc, 0., 0.); acc += m_histos[njet+3]->bin(nb).area(); } } } diff --git a/analyses/pluginLEP/OPAL_2004_I648738.cc b/analyses/pluginLEP/OPAL_2004_I648738.cc --- a/analyses/pluginLEP/OPAL_2004_I648738.cc +++ b/analyses/pluginLEP/OPAL_2004_I648738.cc @@ -1,118 +1,118 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" namespace Rivet { class OPAL_2004_I648738 : public Analysis { public: /// Constructor OPAL_2004_I648738() : Analysis("OPAL_2004_I648738"), _sumW(3), _histo_xE(3) { } /// @name Analysis methods //@{ void init() { declare(FinalState(), "FS"); declare(ChargedFinalState(), "CFS"); unsigned int ih=0; if (inRange(0.5*sqrtS()/GeV, 4.0, 9.0)) { ih = 1; } else if (inRange(0.5*sqrtS()/GeV, 9.0, 19.0)) { ih = 2; } else if (inRange(0.5*sqrtS()/GeV, 19.0, 30.0)) { ih = 3; } else if (inRange(0.5*sqrtS()/GeV, 45.5, 45.7)) { ih = 5; } else if (inRange(0.5*sqrtS()/GeV, 30.0, 70.0)) { ih = 4; } else if (inRange(0.5*sqrtS()/GeV, 91.5, 104.5)) { ih = 6; } assert(ih>0); // book the histograms book(_histo_xE[0], ih+5,1,1); book(_histo_xE[1], ih+5,1,2); if(ih<5) book(_histo_xE[2] ,ih+5,1,3); book(_sumW[0], "sumW_0"); book(_sumW[1], "sumW_1"); book(_sumW[2], "sumW_2"); } /// Perform the per-event analysis void analyze(const Event& event) { // find the initial quarks/gluons ParticleVector initial; for (const GenParticle* p : Rivet::particles(event.genEvent())) { const GenVertex* pv = p->production_vertex(); const PdgId pid = abs(p->pdg_id()); if(!( (pid>=1&&pid<=5) || pid ==21) ) continue; bool passed = false; for (const GenParticle* pp : particles_in(pv)) { const PdgId ppid = abs(pp->pdg_id()); passed = (ppid == PID::ELECTRON || ppid == PID::HIGGS || ppid == PID::ZBOSON || ppid == PID::GAMMA); if(passed) break; } if(passed) initial.push_back(Particle(*p)); } if(initial.size()!=2) { vetoEvent; } // type of event unsigned int itype=2; - if(initial[0].pdgId()==-initial[1].pdgId()) { - PdgId pid = abs(initial[0].pdgId()); + if(initial[0].pid()==-initial[1].pid()) { + PdgId pid = abs(initial[0].pid()); if(pid>=1&&pid<=4) itype=0; else itype=1; } assert(itype<_histo_xE.size()); // fill histograms _sumW[itype]->fill(2.); const Particles& chps = applyProjection(event, "CFS").particles(); for(const Particle& p : chps) { double xE = 2.*p.E()/sqrtS(); _histo_xE[itype]->fill(xE); } } /// Normalise histograms etc., after the run void finalize() { for(unsigned int ix=0;ix<_histo_xE.size();++ix) { if(_sumW[ix]->val()>0.) scale(_histo_xE[ix],1./ *_sumW[ix]); } } //@} private: vector _sumW; /// @name Histograms //@{ vector _histo_xE; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(OPAL_2004_I648738); } diff --git a/analyses/pluginLHCb/LHCB_2014_I1281685.cc b/analyses/pluginLHCb/LHCB_2014_I1281685.cc --- a/analyses/pluginLHCb/LHCB_2014_I1281685.cc +++ b/analyses/pluginLHCb/LHCB_2014_I1281685.cc @@ -1,1177 +1,1177 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" namespace Rivet { /// Charged particle multiplicities and densities in $pp$ collisions at $\sqrt{s} = 7$ TeV class LHCB_2014_I1281685 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor LHCB_2014_I1281685() : Analysis("LHCB_2014_I1281685"), _p_min(2.0), _pt_min(0.2), _eta_min(2.0), _eta_max(4.8), _maxlft(1.0e-11) { } //@} /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { fillMap(_partLftMap); // Projections declare(ChargedFinalState(_eta_min, _eta_max, _pt_min*GeV), "CFS"); // Book histograms book(_h_mult_total ,"d03-x01-y01", 50, 0.5, 50.5); book(_h_mult_eta[0] ,"d04-x01-y01", 21, -0.5, 20.5); //eta=[2.0,2.5] book(_h_mult_eta[1] ,"d04-x01-y02", 21, -0.5, 20.5); //eta=[2.5,3.0] book(_h_mult_eta[2] ,"d04-x01-y03", 21, -0.5, 20.5); //eta=[3.0,3.5] book(_h_mult_eta[3] ,"d04-x01-y04", 21, -0.5, 20.5); //eta=[3.5,4.0] book(_h_mult_eta[4] ,"d04-x01-y05", 21, -0.5, 20.5); //eta=[4.0,4.5] book(_h_mult_pt[0] ,"d05-x01-y01", 21, -0.5, 20.5); //pT=[0.2,0.3]GeV book(_h_mult_pt[1] ,"d05-x01-y02", 21, -0.5, 20.5); //pT=[0.3,0.4]GeV book(_h_mult_pt[2] ,"d05-x01-y03", 21, -0.5, 20.5); //pT=[0.4,0.6]GeV book(_h_mult_pt[3] ,"d05-x01-y04", 21, -0.5, 20.5); //pT=[0.6,1.0]GeV book(_h_mult_pt[4] ,"d05-x01-y05", 21, -0.5, 20.5); //pT=[1.0,2.0]GeV book(_h_dndeta ,"d01-x01-y01", 14, 2.0, 4.8); //eta=[2,4.8] book(_h_dndpt ,"d02-x01-y01", 18, 0.2, 2.0); //pT =[0,2]GeV // Counters book(_sumW, "TMP/sumW"); } /// Perform the per-event analysis void analyze(const Event& event) { // Variable to store multiplicities per event int LHCbcountAll = 0; //count particles fulfiling all requirements int LHCbcountEta[8] = {0,0,0,0,0,0,0,0}; //count per eta-bin int LHCbcountPt[7] = {0,0,0,0,0,0,0}; //count per pT-bin vector val_dNdEta; vector val_dNdPt; val_dNdEta.clear(); val_dNdPt.clear(); const ChargedFinalState& cfs = apply(event, "CFS"); for (const Particle& p : cfs.particles()) { - int id = p.pdgId(); + int id = p.pid(); // continue if particle is not a pion, kaon, proton, muon or electron if ( !( (abs(id) == 211) || (abs(id) == 321) || (abs(id) == 2212) || (abs(id) == 13) || (abs(id) == 11)) ) { continue; } const FourMomentum& qmom = p.momentum(); const double eta = p.momentum().eta(); const double pT = p.momentum().pT(); //minimum momentum if (qmom.p3().mod() < _p_min) continue; //minimum tr. momentum if (pT < _pt_min) continue; //eta range if ((eta < _eta_min) || (eta > _eta_max)) continue; /* Select only prompt particles via lifetime */ //Sum of all mother lifetimes (PDG lifetime) < 10ps double ancestors_sumlft = getAncestorSumLifetime(p); if( (ancestors_sumlft > _maxlft) || (ancestors_sumlft < 0) ) continue; //after all cuts; LHCbcountAll++; //count particles in whole kin. range //in eta bins if( eta >2.0 && eta <= 2.5) LHCbcountEta[0]++; if( eta >2.5 && eta <= 3.0) LHCbcountEta[1]++; if( eta >3.0 && eta <= 3.5) LHCbcountEta[2]++; if( eta >3.5 && eta <= 4.0) LHCbcountEta[3]++; if( eta >4.0 && eta <= 4.5) LHCbcountEta[4]++; if( eta >2.0 && eta <= 4.8) LHCbcountEta[5]++; //cross-check //in pT bins if( pT > 0.2 && pT <= 0.3) LHCbcountPt[0]++; if( pT > 0.3 && pT <= 0.4) LHCbcountPt[1]++; if( pT > 0.4 && pT <= 0.6) LHCbcountPt[2]++; if( pT > 0.6 && pT <= 1.0) LHCbcountPt[3]++; if( pT > 1.0 && pT <= 2.0) LHCbcountPt[4]++; if( pT > 0.2) LHCbcountPt[5]++; //cross-check //particle densities -> need proper normalization (finalize) val_dNdPt.push_back( pT ); val_dNdEta.push_back( eta ); }//end for // Fill histograms only, if at least 1 particle pre event was within the // kinematic range of the analysis! if (LHCbcountAll) { _sumW->fill(); _h_mult_total->fill(LHCbcountAll); _h_mult_eta[0]->fill(LHCbcountEta[0]); _h_mult_eta[1]->fill(LHCbcountEta[1]); _h_mult_eta[2]->fill(LHCbcountEta[2]); _h_mult_eta[3]->fill(LHCbcountEta[3]); _h_mult_eta[4]->fill(LHCbcountEta[4]); _h_mult_pt[0]->fill(LHCbcountPt[0]); _h_mult_pt[1]->fill(LHCbcountPt[1]); _h_mult_pt[2]->fill(LHCbcountPt[2]); _h_mult_pt[3]->fill(LHCbcountPt[3]); _h_mult_pt[4]->fill(LHCbcountPt[4]); for (size_t part = 0; part < val_dNdEta.size(); part++) _h_dndeta->fill(val_dNdEta[part]); for (size_t part = 0; part < val_dNdPt.size(); part++) _h_dndpt->fill(val_dNdPt[part]); } } /// Normalise histograms etc., after the run void finalize() { const double scalefactor = 1.0/_sumW->val(); // normalize multiplicity histograms by nEvents const double scale1k = 1000.; // to match '10^3' scale in reference histograms scale( _h_dndeta, scalefactor ); scale( _h_dndpt, scalefactor*0.1 ); //additional factor 0.1 for [0.1 GeV/c] scale( _h_mult_total, scalefactor*scale1k); _h_mult_eta[0]->scaleW( scalefactor*scale1k ); _h_mult_eta[1]->scaleW( scalefactor*scale1k ); _h_mult_eta[2]->scaleW( scalefactor*scale1k ); _h_mult_eta[3]->scaleW( scalefactor*scale1k ); _h_mult_eta[4]->scaleW( scalefactor*scale1k ); _h_mult_pt[0]->scaleW( scalefactor*scale1k ); _h_mult_pt[1]->scaleW( scalefactor*scale1k ); _h_mult_pt[2]->scaleW( scalefactor*scale1k ); _h_mult_pt[3]->scaleW( scalefactor*scale1k ); _h_mult_pt[4]->scaleW( scalefactor*scale1k ); } //@} private: // Get mean PDG lifetime for particle with PID double getLifetime(int pid) { double lft = 0.; map::iterator pPartLft = _partLftMap.find(pid); if (pPartLft != _partLftMap.end()) { lft = (*pPartLft).second; } else { // allow identifying missing life times only in debug mode MSG_DEBUG("Could not determine lifetime for particle with PID " << pid << "... Assume non-prompt particle"); lft = -1; } return lft; } // Get sum of all ancestor particles const double getAncestorSumLifetime(const Particle& p) { double lftSum = 0.; double plft = 0.; const GenParticle* part = p.genParticle(); if ( 0 == part ) return -1; const GenVertex* ivtx = part->production_vertex(); while(ivtx) { if (ivtx->particles_in_size() < 1) { lftSum = -1.; break; }; const GenVertex::particles_in_const_iterator iPart_invtx = ivtx->particles_in_const_begin(); part = (*iPart_invtx); if ( !(part) ) { lftSum = -1.; break; }; ivtx = part->production_vertex(); if ( (part->pdg_id() == 2212) || !(ivtx) ) break; // reached beam plft = getLifetime(part->pdg_id()); if (plft < 0.) { lftSum = -1.; break; }; lftSum += plft; } return (lftSum); } /// Hard-coded map linking PDG ID with PDG lifetime[s] (converted from ParticleTable.txt) bool fillMap(map& m) { // PDGID = LIFETIME m[22] = 1.000000e+016; m[-11] = 1.000000e+016; m[11] = 1.000000e+016; m[12] = 1.000000e+016; m[-13] = 2.197036e-006; m[13] = 2.197036e-006; m[111] = 8.438618e-017; m[211] = 2.603276e-008; m[-211] = 2.603276e-008; m[130] = 5.174624e-008; m[321] = 1.238405e-008; m[-321] = 1.238405e-008; m[2112] = 885.646128; m[2212] = 1.000000e+016; m[-2212] = 1.000000e+016; m[310] = 8.934603e-011; m[221] = 5.578070e-019; m[3122] = 2.631796e-010; m[3222] = 8.018178e-011; m[3212] = 7.395643e-020; m[3112] = 1.479129e-010; m[3322] = 2.899613e-010; m[3312] = 1.637344e-010; m[3334] = 8.207135e-011; m[-2112] = 885.646128; m[-3122] = 2.631796e-010; m[-3222] = 8.018178e-011; m[-3212] = 7.395643e-020; m[-3112] = 1.479129e-010; m[-3322] = 2.899613e-010; m[-3312] = 1.637344e-010; m[-3334] = 8.207135e-011; m[113] = 4.411610e-024; m[213] = 4.411610e-024; m[-213] = 4.411610e-024; m[223] = 7.798723e-023; m[333] = 1.545099e-022; m[323] = 1.295693e-023; m[-323] = 1.295693e-023; m[313] = 1.298249e-023; m[-313] = 1.298249e-023; m[20213] = 1.500000e-024; m[-20213] = 1.500000e-024; m[450000000] = 1.000000e+015; m[460000000] = 1.000000e+015; m[470000000] = 1.000000e+015; m[480000000] = 1.000000e+015; m[490000000] = 1.000000e+015; m[20022] = 1.000000e+016; m[-15] = 2.906014e-013; m[15] = 2.906014e-013; m[24] = 3.104775e-025; m[-24] = 3.104775e-025; m[23] = 2.637914e-025; m[411] = 1.051457e-012; m[-411] = 1.051457e-012; m[421] = 4.116399e-013; m[-421] = 4.116399e-013; m[431] = 4.904711e-013; m[-431] = 4.904711e-013; m[4122] = 1.994582e-013; m[-4122] = 1.994582e-013; m[443] = 7.565657e-021; m[413] = 6.856377e-021; m[-413] = 6.856377e-021; m[423] = 1.000003e-019; m[-423] = 1.000003e-019; m[433] = 1.000003e-019; m[-433] = 1.000003e-019; m[521] = 1.671000e-012; m[-521] = 1.671000e-012; m[511] = 1.536000e-012; m[-511] = 1.536000e-012; m[531] = 1.461000e-012; m[-531] = 1.461000e-012; m[541] = 4.600000e-013; m[-541] = 4.600000e-013; m[5122] = 1.229000e-012; m[-5122] = 1.229000e-012; m[4112] = 4.388081e-022; m[-4112] = 4.388081e-022; m[4212] = 3.999999e-022; m[-4212] = 3.999999e-022; m[4222] = 3.291060e-022; m[-4222] = 3.291060e-022; m[25] = 9.400000e-026; m[35] = 9.400000e-026; m[36] = 9.400000e-026; m[37] = 9.400000e-026; m[-37] = 9.400000e-026; m[4312] = 9.800002e-014; m[-4312] = 9.800002e-014; m[4322] = 3.500001e-013; m[-4322] = 3.500001e-013; m[4332] = 6.453061e-014; m[-4332] = 6.453061e-014; m[4132] = 9.824063e-014; m[-4132] = 9.824063e-014; m[4232] = 4.417532e-013; m[-4232] = 4.417532e-013; m[5222] = 1.000000e-019; m[-5222] = 1.000000e-019; m[5212] = 1.000000e-019; m[-5212] = 1.000000e-019; m[5112] = 1.000000e-019; m[-5112] = 1.000000e-019; m[5312] = 1.000000e-019; m[-5312] = 1.000000e-019; m[5322] = 1.000000e-019; m[-5322] = 1.000000e-019; m[5332] = 1.550000e-012; m[-5332] = 1.550000e-012; m[5132] = 1.390000e-012; m[-5132] = 1.390000e-012; m[5232] = 1.390000e-012; m[-5232] = 1.390000e-012; m[100443] = 2.194041e-021; m[331] = 3.258476e-021; m[441] = 4.113826e-023; m[10441] = 4.063038e-023; m[20443] = 7.154480e-022; m[445] = 3.164482e-022; m[9000111] = 1.149997e-023; m[9000211] = 1.149997e-023; m[-9000211] = 1.149997e-023; m[20113] = 1.500000e-024; m[115] = 6.151516e-024; m[215] = 6.151516e-024; m[-215] = 6.151516e-024; m[10323] = 7.313469e-024; m[-10323] = 7.313469e-024; m[10313] = 7.313469e-024; m[-10313] = 7.313469e-024; m[20323] = 3.782829e-024; m[-20323] = 3.782829e-024; m[20313] = 3.782829e-024; m[-20313] = 3.782829e-024; m[10321] = 2.238817e-024; m[-10321] = 2.238817e-024; m[10311] = 2.238817e-024; m[-10311] = 2.238817e-024; m[325] = 6.682357e-024; m[-325] = 6.682357e-024; m[315] = 6.038644e-024; m[-315] = 6.038644e-024; m[10411] = 4.380000e-024; m[20413] = 2.630000e-024; m[10413] = 3.290000e-023; m[-415] = 2.632849e-023; m[-10411] = 4.380000e-024; m[-20413] = 2.630000e-024; m[-10413] = 3.290000e-023; m[415] = 2.632849e-023; m[10421] = 4.380000e-024; m[20423] = 2.630000e-024; m[10423] = 3.482604e-023; m[-425] = 2.861792e-023; m[-10421] = 4.380000e-024; m[-20423] = 2.630000e-024; m[-10423] = 3.482604e-023; m[425] = 2.861792e-023; m[10431] = 6.582100e-022; m[20433] = 6.582100e-022; m[10433] = 6.582100e-022; m[435] = 4.388100e-023; m[-10431] = 6.582100e-022; m[-20433] = 6.582100e-022; m[-10433] = 6.582100e-022; m[-435] = 4.388100e-023; m[2224] = 5.485102e-024; m[2214] = 5.485102e-024; m[2114] = 5.485102e-024; m[1114] = 5.485102e-024; m[-2224] = 5.485102e-024; m[-2214] = 5.485102e-024; m[-2114] = 5.485102e-024; m[-1114] = 5.485102e-024; m[-523] = 1.000019e-019; m[523] = 1.000019e-019; m[513] = 1.000019e-019; m[-513] = 1.000019e-019; m[533] = 1.000000e-019; m[-533] = 1.000000e-019; m[10521] = 4.390000e-024; m[20523] = 2.630000e-024; m[10523] = 1.650000e-023; m[525] = 1.310000e-023; m[-10521] = 4.390000e-024; m[-20523] = 2.630000e-024; m[-10523] = 1.650000e-023; m[-525] = 1.310000e-023; m[10511] = 4.390000e-024; m[20513] = 2.630000e-024; m[10513] = 1.650000e-023; m[515] = 1.310000e-023; m[-10511] = 4.390000e-024; m[-20513] = 2.630000e-024; m[-10513] = 1.650000e-023; m[-515] = 1.310000e-023; m[10531] = 4.390000e-024; m[20533] = 2.630000e-024; m[10533] = 1.650000e-023; m[535] = 1.310000e-023; m[-10531] = 4.390000e-024; m[-20533] = 2.630000e-024; m[-10533] = 1.650000e-023; m[-535] = 1.310000e-023; m[14] = 1.000000e+016; m[-14] = 1.000000e+016; m[-12] = 1.000000e+016; m[1] = 0.000000e+000; m[-1] = 0.000000e+000; m[2] = 0.000000e+000; m[-2] = 0.000000e+000; m[3] = 0.000000e+000; m[-3] = 0.000000e+000; m[4] = 0.000000e+000; m[-4] = 0.000000e+000; m[5] = 0.000000e+000; m[-5] = 0.000000e+000; m[6] = 4.707703e-025; m[-6] = 4.707703e-025; m[7] = 0.000000e+000; m[-7] = 0.000000e+000; m[8] = 0.000000e+000; m[-8] = 0.000000e+000; m[16] = 1.000000e+016; m[-16] = 1.000000e+016; m[17] = 0.000000e+000; m[-17] = 0.000000e+000; m[18] = 0.000000e+000; m[-18] = 0.000000e+000; m[21] = 0.000000e+000; m[32] = 0.000000e+000; m[33] = 0.000000e+000; m[34] = 0.000000e+000; m[-34] = 0.000000e+000; m[39] = 0.000000e+000; m[41] = 0.000000e+000; m[-41] = 0.000000e+000; m[42] = 0.000000e+000; m[-42] = 0.000000e+000; m[43] = 0.000000e+000; m[44] = 0.000000e+000; m[-44] = 0.000000e+000; m[81] = 0.000000e+000; m[82] = 0.000000e+000; m[-82] = 0.000000e+000; m[83] = 0.000000e+000; m[84] = 3.335641e-013; m[-84] = 3.335641e-013; m[85] = 1.290893e-012; m[-85] = 1.290893e-012; m[86] = 0.000000e+000; m[-86] = 0.000000e+000; m[87] = 0.000000e+000; m[-87] = 0.000000e+000; m[88] = 0.000000e+000; m[90] = 0.000000e+000; m[91] = 0.000000e+000; m[92] = 0.000000e+000; m[93] = 0.000000e+000; m[94] = 0.000000e+000; m[95] = 0.000000e+000; m[96] = 0.000000e+000; m[97] = 0.000000e+000; m[98] = 0.000000e+000; m[99] = 0.000000e+000; m[117] = 4.088275e-024; m[119] = 1.828367e-024; m[217] = 4.088275e-024; m[-217] = 4.088275e-024; m[219] = 1.828367e-024; m[-219] = 1.828367e-024; m[225] = 3.555982e-024; m[227] = 3.917930e-024; m[229] = 3.392846e-024; m[311] = 1.000000e+016; m[-311] = 1.000000e+016; m[317] = 4.139699e-024; m[-317] = 4.139699e-024; m[319] = 3.324304e-024; m[-319] = 3.324304e-024; m[327] = 4.139699e-024; m[-327] = 4.139699e-024; m[329] = 3.324304e-024; m[-329] = 3.324304e-024; m[335] = 8.660687e-024; m[337] = 7.565657e-024; m[543] = 0.000000e+000; m[-543] = 0.000000e+000; m[545] = 0.000000e+000; m[-545] = 0.000000e+000; m[551] = 0.000000e+000; m[553] = 1.253738e-020; m[555] = 1.000000e+016; m[557] = 0.000000e+000; m[-450000000] = 0.000000e+000; m[-490000000] = 0.000000e+000; m[-460000000] = 0.000000e+000; m[-470000000] = 0.000000e+000; m[1103] = 0.000000e+000; m[-1103] = 0.000000e+000; m[1112] = 4.388081e-024; m[-1112] = 4.388081e-024; m[1116] = 1.880606e-024; m[-1116] = 1.880606e-024; m[1118] = 2.194041e-024; m[-1118] = 2.194041e-024; m[1212] = 4.388081e-024; m[-1212] = 4.388081e-024; m[1214] = 5.485102e-024; m[-1214] = 5.485102e-024; m[1216] = 1.880606e-024; m[-1216] = 1.880606e-024; m[1218] = 1.462694e-024; m[-1218] = 1.462694e-024; m[2101] = 0.000000e+000; m[-2101] = 0.000000e+000; m[2103] = 0.000000e+000; m[-2103] = 0.000000e+000; m[2116] = 4.388081e-024; m[-2116] = 4.388081e-024; m[2118] = 2.194041e-024; m[-2118] = 2.194041e-024; m[2122] = 4.388081e-024; m[-2122] = 4.388081e-024; m[2124] = 5.485102e-024; m[-2124] = 5.485102e-024; m[2126] = 1.880606e-024; m[-2126] = 1.880606e-024; m[2128] = 1.462694e-024; m[-2128] = 1.462694e-024; m[2203] = 0.000000e+000; m[-2203] = 0.000000e+000; m[2216] = 4.388081e-024; m[-2216] = 4.388081e-024; m[2218] = 2.194041e-024; m[-2218] = 2.194041e-024; m[2222] = 4.388081e-024; m[-2222] = 4.388081e-024; m[2226] = 1.880606e-024; m[-2226] = 1.880606e-024; m[2228] = 2.194041e-024; m[-2228] = 2.194041e-024; m[3101] = 0.000000e+000; m[-3101] = 0.000000e+000; m[3103] = 0.000000e+000; m[-3103] = 0.000000e+000; m[3114] = 1.670589e-023; m[-3114] = 1.670589e-023; m[3116] = 5.485102e-024; m[-3116] = 5.485102e-024; m[3118] = 3.656734e-024; m[-3118] = 3.656734e-024; m[3124] = 4.219309e-023; m[-3124] = 4.219309e-023; m[3126] = 8.227653e-024; m[-3126] = 8.227653e-024; m[3128] = 3.291061e-024; m[-3128] = 3.291061e-024; m[3201] = 0.000000e+000; m[-3201] = 0.000000e+000; m[3203] = 0.000000e+000; m[-3203] = 0.000000e+000; m[3214] = 1.828367e-023; m[-3214] = 1.828367e-023; m[3216] = 5.485102e-024; m[-3216] = 5.485102e-024; m[3218] = 3.656734e-024; m[-3218] = 3.656734e-024; m[3224] = 1.838582e-023; m[-3224] = 1.838582e-023; m[3226] = 5.485102e-024; m[-3226] = 5.485102e-024; m[3228] = 3.656734e-024; m[-3228] = 3.656734e-024; m[3303] = 0.000000e+000; m[-3303] = 0.000000e+000; m[3314] = 6.648608e-023; m[-3314] = 6.648608e-023; m[3324] = 7.233101e-023; m[-3324] = 7.233101e-023; m[4101] = 0.000000e+000; m[-4101] = 0.000000e+000; m[4103] = 0.000000e+000; m[-4103] = 0.000000e+000; m[4114] = 0.000000e+000; m[-4114] = 0.000000e+000; m[4201] = 0.000000e+000; m[-4201] = 0.000000e+000; m[4203] = 0.000000e+000; m[-4203] = 0.000000e+000; m[4214] = 3.291061e-022; m[-4214] = 3.291061e-022; m[4224] = 0.000000e+000; m[-4224] = 0.000000e+000; m[4301] = 0.000000e+000; m[-4301] = 0.000000e+000; m[4303] = 0.000000e+000; m[-4303] = 0.000000e+000; m[4314] = 0.000000e+000; m[-4314] = 0.000000e+000; m[4324] = 0.000000e+000; m[-4324] = 0.000000e+000; m[4334] = 0.000000e+000; m[-4334] = 0.000000e+000; m[4403] = 0.000000e+000; m[-4403] = 0.000000e+000; m[4412] = 3.335641e-013; m[-4412] = 3.335641e-013; m[4414] = 3.335641e-013; m[-4414] = 3.335641e-013; m[4422] = 3.335641e-013; m[-4422] = 3.335641e-013; m[4424] = 3.335641e-013; m[-4424] = 3.335641e-013; m[4432] = 3.335641e-013; m[-4432] = 3.335641e-013; m[4434] = 3.335641e-013; m[-4434] = 3.335641e-013; m[4444] = 3.335641e-013; m[-4444] = 3.335641e-013; m[5101] = 0.000000e+000; m[-5101] = 0.000000e+000; m[5103] = 0.000000e+000; m[-5103] = 0.000000e+000; m[5114] = 0.000000e+000; m[-5114] = 0.000000e+000; m[5142] = 1.290893e-012; m[-5142] = 1.290893e-012; m[5201] = 0.000000e+000; m[-5201] = 0.000000e+000; m[5203] = 0.000000e+000; m[-5203] = 0.000000e+000; m[5214] = 0.000000e+000; m[-5214] = 0.000000e+000; m[5224] = 0.000000e+000; m[-5224] = 0.000000e+000; m[5242] = 1.290893e-012; m[-5242] = 1.290893e-012; m[5301] = 0.000000e+000; m[-5301] = 0.000000e+000; m[5303] = 0.000000e+000; m[-5303] = 0.000000e+000; m[5314] = 0.000000e+000; m[-5314] = 0.000000e+000; m[5324] = 0.000000e+000; m[-5324] = 0.000000e+000; m[5334] = 0.000000e+000; m[-5334] = 0.000000e+000; m[5342] = 1.290893e-012; m[-5342] = 1.290893e-012; m[5401] = 0.000000e+000; m[-5401] = 0.000000e+000; m[5403] = 0.000000e+000; m[-5403] = 0.000000e+000; m[5412] = 1.290893e-012; m[-5412] = 1.290893e-012; m[5414] = 1.290893e-012; m[-5414] = 1.290893e-012; m[5422] = 1.290893e-012; m[-5422] = 1.290893e-012; m[5424] = 1.290893e-012; m[-5424] = 1.290893e-012; m[5432] = 1.290893e-012; m[-5432] = 1.290893e-012; m[5434] = 1.290893e-012; m[-5434] = 1.290893e-012; m[5442] = 1.290893e-012; m[-5442] = 1.290893e-012; m[5444] = 1.290893e-012; m[-5444] = 1.290893e-012; m[5503] = 0.000000e+000; m[-5503] = 0.000000e+000; m[5512] = 1.290893e-012; m[-5512] = 1.290893e-012; m[5514] = 1.290893e-012; m[-5514] = 1.290893e-012; m[5522] = 1.290893e-012; m[-5522] = 1.290893e-012; m[5524] = 1.290893e-012; m[-5524] = 1.290893e-012; m[5532] = 1.290893e-012; m[-5532] = 1.290893e-012; m[5534] = 1.290893e-012; m[-5534] = 1.290893e-012; m[5542] = 1.290893e-012; m[-5542] = 1.290893e-012; m[5544] = 1.290893e-012; m[-5544] = 1.290893e-012; m[5554] = 1.290893e-012; m[-5554] = 1.290893e-012; m[10022] = 0.000000e+000; m[10111] = 2.483820e-024; m[10113] = 4.635297e-024; m[10115] = 2.541360e-024; m[10211] = 2.483820e-024; m[-10211] = 2.483820e-024; m[10213] = 4.635297e-024; m[-10213] = 4.635297e-024; m[10215] = 2.541360e-024; m[-10215] = 2.541360e-024; m[9010221] = 1.316424e-023; m[10223] = 1.828367e-024; m[10225] = 0.000000e+000; m[10315] = 3.538775e-024; m[-10315] = 3.538775e-024; m[10325] = 3.538775e-024; m[-10325] = 3.538775e-024; m[10331] = 5.265698e-024; m[10333] = 0.000000e+000; m[10335] = 0.000000e+000; m[10443] = 0.000000e+000; m[10541] = 0.000000e+000; m[-10541] = 0.000000e+000; m[10543] = 0.000000e+000; m[-10543] = 0.000000e+000; m[10551] = 1.000000e+016; m[10553] = 0.000000e+000; m[10555] = 0.000000e+000; m[11112] = 0.000000e+000; m[-11112] = 0.000000e+000; m[11114] = 2.194041e-024; m[-11114] = 2.194041e-024; m[11116] = 1.880606e-024; m[-11116] = 1.880606e-024; m[11212] = 1.880606e-024; m[-11212] = 1.880606e-024; m[11216] = 0.000000e+000; m[-11216] = 0.000000e+000; m[12112] = 1.880606e-024; m[-12112] = 1.880606e-024; m[12114] = 2.194041e-024; m[-12114] = 2.194041e-024; m[12116] = 5.063171e-024; m[-12116] = 5.063171e-024; m[12118] = 0.000000e+000; m[-12118] = 0.000000e+000; m[12122] = 0.000000e+000; m[-12122] = 0.000000e+000; m[12126] = 1.880606e-024; m[-12126] = 1.880606e-024; m[12212] = 1.880606e-024; m[-12212] = 1.880606e-024; m[12214] = 2.194041e-024; m[-12214] = 2.194041e-024; m[12216] = 5.063171e-024; m[-12216] = 5.063171e-024; m[12218] = 0.000000e+000; m[-12218] = 0.000000e+000; m[12222] = 0.000000e+000; m[-12222] = 0.000000e+000; m[12224] = 2.194041e-024; m[-12224] = 2.194041e-024; m[12226] = 1.880606e-024; m[-12226] = 1.880606e-024; m[13112] = 6.582122e-024; m[-13112] = 6.582122e-024; m[13114] = 1.097020e-023; m[-13114] = 1.097020e-023; m[13116] = 5.485102e-024; m[-13116] = 5.485102e-024; m[13122] = 1.316424e-023; m[-13122] = 1.316424e-023; m[13124] = 1.097020e-023; m[-13124] = 1.097020e-023; m[13126] = 6.928549e-024; m[-13126] = 6.928549e-024; m[13212] = 6.582122e-024; m[-13212] = 6.582122e-024; m[13214] = 1.097020e-023; m[-13214] = 1.097020e-023; m[13216] = 5.485102e-024; m[-13216] = 5.485102e-024; m[13222] = 6.582122e-024; m[-13222] = 6.582122e-024; m[13224] = 1.097020e-023; m[-13224] = 1.097020e-023; m[13226] = 5.485102e-024; m[-13226] = 5.485102e-024; m[13314] = 2.742551e-023; m[-13314] = 2.742551e-023; m[13316] = 0.000000e+000; m[-13316] = 0.000000e+000; m[13324] = 2.742551e-023; m[-13324] = 2.742551e-023; m[13326] = 0.000000e+000; m[-13326] = 0.000000e+000; m[14122] = 1.828367e-022; m[-14122] = 1.828367e-022; m[14124] = 0.000000e+000; m[-14124] = 0.000000e+000; m[10221] = 2.194040e-024; m[20223] = 2.742551e-023; m[20315] = 2.384827e-024; m[-20315] = 2.384827e-024; m[20325] = 2.384827e-024; m[-20325] = 2.384827e-024; m[20333] = 1.185968e-023; m[20543] = 0.000000e+000; m[-20543] = 0.000000e+000; m[20553] = 1.000000e+016; m[20555] = 0.000000e+000; m[21112] = 2.632849e-024; m[-21112] = 2.632849e-024; m[21114] = 3.291061e-024; m[-21114] = 3.291061e-024; m[21212] = 2.632849e-024; m[-21212] = 2.632849e-024; m[21214] = 6.582122e-024; m[-21214] = 6.582122e-024; m[22112] = 4.388081e-024; m[-22112] = 4.388081e-024; m[22114] = 3.291061e-024; m[-22114] = 3.291061e-024; m[22122] = 2.632849e-024; m[-22122] = 2.632849e-024; m[22124] = 6.582122e-024; m[-22124] = 6.582122e-024; m[22212] = 4.388081e-024; m[-22212] = 4.388081e-024; m[22214] = 3.291061e-024; m[-22214] = 3.291061e-024; m[22222] = 2.632849e-024; m[-22222] = 2.632849e-024; m[22224] = 3.291061e-024; m[-22224] = 3.291061e-024; m[23112] = 7.313469e-024; m[-23112] = 7.313469e-024; m[23114] = 2.991874e-024; m[-23114] = 2.991874e-024; m[23122] = 4.388081e-024; m[-23122] = 4.388081e-024; m[23124] = 6.582122e-024; m[-23124] = 6.582122e-024; m[23126] = 3.291061e-024; m[-23126] = 3.291061e-024; m[23212] = 7.313469e-024; m[-23212] = 7.313469e-024; m[23214] = 2.991874e-024; m[-23214] = 2.991874e-024; m[23222] = 7.313469e-024; m[-23222] = 7.313469e-024; m[23224] = 2.991874e-024; m[-23224] = 2.991874e-024; m[23314] = 0.000000e+000; m[-23314] = 0.000000e+000; m[23324] = 0.000000e+000; m[-23324] = 0.000000e+000; m[30113] = 2.742551e-024; m[30213] = 2.742551e-024; m[-30213] = 2.742551e-024; m[30223] = 2.991874e-024; m[30313] = 2.056913e-024; m[-30313] = 2.056913e-024; m[30323] = 2.056913e-024; m[-30323] = 2.056913e-024; m[30343] = 0.000000e+000; m[-30343] = 0.000000e+000; m[30353] = 0.000000e+000; m[-30353] = 0.000000e+000; m[30363] = 0.000000e+000; m[-30363] = 0.000000e+000; m[30411] = 0.000000e+000; m[-30411] = 0.000000e+000; m[30413] = 0.000000e+000; m[-30413] = 0.000000e+000; m[30421] = 0.000000e+000; m[-30421] = 0.000000e+000; m[30423] = 0.000000e+000; m[-30423] = 0.000000e+000; m[30443] = 2.789035e-023; m[30553] = 0.000000e+000; m[31114] = 1.880606e-024; m[-31114] = 1.880606e-024; m[31214] = 4.388081e-024; m[-31214] = 4.388081e-024; m[32112] = 4.388081e-024; m[-32112] = 4.388081e-024; m[32114] = 1.880606e-024; m[-32114] = 1.880606e-024; m[32124] = 4.388081e-024; m[-32124] = 4.388081e-024; m[32212] = 4.388081e-024; m[-32212] = 4.388081e-024; m[32214] = 1.880606e-024; m[-32214] = 1.880606e-024; m[32224] = 1.880606e-024; m[-32224] = 1.880606e-024; m[33122] = 1.880606e-023; m[-33122] = 1.880606e-023; m[33314] = 0.000000e+000; m[-33314] = 0.000000e+000; m[33324] = 0.000000e+000; m[-33324] = 0.000000e+000; m[41214] = 0.000000e+000; m[-41214] = 0.000000e+000; m[42112] = 6.582122e-024; m[-42112] = 6.582122e-024; m[42124] = 0.000000e+000; m[-42124] = 0.000000e+000; m[42212] = 6.582122e-024; m[-42212] = 6.582122e-024; m[43122] = 2.194041e-024; m[-43122] = 2.194041e-024; m[52114] = 0.000000e+000; m[-52114] = 0.000000e+000; m[52214] = 0.000000e+000; m[-52214] = 0.000000e+000; m[53122] = 4.388081e-024; m[-53122] = 4.388081e-024; m[100111] = 1.645531e-024; m[100113] = 2.123265e-024; m[100211] = 1.645531e-024; m[-100211] = 1.645531e-024; m[100213] = 2.123265e-024; m[-100213] = 2.123265e-024; m[100221] = 1.196749e-023; m[100223] = 3.871836e-024; m[100225] = 0.000000e+000; m[100311] = 0.000000e+000; m[-100311] = 0.000000e+000; m[100313] = 2.837122e-024; m[-100313] = 2.837122e-024; m[100315] = 0.000000e+000; m[-100315] = 0.000000e+000; m[100321] = 0.000000e+000; m[-100321] = 0.000000e+000; m[100323] = 2.837122e-024; m[-100323] = 2.837122e-024; m[100325] = 0.000000e+000; m[-100325] = 0.000000e+000; m[100331] = 0.000000e+000; m[100333] = 4.388081e-024; m[100335] = 3.291061e-024; m[100441] = 0.000000e+000; m[100551] = 0.000000e+000; m[100553] = 1.495937e-020; m[100555] = 1.000000e+016; m[100557] = 0.000000e+000; m[110551] = 1.000000e+016; m[110553] = 0.000000e+000; m[110555] = 0.000000e+000; m[120553] = 1.000000e+016; m[120555] = 0.000000e+000; m[130553] = 0.000000e+000; m[200111] = 3.134344e-024; m[200211] = 3.134344e-024; m[-200211] = 3.134344e-024; m[200551] = 0.000000e+000; m[200553] = 2.502708e-020; m[200555] = 0.000000e+000; m[210551] = 0.000000e+000; m[210553] = 0.000000e+000; m[220553] = 0.000000e+000; m[300553] = 4.701516e-023; m[9000221] = 0.000000e+000; m[9000443] = 1.265793e-023; m[9000553] = 5.983747e-024; m[9010443] = 8.438618e-024; m[9010553] = 8.331800e-024; m[9020221] = 6.038644e-024; m[9020443] = 1.530726e-023; m[9060225] = 4.388081e-024; m[9070225] = 2.056913e-024; m[1000001] = 0.000000e+000; m[-1000001] = 0.000000e+000; m[1000002] = 0.000000e+000; m[-1000002] = 0.000000e+000; m[1000003] = 0.000000e+000; m[-1000003] = 0.000000e+000; m[1000004] = 0.000000e+000; m[-1000004] = 0.000000e+000; m[1000005] = 0.000000e+000; m[-1000005] = 0.000000e+000; m[1000006] = 0.000000e+000; m[-1000006] = 0.000000e+000; m[1000011] = 0.000000e+000; m[-1000011] = 0.000000e+000; m[1000012] = 0.000000e+000; m[-1000012] = 0.000000e+000; m[1000013] = 0.000000e+000; m[-1000013] = 0.000000e+000; m[1000014] = 0.000000e+000; m[-1000014] = 0.000000e+000; m[1000015] = 0.000000e+000; m[-1000015] = 0.000000e+000; m[1000016] = 0.000000e+000; m[-1000016] = 0.000000e+000; m[1000021] = 0.000000e+000; m[1000022] = 0.000000e+000; m[1000023] = 0.000000e+000; m[1000024] = 0.000000e+000; m[-1000024] = 0.000000e+000; m[1000025] = 0.000000e+000; m[1000035] = 0.000000e+000; m[1000037] = 0.000000e+000; m[-1000037] = 0.000000e+000; m[1000039] = 0.000000e+000; m[2000001] = 0.000000e+000; m[-2000001] = 0.000000e+000; m[2000002] = 0.000000e+000; m[-2000002] = 0.000000e+000; m[2000003] = 0.000000e+000; m[-2000003] = 0.000000e+000; m[2000004] = 0.000000e+000; m[-2000004] = 0.000000e+000; m[2000005] = 0.000000e+000; m[-2000005] = 0.000000e+000; m[2000006] = 0.000000e+000; m[-2000006] = 0.000000e+000; m[2000011] = 0.000000e+000; m[-2000011] = 0.000000e+000; m[2000012] = 0.000000e+000; m[-2000012] = 0.000000e+000; m[2000013] = 0.000000e+000; m[-2000013] = 0.000000e+000; m[2000014] = 0.000000e+000; m[-2000014] = 0.000000e+000; m[2000015] = 0.000000e+000; m[-2000015] = 0.000000e+000; m[2000016] = 0.000000e+000; m[-2000016] = 0.000000e+000; m[3000111] = 0.000000e+000; m[3000113] = 0.000000e+000; m[3000211] = 0.000000e+000; m[-3000211] = 0.000000e+000; m[3000213] = 0.000000e+000; m[-3000213] = 0.000000e+000; m[3000221] = 0.000000e+000; m[3000223] = 0.000000e+000; m[3000331] = 0.000000e+000; m[3100021] = 0.000000e+000; m[3100111] = 0.000000e+000; m[3100113] = 0.000000e+000; m[3200111] = 0.000000e+000; m[3200113] = 0.000000e+000; m[3300113] = 0.000000e+000; m[3400113] = 0.000000e+000; m[4000001] = 0.000000e+000; m[-4000001] = 0.000000e+000; m[4000002] = 0.000000e+000; m[-4000002] = 0.000000e+000; m[4000011] = 0.000000e+000; m[-4000011] = 0.000000e+000; m[4000012] = 0.000000e+000; m[-4000012] = 0.000000e+000; m[5000039] = 0.000000e+000; m[9900012] = 0.000000e+000; m[9900014] = 0.000000e+000; m[9900016] = 0.000000e+000; m[9900023] = 0.000000e+000; m[9900024] = 0.000000e+000; m[-9900024] = 0.000000e+000; m[9900041] = 0.000000e+000; m[-9900041] = 0.000000e+000; m[9900042] = 0.000000e+000; m[-9900042] = 0.000000e+000; m[1027013000] = 0.000000e+000; m[1012006000] = 0.000000e+000; m[1063029000] = 0.000000e+000; m[1014007000] = 0.000000e+000; m[1016008000] = 0.000000e+000; m[1028014000] = 0.000000e+000; m[1065029000] = 0.000000e+000; m[1009004000] = 0.000000e+000; m[1019009000] = 0.000000e+000; m[1056026000] = 0.000000e+000; m[1207082000] = 0.000000e+000; m[1208082000] = 0.000000e+000; m[1029014000] = 0.000000e+000; m[1206082000] = 0.000000e+000; m[1054026000] = 0.000000e+000; m[1018008000] = 0.000000e+000; m[1030014000] = 0.000000e+000; m[1057026000] = 0.000000e+000; m[1204082000] = 0.000000e+000; m[-99000000] = 0.000000e+000; m[1028013000] = 0.000000e+000; m[1040018000] = 0.000000e+000; m[1011005000] = 0.000000e+000; m[1012005000] = 0.000000e+000; m[1013006000] = 0.000000e+000; m[1014006000] = 0.000000e+000; m[1052024000] = 0.000000e+000; m[1024012000] = 0.000000e+000; m[1026012000] = 0.000000e+000; m[1027012000] = 0.000000e+000; m[1015007000] = 0.000000e+000; m[1022010000] = 0.000000e+000; m[1058028000] = 0.000000e+000; m[1060028000] = 0.000000e+000; m[1062028000] = 0.000000e+000; m[1064028000] = 0.000000e+000; m[1007003000] = 0.000000e+000; m[1025012000] = 0.000000e+000; m[1053024000] = 0.000000e+000; m[1055025000] = 0.000000e+000; m[1008004000] = 0.000000e+000; m[1010004000] = 0.000000e+000; m[1010005000] = 0.000000e+000; m[1016007000] = 0.000000e+000; m[1017008000] = 0.000000e+000; m[1019008000] = 0.000000e+000; m[1023010000] = 0.000000e+000; m[1024011000] = 0.000000e+000; m[1031015000] = 0.000000e+000; m[1039017000] = 0.000000e+000; m[1040017000] = 0.000000e+000; m[1036018000] = 0.000000e+000; m[1050024000] = 0.000000e+000; m[1054024000] = 0.000000e+000; m[1059026000] = 0.000000e+000; m[1061028000] = 0.000000e+000; m[1063028000] = 0.000000e+000; m[1092042000] = 0.000000e+000; m[1095042000] = 0.000000e+000; m[1096042000] = 0.000000e+000; m[1097042000] = 0.000000e+000; m[1098042000] = 0.000000e+000; m[1100042000] = 0.000000e+000; m[1108046000] = 0.000000e+000; // Added by hand: m[9902210] = 0.000000e+000; //diffractive p-state -> assume no lifetime return true; } private: /// @name Histograms //@{ Histo1DPtr _h_mult_total; // full kinematic range Histo1DPtr _h_mult_eta[5]; // in eta bins Histo1DPtr _h_mult_pt[5]; // in pT bins Histo1DPtr _h_dndeta; // density dn/deta Histo1DPtr _h_dndpt; // density dn/dpT //@} /// @name Private variables double _p_min; double _pt_min; double _eta_min; double _eta_max; double _maxlft; /// Count selected events CounterPtr _sumW; map _partLftMap; // Map }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(LHCB_2014_I1281685); } diff --git a/analyses/pluginLHCb/LHCB_2015_I1333223.cc b/analyses/pluginLHCb/LHCB_2015_I1333223.cc --- a/analyses/pluginLHCb/LHCB_2015_I1333223.cc +++ b/analyses/pluginLHCb/LHCB_2015_I1333223.cc @@ -1,112 +1,112 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Math/Units.hh" #include using namespace std; namespace Rivet { class LHCB_2015_I1333223 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor LHCB_2015_I1333223() : Analysis("LHCB_2015_I1333223") { } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Charged particles declare(ChargedFinalState(Cuts::eta> 2.0 && Cuts::eta <4.5 && Cuts::pT >0.2*GeV), "CFS"); // Reproducing only measurement for prompt charged particles book(_hInelasticXs ,1, 1, 1); } /// Perform the per-event analysis void analyze(const Event& event) { const ChargedFinalState &cfs = apply (event, "CFS"); // eliminate non-inelastic events and empty events in LHCb if (cfs.particles().size() == 0) vetoEvent; // See if this event has at least one prompt particle for (const Particle &myp : cfs.particles()){ double dPV = getPVDCA(myp); // if IP > 200 microns the particle is not considered prompt if ((dPV < 0.) || (dPV > 0.2 * millimeter)) { - MSG_DEBUG(" Vetoing " << myp.pdgId() << " at " << dPV); + MSG_DEBUG(" Vetoing " << myp.pid() << " at " << dPV); continue; } // histo gets filled only for inelastic events (at least one prompt charged particle) _hInelasticXs->fill(sqrtS()); break; } //end loop on particles } /// Normalise histograms etc., after the run void finalize() { scale(_hInelasticXs, crossSection()/sumOfWeights()/millibarn); } //@} private: /* * Compute Distance of Closest Approach in z range for one particle. * Assuming length unit is mm. * Returns -1. if unable to compute the DCA to PV. */ double getPVDCA(const Particle& p) { const HepMC::GenVertex* vtx = p.genParticle()->production_vertex(); if ( 0 == vtx ) return -1.; // Unit vector of particle's MOMENTUM three vector const Vector3 u = p.momentum().p3().unit(); // The interaction point is always at (0, 0,0,0) hence the // vector pointing from the PV to the particle production vertex is: Vector3 d(vtx->position().x(), vtx->position().y(), vtx->position().z()); // Subtract projection of d onto u from d double proj = d.dot(u); d -= (u * proj); // d should be orthogonal to u and it's length give the distance of // closest approach return d.mod(); } /// @name Histograms //@{ Histo1DPtr _hInelasticXs; //@} // }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(LHCB_2015_I1333223); } diff --git a/configure.ac b/configure.ac --- a/configure.ac +++ b/configure.ac @@ -1,323 +1,322 @@ ## Process this file with autoconf to produce a configure script. AC_PREREQ(2.59) AC_INIT([Rivet],[3.0.0-pre],[rivet@projects.hepforge.org],[Rivet]) ## Check and block installation into the src/build dir if test "$prefix" = "$PWD"; then AC_MSG_ERROR([Installation into the build directory is not supported: use a different --prefix argument]) fi ## Force default prefix to have a path value rather than NONE if test "$prefix" = "NONE"; then prefix=/usr/local fi AC_CONFIG_SRCDIR([src/Core/Analysis.cc]) AC_CONFIG_HEADERS([include/Rivet/Config/DummyConfig.hh include/Rivet/Config/RivetConfig.hh]) AM_INIT_AUTOMAKE([dist-bzip2 -Wall 1.10]) m4_ifdef([AM_SILENT_RULES], [AM_SILENT_RULES([yes])]) m4_ifdef([AM_PROG_AR], [AM_PROG_AR]) AC_CONFIG_MACRO_DIR([m4]) AC_SUBST(LT_OBJDIR) ## Package-specific #defines AC_DEFINE_UNQUOTED(RIVET_VERSION, "$PACKAGE_VERSION", "Rivet version string") AC_DEFINE_UNQUOTED(RIVET_NAME, "$PACKAGE_NAME", "Rivet name string") AC_DEFINE_UNQUOTED(RIVET_STRING, "$PACKAGE_STRING", "Rivet name and version string") AC_DEFINE_UNQUOTED(RIVET_TARNAME, "$PACKAGE_TARNAME", "Rivet short name string") AC_DEFINE_UNQUOTED(RIVET_BUGREPORT, "$PACKAGE_BUGREPORT", "Rivet contact email address") ## OS X AC_CEDAR_OSX ## Work out the LCG platform tag AC_LCG_TAG ## Set default compiler flags if test "x$CXXFLAGS" == "x"; then CXXFLAGS="-O2"; fi ## Compiler setup AC_LANG(C++) AC_PROG_CXX AX_CXX_COMPILE_STDCXX([11], [noext], [mandatory]) ## Store and propagate the compiler identity and flags RIVETCXX="$CXX" AC_SUBST(RIVETCXX) RIVETCXXFLAGS="$CXXFLAGS" AC_SUBST(RIVETCXXFLAGS) ## Checks for programs. AC_PROG_INSTALL AC_PROG_LN_S AC_DISABLE_STATIC AC_LIBTOOL_DLOPEN AC_PROG_LIBTOOL AX_EXECINFO AC_FUNC_STRERROR_R ## YODA histogramming library AC_CEDAR_LIBRARYANDHEADERS([YODA], , , [AC_MSG_ERROR([YODA is required])]) YODABINPATH=$YODALIBPATH/../bin AC_SUBST(YODABINPATH) AC_PATH_PROG(YODACONFIG, yoda-config, [], [$YODALIBPATH/../bin:$PATH]) YODA_PYTHONPATH="" if test -f "$YODACONFIG"; then AC_MSG_CHECKING([YODA version using yoda-config]) YODA_VERSION=`$YODACONFIG --version` AC_MSG_RESULT([$YODA_VERSION]) YODA_VERSION1=[`echo $YODA_VERSION | cut -d. -f1 | sed -e 's/\([0-9]*\).*/\1/g'`] YODA_VERSION2=[`echo $YODA_VERSION | cut -d. -f2 | sed -e 's/\([0-9]*\).*/\1/g'`] YODA_VERSION3=[`echo $YODA_VERSION | cut -d. -f3 | sed -e 's/\([0-9]*\).*/\1/g'`] let YODA_VERSION_INT=YODA_VERSION1*10000+YODA_VERSION2*100+YODA_VERSION3 if test $YODA_VERSION_INT -lt 10500; then AC_MSG_ERROR([YODA version isn't sufficient: at least version 1.5.0 required]) fi AC_MSG_CHECKING([YODA Python path using yoda-config]) YODA_PYTHONPATH=`$YODACONFIG --pythonpath` AC_MSG_RESULT([$YODA_PYTHONPATH]) fi AC_SUBST(YODA_PYTHONPATH) ## HepMC event record library AC_CEDAR_LIBRARYANDHEADERS([HepMC], , , [AC_MSG_ERROR([HepMC is required])]) oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS -I$HEPMCINCPATH" if test -e "$HEPMCINCPATH/HepMC/HepMCDefs.h"; then AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #include "HepMC/HepMCDefs.h" int main() { std::cout << HEPMC_VERSION << std::endl; return 0; }])]) else AC_LANG_CONFTEST([AC_LANG_SOURCE([#include #include "HepMC/defs.h" int main() { std::cout << VERSION << std::endl; return 0; }])]) fi if test -f conftest.cc; then $CXX $CPPFLAGS conftest.cc -o conftest 2>&1 1>&5 else $CXX $CPPFLAGS conftest.cpp -o conftest 2>&1 1>&5 fi hepmc_version=`./conftest` AC_MSG_CHECKING([HepMC version]) AX_COMPARE_VERSION([$hepmc_version], [lt], [2.06.09], AC_MSG_RESULT([no]) AC_MSG_ERROR([Need HepMC 2.06.09 or later]), AC_MSG_RESULT([ok])) if test x$hepmc_version != x; then let hepmc_major=`echo "$hepmc_version" | cut -d. -f1` let hepmc_minor=`echo "$hepmc_version" | cut -d. -f2` fi rm -f conftest conftest.cpp conftest.cc conftest.C HEPMC_VERSION=$hepmc_major$hepmc_minor AC_MSG_NOTICE([HepMC version is $hepmc_version -> $HEPMC_VERSION]) AC_SUBST(HEPMC_VERSION) CPPFLAGS=$oldCPPFLAGS ## FastJet clustering library AC_CEDAR_LIBRARYANDHEADERS([fastjet], , , [AC_MSG_ERROR([FastJet is required])]) AC_PATH_PROG(FJCONFIG, fastjet-config, [], $FASTJETPATH/bin:$PATH) if test -f "$FJCONFIG"; then AC_MSG_CHECKING([FastJet version using fastjet-config]) fjversion=`$FJCONFIG --version` AC_MSG_RESULT([$fjversion]) fjmajor=$(echo $fjversion | cut -f1 -d.) fjminor=$(echo $fjversion | cut -f2 -d.) fjmicro=$(echo $fjversion | cut -f3 -d.) if test "$fjmajor" -lt 3; then AC_MSG_ERROR([FastJet version 3.0.0 or later is required]) fi FASTJETCONFIGLIBADD="$($FJCONFIG --plugins --shared --libs)" else FASTJETCONFIGLIBADD="-L$FASTJETLIBPATH -l$FASTJETLIBNAME" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lSISConePlugin -lsiscone -lsiscone_spherical" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lCDFConesPlugin -lD0RunIIConePlugin -lNestedDefsPlugin" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lTrackJetPlugin -lATLASConePlugin -lCMSIterativeConePlugin" FASTJETCONFIGLIBADD="$FASTJETCONFIGLIBADD -lEECambridgePlugin -lJadePlugin" fi; AC_MSG_NOTICE([FastJet LIBADD = $FASTJETCONFIGLIBADD]) AC_SUBST(FASTJETCONFIGLIBADD) # Check for FastJet headers that require the --enable-all(cxx)plugins option FASTJET_ERRMSG="Required FastJet plugin headers were not found: did you build FastJet with the --enable-allcxxplugins option?" oldCPPFLAGS=$CPPFLAGS CPPFLAGS="$CPPFLAGS -I$FASTJETINCPATH" AC_CHECK_HEADER([fastjet/D0RunIIConePlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) AC_CHECK_HEADER([fastjet/TrackJetPlugin.hh], [], [AC_MSG_ERROR([$FASTJET_ERRMSG])]) CPPFLAGS=$oldCPPFLAGS ## Disable build/install of standard analyses AC_ARG_ENABLE([analyses], [AC_HELP_STRING(--disable-analyses, [don't try to build or install standard analyses])], [], [enable_analyses=yes]) if test x$enable_analyses != xyes; then AC_MSG_WARN([Not building standard Rivet analyses, by request]) fi AM_CONDITIONAL(ENABLE_ANALYSES, [test x$enable_analyses = xyes]) ## Build LaTeX docs if possible... AC_PATH_PROG(PDFLATEX, pdflatex) AM_CONDITIONAL(WITH_PDFLATEX, [test x$PDFLATEX != x]) ## ... unless told otherwise! AC_ARG_ENABLE([pdfmanual], [AC_HELP_STRING(--enable-pdfmanual, [build and install the manual])], [], [enable_pdfmanual=no]) if test x$enable_pdfmanual = xyes; then AC_MSG_WARN([Building Rivet manual, by request]) fi AM_CONDITIONAL(ENABLE_PDFMANUAL, [test x$enable_pdfmanual = xyes]) ## Build Doxygen documentation if possible AC_ARG_ENABLE([doxygen], [AC_HELP_STRING(--disable-doxygen, [don't try to make Doxygen documentation])], [], [enable_doxygen=yes]) if test x$enable_doxygen = xyes; then AC_PATH_PROG(DOXYGEN, doxygen) fi AM_CONDITIONAL(WITH_DOXYGEN, [test x$DOXYGEN != x]) ## Build asciidoc docs if possible AC_PATH_PROG(ASCIIDOC, asciidoc) AM_CONDITIONAL(WITH_ASCIIDOC, [test x$ASCIIDOC != x]) ## Python extension AC_ARG_ENABLE(pyext, [AC_HELP_STRING(--disable-pyext, [don't build Python module (default=build)])], [], [enable_pyext=yes]) ## Basic Python checks if test x$enable_pyext == xyes; then AX_PYTHON_DEVEL([>= '2.7.3']) AC_SUBST(PYTHON_VERSION) RIVET_PYTHONPATH=`$PYTHON -c "from __future__ import print_function; import distutils.sysconfig; print(distutils.sysconfig.get_python_lib(prefix='$prefix', plat_specific=True));"` AC_SUBST(RIVET_PYTHONPATH) if test -z "$PYTHON"; then AC_MSG_ERROR([Can't build Python extension since python can't be found]) enable_pyext=no fi if test -z "$PYTHON_CPPFLAGS"; then AC_MSG_ERROR([Can't build Python extension since Python.h header file cannot be found]) enable_pyext=no fi fi AM_CONDITIONAL(ENABLE_PYEXT, [test x$enable_pyext == xyes]) dnl dnl setup.py puts its build artifacts into a labelled path dnl this helps the test scripts to find them locally instead of dnl having to install first dnl RIVET_SETUP_PY_PATH=$(${PYTHON} -c 'from __future__ import print_function; import distutils.util as u, sys; vi=sys.version_info; print("lib.%s-%s.%s" % (u.get_platform(),vi.major, vi.minor))') AC_SUBST(RIVET_SETUP_PY_PATH) ## Cython checks if test x$enable_pyext == xyes; then AM_CHECK_CYTHON([0.24.0], [:], [:]) if test x$CYTHON_FOUND = xyes; then AC_MSG_NOTICE([Cython >= 0.24 found: Python extension source can be rebuilt (for developers)]) fi AC_CHECK_FILE([pyext/rivet/core.cpp], [], [if test "x$CYTHON_FOUND" != "xyes"; then AC_MSG_ERROR([Cython is required for --enable-pyext, no pre-built core.cpp was found.]) fi]) cython_compiler=$CXX ## Set extra Python extension build flags (to cope with Cython output code oddities) PYEXT_CXXFLAGS="$CXXFLAGS" AC_CEDAR_CHECKCXXFLAG([-Wno-unused-but-set-variable], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-unused-but-set-variable"]) AC_CEDAR_CHECKCXXFLAG([-Wno-sign-compare], [PYEXT_CXXFLAGS="$PYEXT_CXXFLAGS -Wno-sign-compare"]) AC_SUBST(PYEXT_CXXFLAGS) AC_MSG_NOTICE([All Python build checks successful: 'rivet' Python extension will be built]) fi AM_CONDITIONAL(WITH_CYTHON, [test x$CYTHON_FOUND = xyes]) ## Set default build flags AM_CPPFLAGS="-I\$(top_srcdir)/include -I\$(top_builddir)/include" -#AM_CPPFLAGS="$AM_CPPFLAGS -I\$(top_srcdir)/include/eigen3" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(YODAINCPATH)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(HEPMCINCPATH)" AM_CPPFLAGS="$AM_CPPFLAGS -I\$(FASTJETINCPATH)" AC_CEDAR_CHECKCXXFLAG([-pedantic], [AM_CXXFLAGS="$AM_CXXFLAGS -pedantic"]) AC_CEDAR_CHECKCXXFLAG([-Wall], [AM_CXXFLAGS="$AM_CXXFLAGS -Wall"]) AC_CEDAR_CHECKCXXFLAG([-Wno-long-long], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-long-long"]) AC_CEDAR_CHECKCXXFLAG([-Wno-format], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-format"]) dnl AC_CEDAR_CHECKCXXFLAG([-Wno-unused-variable], [AM_CXXFLAGS="$AM_CXXFLAGS -Wno-unused-variable"]) AC_CEDAR_CHECKCXXFLAG([-Werror=uninitialized], [AM_CXXFLAGS="$AM_CXXFLAGS -Werror=uninitialized"]) AC_CEDAR_CHECKCXXFLAG([-Werror=delete-non-virtual-dtor], [AM_CXXFLAGS="$AM_CXXFLAGS -Werror=delete-non-virtual-dtor"]) ## Add OpenMP-enabling flags if possible AX_OPENMP([AM_CXXFLAGS="$AM_CXXFLAGS $OPENMP_CXXFLAGS"]) ## Optional zlib support for gzip-compressed data streams/files AX_CHECK_ZLIB ## Debug flag (default=-DNDEBUG, enabled=-g) AC_ARG_ENABLE([debug], [AC_HELP_STRING(--enable-debug, [build with debugging symbols @<:@default=no@:>@])], [], [enable_debug=no]) if test x$enable_debug == xyes; then AM_CXXFLAGS="$AM_CXXFLAGS -g" fi ## Extra warnings flag (default=none) AC_ARG_ENABLE([extra-warnings], [AC_HELP_STRING(--enable-extra-warnings, [build with extra compiler warnings (recommended for developers) @<:@default=no@:>@])], [], [enable_extra_warnings=no]) if test x$enable_extra_warnings == xyes; then AC_CEDAR_CHECKCXXFLAG([-Wextra], [AM_CXXFLAGS="$AM_CXXFLAGS -Wextra "]) fi AC_SUBST(AM_CPPFLAGS) AC_SUBST(AM_CXXFLAGS) AC_EMPTY_SUBST AC_CONFIG_FILES(Makefile Doxyfile) AC_CONFIG_FILES(include/Makefile include/Rivet/Makefile) AC_CONFIG_FILES(src/Makefile) AC_CONFIG_FILES(src/Core/Makefile src/Core/yamlcpp/Makefile) AC_CONFIG_FILES(src/Tools/Makefile) AC_CONFIG_FILES(src/Tools/Nsubjettiness/Makefile) AC_CONFIG_FILES(src/Projections/Makefile) AC_CONFIG_FILES(src/AnalysisTools/Makefile) AC_CONFIG_FILES(analyses/Makefile) AC_CONFIG_FILES(test/Makefile) AC_CONFIG_FILES(pyext/Makefile pyext/rivet/Makefile pyext/setup.py) AC_CONFIG_FILES(data/Makefile data/texmf/Makefile) AC_CONFIG_FILES(doc/Makefile) AC_CONFIG_FILES(doc/rivetversion.sty) AC_CONFIG_FILES(bin/Makefile bin/rivet-config bin/rivet-buildplugin) AC_CONFIG_FILES(rivetenv.sh rivetenv.csh rivet.pc) AC_OUTPUT if test x$enable_pyrivet == xyes; then cat <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); } - /// This Particle's PDG ID code (alias). - /// @deprecated Prefer the pid/abspid form - PdgId pdgId() const { return _id; } //@} /// @name Charge //@{ /// The charge of this Particle. double charge() const { return PID::charge(pid()); } /// The absolute charge of this Particle. double abscharge() const { return PID::abscharge(pid()); } /// Three times the charge of this Particle (i.e. integer multiple of smallest quark charge). int charge3() const { return PID::charge3(pid()); } /// 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 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/src/Projections/ChargedFinalState.cc b/src/Projections/ChargedFinalState.cc --- a/src/Projections/ChargedFinalState.cc +++ b/src/Projections/ChargedFinalState.cc @@ -1,48 +1,48 @@ // -*- C++ -*- #include "Rivet/Projections/ChargedFinalState.hh" namespace Rivet { ChargedFinalState::ChargedFinalState(const FinalState& fsp) { setName("ChargedFinalState"); declare(fsp, "FS"); } ChargedFinalState::ChargedFinalState(const Cut& c) { setName("ChargedFinalState"); declare(FinalState(c), "FS"); } ChargedFinalState::ChargedFinalState(double mineta, double maxeta, double minpt) { setName("ChargedFinalState"); declare(FinalState(mineta, maxeta, minpt), "FS"); } CmpState ChargedFinalState::compare(const Projection& p) const { return mkNamedPCmp(p, "FS"); } } namespace { inline bool chargedParticleFilter(const Rivet::Particle& p) { - return Rivet::PID::charge3(p.pdgId()) == 0; + return Rivet::PID::charge3(p.pid()) == 0; } } namespace Rivet { void ChargedFinalState::project(const Event& e) { const FinalState& fs = applyProjection(e, "FS"); _theParticles.clear(); std::remove_copy_if(fs.particles().begin(), fs.particles().end(), std::back_inserter(_theParticles), chargedParticleFilter); MSG_DEBUG("Number of charged final-state particles = " << _theParticles.size()); if (getLog().isActive(Log::TRACE)) { for (vector::iterator p = _theParticles.begin(); p != _theParticles.end(); ++p) { - MSG_TRACE("Selected: " << p->pdgId() << ", charge = " << PID::charge3(p->pdgId())/3.0); + MSG_TRACE("Selected: " << p->pid() << ", charge = " << PID::charge3(p->pid())/3.0); } } } }