diff --git a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_037.cc b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_037.cc --- a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_037.cc +++ b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_037.cc @@ -1,204 +1,204 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedMET.hh" #include "Rivet/Tools/Cutflow.hh" namespace Rivet { /// @brief ATLAS 2016 2 -SS-lepton / 3-lepton SUSY search, from 13.2/fb CONF note class ATLAS_2016_CONF_2016_037 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_CONF_2016_037); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections FinalState calofs(Cuts::abseta < 4.9); declare(calofs, "Clusters"); FastJets fj(calofs, FastJets::ANTIKT, 0.4); declare(fj, "TruthJets"); declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, [](const Jet& j) { if (j.abseta() > 2.5) return 0.; return j.bTagged(Cuts::pT > 5*GeV) ? 0.70 : j.cTagged(Cuts::pT > 5*GeV) ? 1/12. : j.tauTagged(Cuts::pT > 5*GeV) ? 1/54. : 1/380.; }), "Jets"); MissingMomentum mm(calofs); declare(mm, "TruthMET"); declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "MET"); FinalState es(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && !Cuts::absetaIn(1.37, 1.52) && Cuts::pT > 10*GeV); declare(es, "TruthElectrons"); - declare(SmearedParticles(es, ELECTRON_EFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons"); + declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons"); FinalState mus(Cuts::abspid == PID::MUON && Cuts::abseta < 2.5 && Cuts::pT > 10*GeV); declare(mus, "TruthMuons"); declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "Muons"); FinalState cfs(Cuts::abseta < 2.5 && Cuts::abscharge > 0); declare(cfs, "TruthTracks"); declare(SmearedParticles(cfs, TRK_EFF_ATLAS_RUN2), "Tracks"); // Book histograms/counters book(_h_3l1,"SR3l1"); book(_h_3l2,"SR3l2"); book(_h_0b1,"SR0b1"); book(_h_0b2,"SR0b2"); book(_h_1b,"SR1b"); book(_h_3b,"SR3b"); book(_h_1bDD,"SR1bDD"); book(_h_3bDD,"SR3bDD"); book(_h_1bGG,"SR1bGG"); } /// Perform the per-event analysis void analyze(const Event& event) { // Get baseline electrons, muons, and jets Particles elecs = apply(event, "Electrons").particlesByPt(); Particles muons = apply(event, "Muons").particlesByPt(); Jets jets = apply(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8); const Jets bjets = filter_select(jets, [&](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); }); // Jet/electron/muon overlap removal and selection // Remove any electron or muon within dR = 0.2 of a b-tagged jet for (const Jet& bj : bjets) { ifilter_discard(elecs, deltaRLess(bj, 0.2, RAPIDITY)); ifilter_discard(muons, deltaRLess(bj, 0.2, RAPIDITY)); } // Remove any untagged jet within dR = 0.2 of an electron or muon for (const Particle& e : elecs) ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY)); for (const Particle& m : muons) ifilter_discard(jets, deltaRLess(m, 0.2, RAPIDITY)); // Remove any untagged low-multiplicity/muon-dominated jet within dR = 0.4 of a muon for (const Particle& m : muons) ifilter_discard(jets, [&](const Jet& j) { if (deltaR(m, j, RAPIDITY) > 0.4) return false; const Particles trks = j.particles(Cuts::abscharge != 0); if (trks.size() < 3) return true; return m.pT()/j.pT() > 0.5 && m.pT()/sum(trks, pT, 0.0) > 0.7; }); // Remove any electron or muon near a remaining jet, with a shrinking cone const auto lcone_iso_fn = [&](const Particle& l) { const double dr = min(0.4, 0.04 + 10*GeV/l.pT()); return any(jets, deltaRLess(l, dr, RAPIDITY)); }; ifilter_discard(elecs, lcone_iso_fn); ifilter_discard(muons, lcone_iso_fn); // Track-sharing e,mu also filtered, but that decision can't be made here const Jets& sigjets = jets; const Jets& sigbjets = bjets; // Lepton isolation Particles sigelecs = filter_select(elecs, Cuts::abseta < 2); Particles sigmuons = muons; - ifilter_select(sigelecs, ParticleEffFilter(ELECTRON_IDEFF_ATLAS_RUN2_MEDIUM)); + ifilter_select(sigelecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_MEDIUM)); const Particles trks = apply(event, "Tracks").particles(); const Particles clus = apply(event, "Clusters").particles(); ifilter_discard(sigelecs, [&](const Particle& e) { const double R = min(0.2, 10*GeV/e.pT()); double ptsum = -e.pT(), etsum = -e.Et(); for (const Particle& t : trks) if (deltaR(t,e) < R) ptsum += t.pT(); for (const Particle& c : clus) if (deltaR(c,e) < 0.2) etsum += c.pT(); ///< @todo Bit vague about "energy" return ptsum / e.pT() > 0.06 || etsum / e.pT() > 0.06; }); ifilter_discard(sigmuons, [&](const Particle& m) { const double R = min(0.3, 10*GeV/m.pT()); double ptsum = -m.pT(); for (const Particle& t : trks) if (deltaR(t,m) < R) ptsum += t.pT(); return ptsum / m.pT() > 0.06; }); /// @todo Note is vague about whether "signal lepton" defn includes pT > 20? ifilter_discard(sigelecs, Cuts::pT > 20*GeV); ifilter_discard(sigmuons, Cuts::pT > 20*GeV); // MET calculation (NB. done generically, with smearing, rather than via explicit physics objects) const Vector3 vmet = -apply(event, "MET").vectorEt(); const double etmiss = vmet.mod(); ////////////////// // Event selection cuts const Particles sigleptons = sigelecs + sigmuons; if (sigleptons.size() < 2) vetoEvent; if (sigleptons.size() == 2 && sigleptons[0].charge() != sigleptons[1].charge()) vetoEvent; // Jet sub-selections and meff calculation const Jets sigjets25 = filter_select(sigjets, Cuts::pT > 25*GeV); const Jets sigjets40 = filter_select(sigjets25, Cuts::pT > 40*GeV); const Jets sigjets50 = filter_select(sigjets40, Cuts::pT > 50*GeV); /// @todo Is meff specific to the jet pT cut? const double meff = sum(sigjets, pT, 0.0) + sum(sigleptons, pT, 0.0); // Fill counters if (sigleptons.size() >= 3 && sigbjets.empty() && sigjets40.size() >= 4 && etmiss > 150*GeV) _h_3l1->fill(); if (sigleptons.size() >= 3 && sigbjets.empty() && sigjets40.size() >= 4 && etmiss > 200*GeV && meff > 1500*GeV) _h_3l2->fill(); if (sigleptons.size() >= 2 && sigbjets.empty() && sigjets25.size() >= 6 && etmiss > 150*GeV && meff > 500*GeV) _h_0b1->fill(); if (sigleptons.size() >= 2 && sigbjets.empty() && sigjets40.size() >= 6 && etmiss > 150*GeV && meff > 900*GeV) _h_0b2->fill(); if (sigleptons.size() >= 2 && sigbjets.size() >= 1 && sigjets25.size() >= 6 && etmiss > 200*GeV && meff > 650*GeV) _h_1b->fill(); if (sigleptons.size() >= 2 && sigbjets.size() >= 3 && sigjets25.size() >= 6 && etmiss > 150*GeV && meff > 600*GeV) _h_3b->fill(); if (filter_select(sigleptons, Cuts::charge < 0).size() >= 2) { if (sigleptons.size() >= 2 && sigbjets.size() >= 1 && sigjets50.size() >= 6 && meff > 1200*GeV) _h_1bDD->fill(); if (sigleptons.size() >= 2 && sigbjets.size() >= 3 && sigjets50.size() >= 6 && meff > 1000*GeV) _h_3bDD->fill(); if (sigleptons.size() >= 2 && sigbjets.size() >= 1 && sigjets50.size() >= 6 && meff > 1800*GeV) _h_1bGG->fill(); } } /// Normalise counters after the run void finalize() { const double sf = 13.2*crossSection()/femtobarn/sumOfWeights(); scale(_h_3l1, sf); scale(_h_3l2, sf); scale(_h_0b1, sf); scale(_h_0b2, sf); scale(_h_1b, sf); scale(_h_3b, sf); scale(_h_1bDD, sf); scale(_h_3bDD, sf); scale(_h_1bGG, sf); } //@} private: /// @name Histograms //@{ CounterPtr _h_3l1, _h_3l2, _h_0b1, _h_0b2, _h_1b, _h_3b, _h_1bDD, _h_3bDD, _h_1bGG; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_CONF_2016_037); } diff --git a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_054.cc b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_054.cc --- a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_054.cc +++ b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_054.cc @@ -1,218 +1,218 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedMET.hh" #include "Rivet/Tools/Cutflow.hh" namespace Rivet { /// @brief ATLAS 2016 1-lepton SUSY search, from 14.8/fb CONF note class ATLAS_2016_CONF_2016_054 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_CONF_2016_054); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections FinalState calofs(Cuts::abseta < 4.9); FastJets fj(calofs, FastJets::ANTIKT, 0.4); declare(fj, "TruthJets"); declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, [](const Jet& j) { if (j.abseta() > 2.5) return 0.; return j.bTagged(Cuts::pT > 5*GeV) ? 0.77 : j.cTagged(Cuts::pT > 5*GeV) ? 1/6.2 : 1/134.; }), "Jets"); MissingMomentum mm(calofs); declare(mm, "TruthMET"); declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "MET"); FinalState es(Cuts::abseta < 2.47 && Cuts::pT > 7*GeV && Cuts::abspid == PID::ELECTRON); declare(es, "TruthElectrons"); - declare(SmearedParticles(es, ELECTRON_EFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons"); + declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons"); FinalState mus(Cuts::abseta < 2.5 && Cuts::pT > 6*GeV && Cuts::abspid == PID::MUON); declare(mus, "TruthMuons"); declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "Muons"); // Book histograms/counters book(_h_gg2j,"GG-2j"); book(_h_gg6j0,"GG-6j-0bulk"); book(_h_gg6j1,"GG-6j-1highmass"); book(_h_gg4j0,"GG-4j-0lowx"); book(_h_gg4j1,"GG-4j-1lowxbveto"); book(_h_gg4j2,"GG-4j-2highx"); book(_h_ss4j0,"SS-4j-0x12"); book(_h_ss4j1,"SS-4j-1lowx"); book(_h_ss5j0,"SS-5j-0x12"); book(_h_ss5j1,"SS-5j-1highx"); } /// Perform the per-event analysis void analyze(const Event& event) { // Get baseline electrons, muons, and jets Particles elecs = apply(event, "Electrons").particles(); Particles muons = apply(event, "Muons").particles(); Jets jets = apply(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 4.5); // Jet/electron/muons overlap removal and selection // Remove any jet within dR = 0.2 of an electron for (const Particle& e : elecs) ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY)); // Remove any electron within dR = 0.01 of a muon for (const Particle& m : muons) ifilter_discard(elecs, deltaRLess(m, 0.01, RAPIDITY)); // Assemble b-jets collection, and remove muons within dR = 0.2 of a b-tagged jet Jets bjets; for (const Jet& j : jets) { if (j.abseta() < 2.5 && j.pT() > 30*GeV && j.bTagged(Cuts::pT > 5*GeV)) { bjets += j; ifilter_discard(muons, deltaRLess(j, 0.2, RAPIDITY)); } } // Remove any jet within dR = 0.2 of a muon if track conditions are met for (const Particle& m : muons) ifilter_discard(jets, [&](const Jet& j){ if (deltaR(j,m) > 0.2) return false; /// @todo Add track efficiency random filtering const Particles trks = j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.5*GeV); return trks.size() < 4 || m.pT()/j.pT() > 0.7; }); // Remove any muon within dR = 0.2 of a remaining jet if the same track conditions are *not* met /// @todo There must be nicer way to do complementary removal... for (const Jet& j : jets) { /// @todo Add track efficiency random filtering const size_t ntrks = j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.5*GeV).size(); ifilter_discard(muons, [&](const Particle& m){ if (deltaR(j,m) > 0.2) return false; return ntrks > 3 && m.pT()/j.pT() < 0.7; }); } // Remove any muon with dR close to a remaining jet, via a functional form for (const Jet& j : jets) ifilter_discard(muons, [&](const Particle& m) { return deltaR(m,j, RAPIDITY) < min(0.4, 0.04 + 10*GeV/m.pT()); }); // Signal jet selection const Jets sigjets = filter_select(jets, Cuts::pT > 30*GeV && Cuts::abseta < 2.8); const Jets sigbjets = bjets; // "Gradient-loose" signal lepton selection const ParticleEffFilter grad_loose_filter([](const Particle& e) { return e.pT() > 60*GeV ? 0.98 : 0.95; }); Particles sigelecs = filter_select(elecs, grad_loose_filter); Particles sigmuons = filter_select(muons, grad_loose_filter); // Tight electron selection (NB. assuming independent eff to gradient-loose... hmm) - ifilter_select(sigelecs, ParticleEffFilter(ELECTRON_IDEFF_ATLAS_RUN2_TIGHT)); + ifilter_select(sigelecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_TIGHT)); // MET calculation (NB. done generically, with smearing, rather than via explicit physics objects) const Vector3 vmet = -apply(event, "MET").vectorEt(); const double etmiss = vmet.mod(); ////////////////// // Event selection cuts if (sigelecs.size() + sigmuons.size() != 1) vetoEvent; const Particle siglepton = sigelecs.empty() ? sigmuons.front() : sigelecs.front(); // mT and m_eff const double mT = sqrt(2*siglepton.pT()*etmiss*(1-cos(deltaPhi(siglepton,vmet)))); const double meff = siglepton.pT() + sum(sigjets, pT, 0.0) + etmiss; // Aplanarities Sphericity sph; vector vecs; transform(sigjets, vecs, mom); sph.calc(vecs); const double jet_aplanarity = sph.aplanarity(); vecs += siglepton.mom(); sph.calc(vecs); const double lepton_aplanarity = sph.aplanarity(); ////////////////// // Fill counters // GG if (siglepton.pT() < 35*GeV && sigjets.size() >= 2 && sigjets[0].pT() > 200*GeV && sigjets[1].pT() > 30*GeV && mT > 100*GeV && etmiss > 460*GeV && etmiss/meff > 0.35) _h_gg2j->fill(); if (siglepton.pT() > 35*GeV && sigjets.size() >= 6 && sigjets[0].pT() > 125*GeV && sigjets[5].pT() > 30*GeV && mT > 225*GeV && etmiss > 250*GeV && meff > 1000*GeV && etmiss/meff > 0.2 && jet_aplanarity > 0.04) _h_gg6j0->fill(); if (siglepton.pT() > 35*GeV && sigjets.size() >= 6 && sigjets[0].pT() > 125*GeV && sigjets[5].pT() > 30*GeV && mT > 225*GeV && etmiss > 250*GeV && meff > 2000*GeV && etmiss/meff > 0.1 && jet_aplanarity > 0.04) _h_gg6j1->fill(); if (sigjets.size() >= 4 && sigjets[3].pT() > 100*GeV && mT > 125*GeV && etmiss > 250*GeV && meff > 2000*GeV && jet_aplanarity > 0.06) _h_gg4j0->fill(); if (sigjets.size() >= 4 && sigjets[3].pT() > 100*GeV && sigbjets.empty() && mT > 125*GeV && etmiss > 250*GeV && meff > 2000*GeV && jet_aplanarity > 0.03) _h_gg4j1->fill(); if (siglepton.pT() > 35*GeV && sigjets.size() >= 4 && sigjets[0].pT() > 400*GeV && inRange(sigjets[3].pT(), 30*GeV, 100*GeV) && mT > 475*GeV && etmiss > 250*GeV && meff > 1600*GeV && etmiss/meff > 0.3) _h_gg4j2->fill(); // SS if (siglepton.pT() > 35*GeV && sigjets.size() >= 4 && sigjets[3].pT() > 50*GeV && mT > 175*GeV && etmiss > 300*GeV && meff > 1200*GeV && lepton_aplanarity > 0.08) _h_ss4j0->fill(); if (siglepton.pT() > 35*GeV && sigjets.size() >= 5 && sigjets[4].pT() > 50*GeV && sigbjets.empty() && mT > 175*GeV && etmiss > 300*GeV && etmiss/meff > 0.2) _h_ss5j0->fill(); if (siglepton.pT() > 35*GeV && sigjets.size() >= 4 && sigjets[0].pT() > 250*GeV && sigjets[3].pT() > 30*GeV && inRange(mT, 150*GeV, 400*GeV) && etmiss > 250*GeV && lepton_aplanarity > 0.03) _h_ss4j1->fill(); if (siglepton.pT() > 35*GeV && sigjets.size() >= 5 && sigjets[4].pT() > 30*GeV && mT > 400*GeV && etmiss > 400*GeV && lepton_aplanarity > 0.03) _h_ss5j1->fill(); } /// Normalise counters after the run void finalize() { const double sf = 14.8*crossSection()/femtobarn/sumOfWeights(); scale(_h_gg2j, sf); scale(_h_gg6j0, sf); scale(_h_gg6j1, sf); scale(_h_gg4j0, sf); scale(_h_gg4j1, sf); scale(_h_gg4j2, sf); scale(_h_ss4j0, sf); scale(_h_ss4j1, sf); scale(_h_ss5j0, sf); scale(_h_ss5j1, sf); } //@} private: /// @name Histograms //@{ CounterPtr _h_gg2j, _h_gg6j0, _h_gg6j1, _h_gg4j0, _h_gg4j1, _h_gg4j2; CounterPtr _h_ss4j0, _h_ss4j1, _h_ss5j0,_h_ss5j1; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_CONF_2016_054); } diff --git a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_078.cc b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_078.cc --- a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_078.cc +++ b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_078.cc @@ -1,267 +1,267 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedMET.hh" #include "Rivet/Tools/Cutflow.hh" namespace Rivet { /// @brief ATLAS 2016 0-lepton SUSY search, from 13/fb ICHEP'16 CONF note class ATLAS_2016_CONF_2016_078 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_CONF_2016_078); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections FinalState calofs(Cuts::abseta < 3.2); FastJets fj(calofs, FastJets::ANTIKT, 0.4); declare(fj, "TruthJets"); declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, //JET_BTAG_ATLAS_RUN2_MV2C10 [](const Jet& j) { if (j.abseta() > 2.5) return 0.; return j.bTagged(Cuts::pT > 5*GeV) ? 0.77 : j.cTagged(Cuts::pT > 5*GeV) ? 1/6. : 1/134.; }), "RecoJets"); MissingMomentum mm(calofs); declare(mm, "TruthMET"); declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "RecoMET"); PromptFinalState es(Cuts::abseta < 2.47 && Cuts::abspid == PID::ELECTRON, true, true); declare(es, "TruthElectrons"); - declare(SmearedParticles(es, ELECTRON_EFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "RecoElectrons"); + declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "RecoElectrons"); PromptFinalState mus(Cuts::abseta < 2.7 && Cuts::abspid == PID::MUON, true); declare(mus, "TruthMuons"); declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "RecoMuons"); // Book histograms/counters book(_h_2j_0800,"2j-0800"); book(_h_2j_1200,"2j-1200"); book(_h_2j_1600,"2j-1600"); book(_h_2j_2000,"2j-2000"); book(_h_3j_1200,"2j-2000"); book(_h_4j_1000,"4j-1000"); book(_h_4j_1400,"4j-1400"); book(_h_4j_1800,"4j-1800"); book(_h_4j_2200,"4j-2200"); book(_h_4j_2600,"4j-2600"); book(_h_5j_1400,"5j-1400"); book(_h_6j_1800,"6j-1800"); book(_h_6j_2200,"6j-2200"); // Book cut-flows const vector cuts23j = {"Pre-sel+MET+pT1+meff", "Njet", "Dphi_min(j123,MET)", "Dphi_min(j4+,MET)", "pT2", "eta_j12", "MET/sqrtHT", "m_eff(incl)"}; _flows.addCutflow("2j-0800", cuts23j); _flows.addCutflow("2j-1200", cuts23j); _flows.addCutflow("2j-1600", cuts23j); _flows.addCutflow("2j-2000", cuts23j); _flows.addCutflow("3j-1200", cuts23j); const vector cuts456j = {"Pre-sel+MET+pT1+meff", "Njet", "Dphi_min(j123,MET)", "Dphi_min(j4+,MET)", "pT4", "eta_j1234", "Aplanarity", "MET/m_eff(Nj)", "m_eff(incl)"}; _flows.addCutflow("4j-1000", cuts456j); _flows.addCutflow("4j-1400", cuts456j); _flows.addCutflow("4j-1800", cuts456j); _flows.addCutflow("4j-2200", cuts456j); _flows.addCutflow("4j-2600", cuts456j); _flows.addCutflow("5j-1400", cuts456j); _flows.addCutflow("6j-1800", cuts456j); _flows.addCutflow("6j-2200", cuts456j); } /// Perform the per-event analysis void analyze(const Event& event) { _flows.fillinit(); // Same MET cut for all signal regions const Vector3 vmet = -apply(event, "RecoMET").vectorEt(); const double met = vmet.mod(); if (met < 250*GeV) vetoEvent; // Get baseline electrons, muons, and jets Particles elecs = apply(event, "RecoElectrons").particles(Cuts::pT > 10*GeV); Particles muons = apply(event, "RecoMuons").particles(Cuts::pT > 10*GeV); Jets jets = apply(event, "RecoJets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8); ///< @todo Pile-up subtraction // Jet/electron/muons overlap removal and selection // Remove electrons within dR = 0.2 of a b-tagged jet for (const Jet& j : jets) if (j.abseta() < 2.5 && j.pT() > 50*GeV && j.bTagged(Cuts::pT > 5*GeV)) ifilter_discard(elecs, deltaRLess(j, 0.2, RAPIDITY)); // Remove any |eta| < 2.8 jet within dR = 0.2 of a remaining electron for (const Particle& e : elecs) ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY)); // Remove any electron with dR in [0.2, 0.4] of a remaining jet for (const Jet& j : jets) ifilter_discard(elecs, [&](const Particle& e) { return inRange(deltaR(e,j, RAPIDITY), 0.2, 0.4); }); // Remove any muon with dR close to a remaining jet, via a functional form for (const Jet& j : jets) ifilter_discard(muons, [&](const Particle& m) { return deltaR(m,j, RAPIDITY) < min(0.4, 0.04 + 10*GeV/m.pT()); }); // Remove any |eta| < 2.8 jet within dR = 0.2 of a remaining muon if track conditions are met for (const Particle& m : muons) /// @todo Add track efficiency random filtering ifilter_discard(jets, [&](const Jet& j) { if (deltaR(j,m, RAPIDITY) > 0.2) return false; const Particles trks = j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.5*GeV); return trks.size() < 3 || (m.pT() > 2*j.pT() && m.pT() > 0.7*sum(trks, pT, 0.0)); }); // Loose electron selection - ifilter_select(elecs, ParticleEffFilter(ELECTRON_IDEFF_ATLAS_RUN2_LOOSE)); + ifilter_select(elecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_LOOSE)); // Veto the event if there are any remaining baseline leptons if (!elecs.empty()) vetoEvent; if (!muons.empty()) vetoEvent; // Passed presel & MET _flows.fill(1); // Get jets and their pTs const Jets jets20 = jets; const Jets jets50 = filterBy(jets, Cuts::pT > 50*GeV); const size_t njets50 = jets50.size(), njets20 = jets20.size(); if (jets50.size() < 2) vetoEvent; vector jetpts20, jetpts50; transform(jets20, jetpts20, pT); transform(jets50, jetpts50, pT); // Construct multi-jet observables const double ht = sum(jetpts20, 0.0); const double met_sqrtHT = met / sqrt(ht); const double meff_incl = sum(jetpts50, met); const double meff_4 = (njets50 >= 4) ? sum(head(jetpts50, 4), met) : -1; const double meff_5 = (njets50 >= 5) ? sum(head(jetpts50, 5), met) : -1; const double meff_6 = (njets50 >= 6) ? sum(head(jetpts50, 6), met) : -1; const double met_meff_4 = met / meff_4; const double met_meff_5 = met / meff_5; const double met_meff_6 = met / meff_6; // Jet |eta|s vector jetetas20; transform(jets20, jetetas20, abseta); const double etamax_2 = (njets20 >= 2) ? max(head(jetetas20, 2)) : -1; const double etamax_4 = (njets20 >= 4) ? max(head(jetetas20, 4)) : -1; const double etamax_6 = (njets20 >= 6) ? max(head(jetetas20, 6)) : -1; // Get dphis between MET and jets vector dphimets50; transform(jets50, dphimets50, deltaPhiWRT(vmet)); const vector dphimets50_123 = head(dphimets50, 3); const vector dphimets50_more = tail(dphimets50, -3); const double dphimin_123 = !dphimets50_123.empty() ? min(dphimets50_123) : -1; const double dphimin_more = !dphimets50_more.empty() ? min(dphimets50_more) : -1; // Jet aplanarity Sphericity sph; sph.calc(jets50); const double aplanarity = sph.aplanarity(); ////////////////// // 2 jet regions if (dphimin_123 > 0.8 && dphimin_more > 0.4) { if (jetpts50[1] > 200*GeV && etamax_2 < 0.8) { //< implicit pT[0] cut if (met_sqrtHT > 14*sqrt(GeV) && meff_incl > 800*GeV) _h_2j_0800->fill(); } if (jetpts50[1] > 250*GeV && etamax_2 < 1.2) { //< implicit pT[0] cut if (met_sqrtHT > 16*sqrt(GeV) && meff_incl > 1200*GeV) _h_2j_1200->fill(); if (met_sqrtHT > 18*sqrt(GeV) && meff_incl > 1600*GeV) _h_2j_1600->fill(); if (met_sqrtHT > 20*sqrt(GeV) && meff_incl > 2000*GeV) _h_2j_2000->fill(); } } // 3 jet region if (njets50 >= 3 && dphimin_123 > 0.4 && dphimin_more > 0.2) { if (jetpts50[0] > 600*GeV && jetpts50[2] > 50*GeV) { //< implicit pT[1] cut if (met_sqrtHT > 16*sqrt(GeV) && meff_incl > 1200*GeV) _h_3j_1200->fill(); } } // 4 jet regions (note implicit pT[1,2] cuts) if (njets50 >= 4 && dphimin_123 > 0.4 && dphimin_more > 0.4 && jetpts50[0] > 200*GeV && aplanarity > 0.04) { if (jetpts50[3] > 100*GeV && etamax_4 < 1.2 && met_meff_4 > 0.25*sqrt(GeV) && meff_incl > 1000*GeV) _h_4j_1000->fill(); if (jetpts50[3] > 100*GeV && etamax_4 < 2.0 && met_meff_4 > 0.25*sqrt(GeV) && meff_incl > 1400*GeV) _h_4j_1400->fill(); if (jetpts50[3] > 100*GeV && etamax_4 < 2.0 && met_meff_4 > 0.20*sqrt(GeV) && meff_incl > 1800*GeV) _h_4j_1800->fill(); if (jetpts50[3] > 150*GeV && etamax_4 < 2.0 && met_meff_4 > 0.20*sqrt(GeV) && meff_incl > 2200*GeV) _h_4j_2200->fill(); if (jetpts50[3] > 150*GeV && met_meff_4 > 0.20*sqrt(GeV) && meff_incl > 2600*GeV) _h_4j_2600->fill(); } // 5 jet region (note implicit pT[1,2,3] cuts) if (njets50 >= 5 && dphimin_123 > 0.4 && dphimin_more > 0.2 && jetpts50[0] > 500*GeV) { if (jetpts50[4] > 50*GeV && met_meff_5 > 0.3*sqrt(GeV) && meff_incl > 1400*GeV) _h_5j_1400->fill(); } // 6 jet regions (note implicit pT[1,2,3,4] cuts) if (njets50 >= 6 && dphimin_123 > 0.4 && dphimin_more > 0.2 && jetpts50[0] > 200*GeV && aplanarity > 0.08) { if (jetpts50[5] > 50*GeV && etamax_6 < 2.0 && met_meff_6*sqrt(GeV) > 0.20 && meff_incl > 1800*GeV) _h_6j_1800->fill(); if (jetpts50[5] > 100*GeV && met_meff_6*sqrt(GeV) > 0.15 && meff_incl > 2200*GeV) _h_6j_2200->fill(); } // Cutflows _flows["2j-0800"].filltail({true, dphimin_123 > 0.8, dphimin_more > 0.4, jetpts50[1] > 200*GeV, etamax_2 < 0.8, met_sqrtHT > 14*sqrt(GeV), meff_incl > 800*GeV}); _flows["2j-1200"].filltail({true, dphimin_123 > 0.8, dphimin_more > 0.4, jetpts50[1] > 250*GeV, etamax_2 < 1.2, met_sqrtHT > 16*sqrt(GeV), meff_incl > 1200*GeV}); _flows["2j-1600"].filltail({true, dphimin_123 > 0.8, dphimin_more > 0.4, jetpts50[1] > 250*GeV, etamax_2 < 1.2, met_sqrtHT > 18*sqrt(GeV), meff_incl > 1600*GeV}); _flows["2j-2000"].filltail({true, dphimin_123 > 0.8, dphimin_more > 0.4, jetpts50[1] > 250*GeV, etamax_2 < 1.2, met_sqrtHT > 20*sqrt(GeV), meff_incl > 2000*GeV}); _flows["3j-1200"].filltail({njets50 >= 3, dphimin_123 > 0.4, dphimin_more > 0.2, jetpts50[0] > 600*GeV && jetpts50[2] > 50*GeV, true, met_sqrtHT > 16*sqrt(GeV), meff_incl > 1200*GeV}); _flows["4j-1000"].filltail({njets50 >= 4, dphimin_123 > 0.4, dphimin_more > 0.4, jetpts50[0] > 200*GeV && jetpts50[3] > 100*GeV, etamax_4 < 1.2, aplanarity > 0.04, met_meff_4 > 0.25*sqrt(GeV), meff_incl > 1000*GeV}); _flows["4j-1400"].filltail({njets50 >= 4, dphimin_123 > 0.4, dphimin_more > 0.4, jetpts50[0] > 200*GeV && jetpts50[3] > 100*GeV, etamax_4 < 2.0, aplanarity > 0.04, met_meff_4 > 0.25*sqrt(GeV), meff_incl > 1400*GeV}); _flows["4j-1800"].filltail({njets50 >= 4, dphimin_123 > 0.4, dphimin_more > 0.4, jetpts50[0] > 200*GeV && jetpts50[3] > 100*GeV, etamax_4 < 2.0, aplanarity > 0.04, met_meff_4 > 0.20*sqrt(GeV), meff_incl > 1800*GeV}); _flows["4j-2200"].filltail({njets50 >= 4, dphimin_123 > 0.4, dphimin_more > 0.4, jetpts50[0] > 200*GeV && jetpts50[3] > 150*GeV, etamax_4 < 2.0, aplanarity > 0.04, met_meff_4 > 0.20*sqrt(GeV), meff_incl > 2200*GeV}); _flows["4j-2600"].filltail({njets50 >= 4, dphimin_123 > 0.4, dphimin_more > 0.4, jetpts50[0] > 200*GeV && jetpts50[3] > 150*GeV, true, aplanarity > 0.04, met_meff_4 > 0.20*sqrt(GeV), meff_incl > 2600*GeV}); _flows["5j-1400"].filltail({njets50 >= 5, dphimin_123 > 0.4, dphimin_more > 0.2, jetpts50[0] > 500*GeV && jetpts50[4] > 50*GeV, true, true, met_meff_5 > 0.3*sqrt(GeV), meff_incl > 1400*GeV}); _flows["6j-1800"].filltail({njets50 >= 6, dphimin_123 > 0.4, dphimin_more > 0.2, jetpts50[0] > 200*GeV && jetpts50[5] > 50*GeV, etamax_6 < 2.0, aplanarity > 0.08, met_meff_6 > 0.20*sqrt(GeV), meff_incl > 1800*GeV}); _flows["6j-2200"].filltail({njets50 >= 6, dphimin_123 > 0.4, dphimin_more > 0.2, jetpts50[0] > 200*GeV && jetpts50[5] > 100*GeV, true, aplanarity > 0.08, met_meff_6 > 0.15*sqrt(GeV), meff_incl > 2200*GeV}); } /// Normalise counters after the run void finalize() { const double sf = 13.3*crossSection()/femtobarn/sumOfWeights(); scale(_h_2j_0800, sf); scale(_h_2j_1200, sf); scale(_h_2j_1600, sf); scale(_h_2j_2000, sf); scale(_h_3j_1200, sf); scale(_h_4j_1000, sf); scale(_h_4j_1400, sf); scale(_h_4j_1800, sf); scale(_h_4j_2200, sf); scale(_h_4j_2600, sf); scale(_h_5j_1400, sf); scale(_h_6j_1800, sf); scale(_h_6j_2200, sf); _flows.scale(sf); MSG_INFO("CUTFLOWS:\n\n" << _flows); } //@} private: /// @name Histograms //@{ CounterPtr _h_2j_0800, _h_2j_1200, _h_2j_1600, _h_2j_2000, _h_3j_1200; CounterPtr _h_4j_1000, _h_4j_1400, _h_4j_1800, _h_4j_2200, _h_4j_2600; CounterPtr _h_5j_1400, _h_6j_1800, _h_6j_2200; //@} /// Cut-flows Cutflows _flows; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_CONF_2016_078); } diff --git a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_094.cc b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_094.cc --- a/analyses/pluginATLAS/ATLAS_2016_CONF_2016_094.cc +++ b/analyses/pluginATLAS/ATLAS_2016_CONF_2016_094.cc @@ -1,167 +1,167 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedMET.hh" #include "Rivet/Tools/Cutflow.hh" namespace Rivet { /// @brief ATLAS 2016 1-lepton + many jets SUSY search, from 14.8/fb CONF note class ATLAS_2016_CONF_2016_094 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_CONF_2016_094); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections FinalState calofs(Cuts::abseta < 4.9); FastJets fj(calofs, FastJets::ANTIKT, 0.4); declare(fj, "TruthJets"); declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, [](const Jet& j) { if (j.abseta() > 2.5) return 0.; return j.bTagged(Cuts::pT > 5*GeV) ? 0.80 : j.cTagged(Cuts::pT > 5*GeV) ? 1/6. : 1/106.; }), "Jets"); // MissingMomentum mm(calofs); // declare(mm, "TruthMET"); // declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "MET"); FinalState es(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && !Cuts::absetaIn(1.37, 1.52) && Cuts::pT > 10*GeV); declare(es, "TruthElectrons"); - declare(SmearedParticles(es, ELECTRON_EFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons"); + declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "Electrons"); FinalState mus(Cuts::abspid == PID::MUON && Cuts::abseta < 2.4 && Cuts::pT > 10*GeV); declare(mus, "TruthMuons"); declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "Muons"); // Book histograms/counters book(_h_08j40_0b,"08j40_0b"); book(_h_09j40_0b,"09j40_0b"); book(_h_10j40_0b,"10j40_0b"); book(_h_08j40_3b,"08j40_3b"); book(_h_09j40_3b,"09j40_3b"); book(_h_10j40_3b,"10j40_3b"); book(_h_08j60_0b,"08j60_0b"); book(_h_09j60_0b,"09j60_0b"); book(_h_10j60_0b,"10j60_0b"); book(_h_08j60_3b,"08j60_3b"); book(_h_09j60_3b,"09j60_3b"); book(_h_10j60_3b,"10j60_3b"); } /// Perform the per-event analysis void analyze(const Event& event) { // Get baseline electrons, muons, and jets // NB. for electrons, we don't apply the loose ID here, since we don't want to double-count effs with later use of tight ID Particles elecs = apply(event, "Electrons").particles(); Particles muons = apply(event, "Muons").particles(); Jets jets = apply(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.4); ifilter_select(jets, JetEffFilter([](const Jet& j) { return j.pT() > 60*GeV ? 1.0 : 0.94; })); // Jet/electron/muon overlap removal and selection // Remove any untagged jet within dR = 0.2 of an electron for (const Particle& e : elecs) ifilter_discard(jets, [&](const Jet& j) { return !j.bTagged(Cuts::pT > 5*GeV) && deltaR(e, j, RAPIDITY) < 0.2; }); // Remove any untagged low-multiplicity/muon-dominated jet within dR = 0.4 of a muon for (const Particle& m : muons) ifilter_discard(jets, [&](const Jet& j) { if (!j.bTagged(Cuts::pT > 5*GeV)) return false; /// @note A different b-tag working point, 85%, was actually used here *sigh* if (deltaR(m, j, RAPIDITY) > 0.4) return false; if (j.particles(Cuts::abscharge != 0).size() < 3) return true; return m.pT()/j.pT() > 0.5; }); // Removing leptons within dR = 0.4 of remaining jets for (const Jet& j : jets) { ifilter_discard(elecs, deltaRLess(j, 0.4, RAPIDITY)); ifilter_discard(muons, deltaRLess(j, 0.4, RAPIDITY)); } // Signal jet and lepton selection const Jets sigjets40 = filter_select(jets, Cuts::pT > 40*GeV); const Jets sigjets60 = filter_select(sigjets40, Cuts::pT > 60*GeV); const Jets sigbjets40 = filter_select(sigjets40, [](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); }); const Jets sigbjets60 = filter_select(sigjets60, [](const Jet& j) { return j.bTagged(Cuts::pT > 5*GeV); }); const Particles sigmuons = filter_select(muons, Cuts::pT > 35*GeV); Particles sigelecs = filter_select(elecs, Cuts::pT > 35*GeV); - ifilter_select(sigelecs, ParticleEffFilter(ELECTRON_IDEFF_ATLAS_RUN2_TIGHT)); + ifilter_select(sigelecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_TIGHT)); ////////////////// // Event selection cuts if (sigelecs.size() + sigmuons.size() != 1) vetoEvent; const Particle siglepton = sigelecs.empty() ? sigmuons.front() : sigelecs.front(); /// @note The note describes Nj = 5, 6, 7, 8, 9, >= 10 and Nb = 0, 1, 2, 3, >= 4 = 30 2D bins /// for each jet cut... but only provides data for six Nj = >= 8, 9, 10, Nb = 0, >= 3 bins. /// We just implement the latter for now. // Fill counters if (sigjets40.size() >= 8 && sigbjets40.empty()) _h_08j40_0b->fill(); if (sigjets40.size() >= 9 && sigbjets40.empty()) _h_09j40_0b->fill(); if (sigjets40.size() >= 10 && sigbjets40.empty()) _h_10j40_0b->fill(); if (sigjets40.size() >= 8 && sigbjets40.size() >= 3) _h_08j40_3b->fill(); if (sigjets40.size() >= 9 && sigbjets40.size() >= 3) _h_09j40_3b->fill(); if (sigjets40.size() >= 10 && sigbjets40.size() >= 3) _h_10j40_3b->fill(); if (sigjets60.size() >= 8 && sigbjets60.empty()) _h_08j60_0b->fill(); if (sigjets60.size() >= 9 && sigbjets60.empty()) _h_09j60_0b->fill(); if (sigjets60.size() >= 10 && sigbjets60.empty()) _h_10j60_0b->fill(); if (sigjets60.size() >= 8 && sigbjets60.size() >= 3) _h_08j60_3b->fill(); if (sigjets60.size() >= 9 && sigbjets60.size() >= 3) _h_09j60_3b->fill(); if (sigjets60.size() >= 10 && sigbjets60.size() >= 3) _h_10j60_3b->fill(); } /// Normalise counters after the run void finalize() { const double sf = 14.8*crossSection()/femtobarn/sumOfWeights(); scale(_h_08j40_0b, sf); scale(_h_09j40_0b, sf); scale(_h_10j40_0b, sf); scale(_h_08j40_3b, sf); scale(_h_09j40_3b, sf); scale(_h_10j40_3b, sf); scale(_h_08j60_0b, sf); scale(_h_09j60_0b, sf); scale(_h_10j60_0b, sf); scale(_h_08j60_3b, sf); scale(_h_09j60_3b, sf); scale(_h_10j60_3b, sf); } //@} private: /// @name Histograms //@{ CounterPtr _h_08j40_0b, _h_09j40_0b, _h_10j40_0b, _h_08j40_3b, _h_09j40_3b, _h_10j40_3b; CounterPtr _h_08j60_0b, _h_09j60_0b, _h_10j60_0b, _h_08j60_3b, _h_09j60_3b, _h_10j60_3b; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_CONF_2016_094); } diff --git a/analyses/pluginATLAS/ATLAS_2016_I1452559.cc b/analyses/pluginATLAS/ATLAS_2016_I1452559.cc --- a/analyses/pluginATLAS/ATLAS_2016_I1452559.cc +++ b/analyses/pluginATLAS/ATLAS_2016_I1452559.cc @@ -1,130 +1,130 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedMET.hh" namespace Rivet { /// ATLAS 13 TeV monojet search with 3.2/fb of pp data class ATLAS_2016_I1452559 : public Analysis { public: DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_I1452559); void init() { FastJets jets(FinalState(Cuts::abseta < 4.9), FastJets::ANTIKT, 0.4); SmearedJets recojets(jets, JET_SMEAR_ATLAS_RUN1); declare(recojets, "Jets"); FinalState electrons(Cuts::abspid == PID::ELECTRON && Cuts::abseta < 2.47 && Cuts::pT > 20*GeV); - SmearedParticles recoelectrons(electrons, ELECTRON_EFF_ATLAS_RUN1); + SmearedParticles recoelectrons(electrons, ELECTRON_RECOEFF_ATLAS_RUN1); declare(recoelectrons, "Electrons"); FinalState muons(Cuts::abspid == PID::MUON && Cuts::abseta < 2.50 && Cuts::pT > 10*GeV); - SmearedParticles recomuons(muons, MUON_EFF_ATLAS_RUN1); + SmearedParticles recomuons(muons, MUON_RECOEFF_ATLAS_RUN1); declare(recomuons, "Muons"); VisibleFinalState calofs(Cuts::abseta < 4.9 && Cuts::abspid != PID::MUON); MissingMomentum met(calofs); SmearedMET recomet(met, MET_SMEAR_ATLAS_RUN1); declare(recomet, "MET"); /// Book histograms for (size_t i = 0; i < 7; ++i) book(_count_IM[i], "count_IM" + toString(i+1)); for (size_t i = 0; i < 6; ++i) book(_count_EM[i], "count_EM" + toString(i+1)); } void analyze(const Event& event) { const Jets jets = apply(event, "Jets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8); const Particles elecs = apply(event, "Electrons").particlesByPt(); const Particles mus = apply(event, "Muons").particlesByPt(); MSG_DEBUG("Number of raw jets, electrons, muons = " << jets.size() << ", " << elecs.size() << ", " << mus.size()); // Discard jets very close to electrons, or with low track multiplicity and close to muons const Jets isojets = filter_discard(jets, [&](const Jet& j) { /// @todo Add track efficiency random filtering if (any(elecs, deltaRLess(j, 0.2))) return true; if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.4*GeV).size() < 3 && any(mus, deltaRLess(j, 0.4))) return true; return false; }); // Discard electrons close to remaining jets const Particles isoelecs = filter_discard(elecs, [&](const Particle& e) { return any(isojets, deltaRLess(e, 0.4)); }); // Discard muons close to remaining jets const Particles isomus = filter_discard(mus, [&](const Particle& m) { for (const Jet& j : isojets) { if (deltaR(j,m) > 0.4) continue; if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 0.4*GeV).size() > 3) return true; } return false; }); // Calculate ETmiss //const Vector3& vet = apply(event, "MET").vectorEt(); const Vector3& vet = apply(event, "MET").vectorEt(); const double etmiss = vet.perp(); // Event selection cuts if (etmiss < 250*GeV) vetoEvent; // Require at least one jet with pT > 250 GeV and |eta| < 2.4 if (filter_select(isojets, Cuts::pT > 250*GeV && Cuts::abseta < 2.4).empty()) vetoEvent; // Require at most 4 jets with pT > 30 GeV and |eta| < 2.8 if (filter_select(isojets, Cuts::pT > 30*GeV).size() > 4) vetoEvent; // Require no isolated jets within |dphi| < 0.4 of the MET vector if (any(isojets, deltaPhiLess(-vet, 0.4))) vetoEvent; // Require no isolated electrons or muons if (!isoelecs.empty() || !isomus.empty()) vetoEvent; //////////////////// // Get ETmiss bin number and fill counters const int i_etmiss = binIndex(etmiss/GeV, ETMISS_CUTS); // Inclusive ETmiss bins for (int ibin = 0; ibin < 7; ++ibin) if (i_etmiss >= ibin) _count_IM[ibin]->fill(); // Exclusive ETmiss bins if (inRange(i_etmiss, 0, 6)) _count_EM[i_etmiss]->fill(); } void finalize() { const double norm = 3.2*crossSection()/femtobarn; scale(_count_IM, norm/sumOfWeights()); scale(_count_EM, norm/sumOfWeights()); } private: const vector ETMISS_CUTS = { 250, 300, 350, 400, 500, 600, 700, 13000 }; CounterPtr _count_IM[7], _count_EM[6]; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_I1452559); } diff --git a/analyses/pluginATLAS/ATLAS_2016_I1458270.cc b/analyses/pluginATLAS/ATLAS_2016_I1458270.cc --- a/analyses/pluginATLAS/ATLAS_2016_I1458270.cc +++ b/analyses/pluginATLAS/ATLAS_2016_I1458270.cc @@ -1,209 +1,209 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/Sphericity.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedMET.hh" #include "Rivet/Tools/Cutflow.hh" namespace Rivet { /// @brief ATLAS 0-lepton SUSY search with 3.2/fb of 13 TeV pp data class ATLAS_2016_I1458270 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_I1458270); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections FinalState calofs(Cuts::abseta < 4.8); FastJets fj(calofs, FastJets::ANTIKT, 0.4); declare(fj, "TruthJets"); declare(SmearedJets(fj, JET_SMEAR_ATLAS_RUN2, JET_BTAG_ATLAS_RUN2_MV2C20), "RecoJets"); MissingMomentum mm(calofs); declare(mm, "TruthMET"); declare(SmearedMET(mm, MET_SMEAR_ATLAS_RUN2), "RecoMET"); PromptFinalState es(Cuts::abseta < 2.47 && Cuts::abspid == PID::ELECTRON, true, true); declare(es, "TruthElectrons"); - declare(SmearedParticles(es, ELECTRON_EFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "RecoElectrons"); + declare(SmearedParticles(es, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2), "RecoElectrons"); PromptFinalState mus(Cuts::abseta < 2.7 && Cuts::abspid == PID::MUON, true); declare(mus, "TruthMuons"); declare(SmearedParticles(mus, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2), "RecoMuons"); // Book histograms/counters book(_h_2jl, "2jl"); book(_h_2jm, "2jm"); book(_h_2jt, "2jt"); book(_h_4jt, "4jt"); book(_h_5j , "5j"); book(_h_6jm, "6jm"); book(_h_6jt, "6jt"); // Book cut-flows const vector cuts2j = {"Pre-sel+MET+pT1", "Njet", "Dphi_min(j,MET)", "pT2", "MET/sqrtHT", "m_eff(incl)"}; _flows.addCutflow("2jl", cuts2j); _flows.addCutflow("2jm", cuts2j); _flows.addCutflow("2jt", cuts2j); const vector cutsXj = {"Pre-sel+MET+pT1", "Njet", "Dphi_min(j,MET)", "pT2", "pT4", "Aplanarity", "MET/m_eff(Nj)", "m_eff(incl)"}; _flows.addCutflow("4jt", cutsXj); _flows.addCutflow("5j", cutsXj); _flows.addCutflow("6jm", cutsXj); _flows.addCutflow("6jt", cutsXj); } /// Perform the per-event analysis void analyze(const Event& event) { _flows.fillinit(); // Same MET cut for all signal regions //const Vector3 vmet = -apply(event, "TruthMET").vectorEt(); const Vector3 vmet = -apply(event, "RecoMET").vectorEt(); const double met = vmet.mod(); if (met < 200*GeV) vetoEvent; // Get baseline electrons, muons, and jets Particles elecs = apply(event, "RecoElectrons").particles(Cuts::pT > 10*GeV); Particles muons = apply(event, "RecoMuons").particles(Cuts::pT > 10*GeV); Jets jets = apply(event, "RecoJets").jetsByPt(Cuts::pT > 20*GeV && Cuts::abseta < 2.8); ///< @todo Pile-up subtraction // Jet/electron/muons overlap removal and selection // Remove any |eta| < 2.8 jet within dR = 0.2 of a baseline electron for (const Particle& e : elecs) ifilter_discard(jets, deltaRLess(e, 0.2, RAPIDITY)); // Remove any electron or muon with dR < 0.4 of a remaining (Nch > 3) jet for (const Jet& j : jets) { /// @todo Add track efficiency random filtering ifilter_discard(elecs, deltaRLess(j, 0.4, RAPIDITY)); if (j.particles(Cuts::abscharge > 0 && Cuts::pT > 500*MeV).size() >= 3) ifilter_discard(muons, deltaRLess(j, 0.4, RAPIDITY)); } // Discard the softer of any electrons within dR < 0.05 for (size_t i = 0; i < elecs.size(); ++i) { const Particle& e1 = elecs[i]; /// @todo Would be nice to pass a "tail view" for the filtering, but awkward without range API / iterator guts ifilter_discard(elecs, [&](const Particle& e2){ return e2.pT() < e1.pT() && deltaR(e1,e2) < 0.05; }); } // Loose electron selection - ifilter_select(elecs, ParticleEffFilter(ELECTRON_IDEFF_ATLAS_RUN2_LOOSE)); + ifilter_select(elecs, ParticleEffFilter(ELECTRON_EFF_ATLAS_RUN2_LOOSE)); // Veto the event if there are any remaining baseline leptons if (!elecs.empty()) vetoEvent; if (!muons.empty()) vetoEvent; // Signal jets have pT > 50 GeV const Jets jets50 = filter_select(jets, Cuts::pT > 50*GeV); if (jets50.size() < 2) vetoEvent; vector jetpts; transform(jets, jetpts, pT); vector jetpts50; transform(jets50, jetpts50, pT); const double j1pt = jetpts50[0]; const double j2pt = jetpts50[1]; if (j1pt < 200*GeV) vetoEvent; // Construct multi-jet observables const double ht = sum(jetpts, 0.0); const double met_sqrt_ht = met / sqrt(ht); const double meff_incl = sum(jetpts50, met); // Get dphis between MET and jets vector dphimets50; transform(jets50, dphimets50, deltaPhiWRT(vmet)); const double min_dphi_met_3 = min(head(dphimets50, 3)); MSG_DEBUG(dphimets50 << ", " << min_dphi_met_3); // Jet aplanarity Sphericity sph; sph.calc(jets); const double aplanarity = sph.aplanarity(); // Fill SR counters // 2-jet SRs if (_flows["2jl"].filltail({true, true, min_dphi_met_3 > 0.8, j2pt > 200*GeV, met_sqrt_ht > 15*sqrt(GeV), meff_incl > 1200*GeV})) _h_2jl->fill(); if (_flows["2jm"].filltail({j1pt > 300*GeV, true, min_dphi_met_3 > 0.4, j2pt > 50*GeV, met_sqrt_ht > 15*sqrt(GeV), meff_incl > 1600*GeV})) _h_2jm->fill(); if (_flows["2jt"].filltail({true, true, min_dphi_met_3 > 0.8, j2pt > 200*GeV, met_sqrt_ht > 20*sqrt(GeV), meff_incl > 2000*GeV})) _h_2jt->fill(); // Upper multiplicity SRs const double j4pt = jets50.size() > 3 ? jetpts50[3] : -1; const double j5pt = jets50.size() > 4 ? jetpts50[4] : -1; const double j6pt = jets50.size() > 5 ? jetpts50[5] : -1; const double meff_4 = jets50.size() > 3 ? sum(head(jetpts50, 4), met) : -1; const double meff_5 = jets50.size() > 4 ? meff_4 + jetpts50[4] : -1; const double meff_6 = jets50.size() > 5 ? meff_5 + jetpts50[5] : -1; const double met_meff_4 = met / meff_4; const double met_meff_5 = met / meff_5; const double met_meff_6 = met / meff_6; const double min_dphi_met_more = jets50.size() > 3 ? min(tail(dphimets50, -3)) : -1; if (_flows["4jt"].filltail({true, jets50.size() >= 4, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2, jetpts[1] > 100*GeV, j4pt > 100*GeV, aplanarity > 0.04, met_meff_4 > 0.20, meff_incl > 2200*GeV})) _h_4jt->fill(); if (_flows["5j"].filltail({true, jets50.size() >= 5, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2, jetpts[1] > 100*GeV, j4pt > 100*GeV && j5pt > 50*GeV, aplanarity > 0.04, met_meff_5 > 0.25, meff_incl > 1600*GeV})) _h_5j->fill(); if (_flows["6jm"].filltail({true, jets50.size() >= 6, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2, jetpts[1] > 100*GeV, j4pt > 100*GeV && j6pt > 50*GeV, aplanarity > 0.04, met_meff_6 > 0.25, meff_incl > 1600*GeV})) _h_6jm->fill(); if (_flows["6jt"].filltail({true, jets50.size() >= 6, min_dphi_met_3 > 0.4 && min_dphi_met_more > 0.2, jetpts[1] > 100*GeV, j4pt > 100*GeV && j6pt > 50*GeV, aplanarity > 0.04, met_meff_6 > 0.20, meff_incl > 2000*GeV})) _h_6jt->fill(); } /// Normalise histograms etc., after the run void finalize() { const double sf = 3.2*crossSection()/femtobarn/sumOfWeights(); scale(_h_2jl, sf); scale(_h_2jl, sf); scale(_h_2jl, sf); scale(_h_4jt, sf); scale(_h_5j, sf); scale(_h_6jm, sf); scale(_h_6jt, sf); MSG_INFO("CUTFLOWS:\n\n" << _flows); } //@} private: /// @name Histograms //@{ CounterPtr _h_2jl, _h_2jm, _h_2jt; CounterPtr _h_4jt, _h_5j; CounterPtr _h_6jm, _h_6jt; //@} /// Cut-flows Cutflows _flows; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_I1458270); } diff --git a/analyses/pluginMC/EXAMPLE_SMEAR.cc b/analyses/pluginMC/EXAMPLE_SMEAR.cc --- a/analyses/pluginMC/EXAMPLE_SMEAR.cc +++ b/analyses/pluginMC/EXAMPLE_SMEAR.cc @@ -1,251 +1,251 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/TauFinder.hh" #include "Rivet/Projections/SmearedJets.hh" #include "Rivet/Projections/SmearedParticles.hh" #include "Rivet/Projections/SmearedMET.hh" namespace Rivet { class EXAMPLE_SMEAR : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(EXAMPLE_SMEAR); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { MissingMomentum mm(Cuts::abseta < 5); declare(mm, "MET0"); SmearedMET smm1(mm, MET_SMEAR_ATLAS_RUN2); declare(smm1, "MET1"); SmearedMET smm2(mm, [](const Vector3& met, double){ return P3_SMEAR_LEN_GAUSS(met, 0.1*met.mod()); }); declare(smm2, "MET2"); FastJets fj(FinalState(Cuts::abseta < 5), FastJets::ANTIKT, 0.4); declare(fj, "Jets0"); SmearedJets sj1(fj, JET_SMEAR_IDENTITY); declare(sj1, "Jets1"); SmearedJets sj2(fj, JET_SMEAR_ATLAS_RUN2, [](const Jet& j){ return j.bTagged() ? 0.7*(1 - exp(-j.pT()/(10*GeV))) : 0.01; } ); declare(sj2, "Jets2"); SmearedJets sj3(fj, JET_SMEAR_CMS_RUN2, JET_BTAG_EFFS(0.7, 0.1, 0.01), JET_CTAG_PERFECT, JET_EFF_CONST(0.8)); declare(sj3, "Jets3"); IdentifiedFinalState photons(Cuts::abseta < 5, PID::PHOTON); IdentifiedFinalState truthelectrons(Cuts::abseta < 5 && Cuts::pT > 10*GeV, {{PID::ELECTRON, PID::POSITRON}}); declare(truthelectrons, "Electrons0"); DressedLeptons dressedelectrons(photons, truthelectrons, 0.2); declare(dressedelectrons, "Electrons1"); - SmearedParticles recoelectrons(dressedelectrons, ELECTRON_EFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2); + SmearedParticles recoelectrons(dressedelectrons, ELECTRON_RECOEFF_ATLAS_RUN2, ELECTRON_SMEAR_ATLAS_RUN2); declare(recoelectrons, "Electrons2"); IdentifiedFinalState truthmuons(Cuts::abseta < 5 && Cuts::pT > 10*GeV, {{PID::MUON, PID::ANTIMUON}}); declare(truthmuons, "Muons0"); DressedLeptons dressedmuons(photons, truthmuons, 0.2); declare(dressedmuons, "Muons1"); SmearedParticles recomuons(dressedmuons, MUON_EFF_ATLAS_RUN2, MUON_SMEAR_ATLAS_RUN2); declare(recomuons, "Muons2"); TauFinder truthtaus(TauFinder::DecayMode::ANY, Cuts::abseta < 5 && Cuts::pT > 10*GeV); declare(truthtaus, "Taus0"); DressedLeptons dressedtaus(photons, truthtaus, 0.2); declare(dressedtaus, "Taus1"); SmearedParticles recotaus(dressedtaus, TAU_EFF_ATLAS_RUN2, TAU_SMEAR_ATLAS_RUN2); declare(recotaus, "Taus2"); book(_h_met_true ,"met_true", 30, 0.0, 120); book(_h_met_reco ,"met_reco", 30, 0.0, 120); book(_h_nj_true ,"jet_N_true", 10, -0.5, 9.5); book(_h_nj_reco ,"jet_N_reco", 10, -0.5, 9.5); book(_h_j1pt_true ,"jet_pt1_true", 30, 0.0, 120); book(_h_j1pt_reco ,"jet_pt1_reco", 30, 0.0, 120); book(_h_j1eta_true ,"jet_eta1_true", 20, -5.0, 5.0); book(_h_j1eta_reco ,"jet_eta1_reco", 20, -5.0, 5.0); book(_h_ne_true ,"elec_N_true", 5, -0.5, 4.5); book(_h_ne_reco ,"elec_N_reco", 5, -0.5, 4.5); book(_h_e1pt_true ,"elec_pt1_true", 30, 0, 120); book(_h_e1pt_reco ,"elec_pt1_reco", 30, 0, 120); book(_h_e1eta_true ,"elec_eta1_true", 20, -5.0, 5.0); book(_h_e1eta_reco ,"elec_eta1_reco", 20, -5.0, 5.0); book(_h_nm_true ,"muon_N_true", 5, -0.5, 4.5); book(_h_nm_reco ,"muon_N_reco", 5, -0.5, 4.5); book(_h_m1pt_true ,"muon_pt1_true", 30, 0, 120); book(_h_m1pt_reco ,"muon_pt1_reco", 30, 0, 120); book(_h_m1eta_true ,"muon_eta1_true", 20, -5.0, 5.0); book(_h_m1eta_reco ,"muon_eta1_reco", 20, -5.0, 5.0); book(_h_nt_true ,"tau_N_true", 5, -0.5, 4.5); book(_h_nt_reco ,"tau_N_reco", 5, -0.5, 4.5); book(_h_t1pt_true ,"tau_pt1_true", 30, 0, 120); book(_h_t1pt_reco ,"tau_pt1_reco", 30, 0, 120); book(_h_t1eta_true ,"tau_eta1_true", 20, -5.0, 5.0); book(_h_t1eta_reco ,"tau_eta1_reco", 20, -5.0, 5.0); } /// Perform the per-event analysis void analyze(const Event& event) { const double weight = 1.0; const Vector3 met0 = apply(event, "MET0").vectorEt(); const Vector3 met1 = apply(event, "MET1").vectorEt(); const Vector3 met2 = apply(event, "MET2").vectorEt(); MSG_DEBUG("MET = " << met0.mod()/GeV << ", " << met1.mod()/GeV << ", " << met2.mod()/GeV << " GeV"); _h_met_true->fill(met0.mod()/GeV, weight); _h_met_reco->fill(met1.mod()/GeV, weight); if (met0.perp() > 0 && met1.perp() > 0 && deltaPhi(met0, met1) > 0.1) { MSG_WARNING("Large MET phi change: " << met0.phi() << " -> " << met1.phi() << "; dphi = " << deltaPhi(met0, met1)); } const Jets jets0 = apply(event, "Jets0").jetsByPt(Cuts::pT > 10*GeV); const Jets jets1 = apply(event, "Jets1").jetsByPt(Cuts::pT > 10*GeV); const Jets jets2 = apply(event, "Jets2").jetsByPt(Cuts::pT > 10*GeV); const Jets jets3 = apply(event, "Jets3").jetsByPt(Cuts::pT > 10*GeV); MSG_DEBUG("Numbers of jets = " << jets0.size() << " true; " << jets1.size() << ", " << jets2.size() << ", " << jets3.size()); if (!jets0.empty() && !jets2.empty() && deltaPhi(jets0[0], jets2[0]) > 0.1) { MSG_DEBUG("Large jet1 phi change (could be a different truth-jet): " << jets0[0].phi() << " -> " << jets2[0].phi() << "; dphi = " << deltaPhi(jets0[0], jets2[0]) << "; pT = " << jets0[0].pT()/GeV << " -> " << jets2[0].pT()/GeV); } _h_nj_true->fill(jets0.size(), weight); _h_nj_reco->fill(jets2.size(), weight); if (!jets0.empty()) { _h_j1pt_true->fill(jets0.front().pT()/GeV, weight); _h_j1eta_true->fill(jets0.front().eta(), weight); } if (!jets2.empty()) { _h_j1pt_reco->fill(jets2.front().pT()/GeV, weight); _h_j1eta_reco->fill(jets2.front().eta(), weight); } const Particles& elecs1 = apply(event, "Electrons1").particlesByPt(); const Particles& elecs2 = apply(event, "Electrons2").particlesByPt(); MSG_DEBUG("Numbers of electrons = " << elecs1.size() << " true; " << elecs2.size() << " reco"); _h_ne_true->fill(elecs1.size(), weight); _h_ne_reco->fill(elecs2.size(), weight); if (!elecs1.empty()) { _h_e1pt_true->fill(elecs1.front().pT()/GeV, weight); _h_e1eta_true->fill(elecs1.front().eta(), weight); } if (!elecs2.empty()) { _h_e1pt_reco->fill(elecs2.front().pT()/GeV, weight); _h_e1eta_reco->fill(elecs2.front().eta(), weight); } const Particles& muons1 = apply(event, "Muons1").particlesByPt(); const Particles& muons2 = apply(event, "Muons2").particlesByPt(); MSG_DEBUG("Numbers of muons = " << muons1.size() << " true; " << muons2.size() << " reco"); _h_nm_true->fill(muons1.size(), weight); _h_nm_reco->fill(muons2.size(), weight); if (!muons1.empty()) { _h_m1pt_true->fill(muons1.front().pT()/GeV, weight); _h_m1eta_true->fill(muons1.front().eta(), weight); } if (!muons2.empty()) { _h_m1pt_reco->fill(muons2.front().pT()/GeV, weight); _h_m1eta_reco->fill(muons2.front().eta(), weight); } const Particles& taus1 = apply(event, "Taus1").particlesByPt(); const Particles& taus2 = apply(event, "Taus2").particlesByPt(); MSG_DEBUG("Numbers of taus = " << taus1.size() << " true; " << taus2.size() << " reco"); _h_nt_true->fill(taus1.size(), weight); _h_nt_reco->fill(taus2.size(), weight); if (!taus1.empty()) { _h_t1pt_true->fill(taus1.front().pT()/GeV, weight); _h_t1eta_true->fill(taus1.front().eta(), weight); } if (!taus2.empty()) { _h_t1pt_reco->fill(taus2.front().pT()/GeV, weight); _h_t1eta_reco->fill(taus2.front().eta(), weight); } } /// Normalise histograms etc., after the run void finalize() { normalize(_h_met_true); normalize(_h_met_reco); normalize(_h_nj_true); normalize(_h_nj_reco); normalize(_h_j1pt_true, 1-_h_nj_true->bin(0).area()); normalize(_h_j1pt_reco, 1-_h_nj_reco->bin(0).area()); normalize(_h_j1eta_true, 1-_h_nj_true->bin(0).area()); normalize(_h_j1eta_reco, 1-_h_nj_reco->bin(0).area()); normalize(_h_ne_true); normalize(_h_ne_reco); normalize(_h_e1pt_true, 1-_h_ne_true->bin(0).area()); normalize(_h_e1pt_reco, 1-_h_ne_reco->bin(0).area()); normalize(_h_e1eta_true, 1-_h_ne_true->bin(0).area()); normalize(_h_e1eta_reco, 1-_h_ne_reco->bin(0).area()); normalize(_h_nm_true); normalize(_h_nm_reco); normalize(_h_m1pt_true, 1-_h_nm_true->bin(0).area()); normalize(_h_m1pt_reco, 1-_h_nm_reco->bin(0).area()); normalize(_h_m1eta_true, 1-_h_nm_true->bin(0).area()); normalize(_h_m1eta_reco, 1-_h_nm_reco->bin(0).area()); normalize(_h_nt_true); normalize(_h_nt_reco); normalize(_h_t1pt_true, 1-_h_nt_true->bin(0).area()); normalize(_h_t1pt_reco, 1-_h_nt_reco->bin(0).area()); normalize(_h_t1eta_true, 1-_h_nt_true->bin(0).area()); normalize(_h_t1eta_reco, 1-_h_nt_reco->bin(0).area()); } //@} private: /// @name Histograms //@{ Histo1DPtr _h_met_true, _h_met_reco; Histo1DPtr _h_nj_true, _h_nj_reco, _h_ne_true, _h_ne_reco, _h_nm_true, _h_nm_reco, _h_nt_true, _h_nt_reco; Histo1DPtr _h_j1pt_true, _h_j1pt_reco, _h_e1pt_true, _h_e1pt_reco, _h_m1pt_true, _h_m1pt_reco, _h_t1pt_true, _h_t1pt_reco; Histo1DPtr _h_j1eta_true, _h_j1eta_reco, _h_e1eta_true, _h_e1eta_reco, _h_m1eta_true, _h_m1eta_reco, _h_t1eta_true, _h_t1eta_reco; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(EXAMPLE_SMEAR); } diff --git a/src/Projections/FastJets.cc b/src/Projections/FastJets.cc --- a/src/Projections/FastJets.cc +++ b/src/Projections/FastJets.cc @@ -1,221 +1,221 @@ // -*- C++ -*- #include "Rivet/Config/RivetCommon.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/HeavyHadrons.hh" #include "Rivet/Projections/TauFinder.hh" namespace Rivet { void FastJets::_initBase() { setName("FastJets"); declare(HeavyHadrons(), "HFHadrons"); declare(TauFinder(TauFinder::DecayMode::HADRONIC), "Taus"); // Print/hide FJ banner - cout.setstate(std::ios_base::badbit); + std::cout.setstate(std::ios_base::badbit); fastjet::ClusterSequence::print_banner(); std::cout.clear(); } void FastJets::_initJdef(Algo alg, double rparameter, double seed_threshold) { MSG_DEBUG("JetAlg = " << static_cast(alg)); MSG_DEBUG("R parameter = " << rparameter); MSG_DEBUG("Seed threshold = " << seed_threshold); if (alg == KT) { _jdef = fastjet::JetDefinition(fastjet::kt_algorithm, rparameter, fastjet::E_scheme); } else if (alg == CAM) { _jdef = fastjet::JetDefinition(fastjet::cambridge_algorithm, rparameter, fastjet::E_scheme); } else if (alg == ANTIKT) { _jdef = fastjet::JetDefinition(fastjet::antikt_algorithm, rparameter, fastjet::E_scheme); } else if (alg == DURHAM) { _jdef = fastjet::JetDefinition(fastjet::ee_kt_algorithm, fastjet::E_scheme); } else if (alg == GENKTEE) { _jdef = fastjet::JetDefinition(fastjet::ee_genkt_algorithm, rparameter, -1); } else { // Plugins: if (alg == SISCONE) { const double OVERLAP_THRESHOLD = 0.75; _plugin.reset(new fastjet::SISConePlugin(rparameter, OVERLAP_THRESHOLD)); } else if (alg == PXCONE) { string msg = "Using own c++ version of PxCone, since FastJet doesn't install it by default. "; msg += "Please notify the Rivet authors if this behaviour should be changed."; MSG_WARNING(msg); // _plugin.reset(new fastjet::PxConePlugin(rparameter)); _plugin.reset(new Rivet::PxConePlugin(rparameter)); } else if (alg == ATLASCONE) { const double OVERLAP_THRESHOLD = 0.5; _plugin.reset(new fastjet::ATLASConePlugin(rparameter, seed_threshold, OVERLAP_THRESHOLD)); } else if (alg == CMSCONE) { _plugin.reset(new fastjet::CMSIterativeConePlugin(rparameter, seed_threshold)); } else if (alg == CDFJETCLU) { const double OVERLAP_THRESHOLD = 0.75; _plugin.reset(new fastjet::CDFJetCluPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == CDFMIDPOINT) { const double OVERLAP_THRESHOLD = 0.5; _plugin.reset(new fastjet::CDFMidPointPlugin(rparameter, OVERLAP_THRESHOLD, seed_threshold)); } else if (alg == D0ILCONE) { const double min_jet_Et = 6.0; _plugin.reset(new fastjet::D0RunIIConePlugin(rparameter, min_jet_Et)); } else if (alg == JADE) { _plugin.reset(new fastjet::JadePlugin()); } else if (alg == TRACKJET) { _plugin.reset(new fastjet::TrackJetPlugin(rparameter)); } _jdef = fastjet::JetDefinition(_plugin.get()); } } CmpState FastJets::compare(const Projection& p) const { const FastJets& other = dynamic_cast(p); return \ cmp(_useMuons, other._useMuons) || cmp(_useInvisibles, other._useInvisibles) || mkNamedPCmp(other, "FS") || cmp(_jdef.jet_algorithm(), other._jdef.jet_algorithm()) || cmp(_jdef.recombination_scheme(), other._jdef.recombination_scheme()) || cmp(_jdef.plugin(), other._jdef.plugin()) || cmp(_jdef.R(), other._jdef.R()) || cmp(_adef, other._adef); } // STATIC PseudoJets FastJets::mkClusterInputs(const Particles& fsparticles, const Particles& tagparticles) { PseudoJets pjs; /// @todo Use FastJet3's UserInfo system to store Particle pointers directly? // Store 4 vector data about each particle into FastJet's PseudoJets for (size_t i = 0; i < fsparticles.size(); ++i) { fastjet::PseudoJet pj = fsparticles[i]; pj.set_user_index(i+1); pjs.push_back(pj); } // And the same for ghost tagging particles (with negative user indices) for (size_t i = 0; i < tagparticles.size(); ++i) { fastjet::PseudoJet pj = tagparticles[i]; pj *= 1e-20; ///< Ghostify the momentum pj.set_user_index(-i-1); pjs.push_back(pj); } return pjs; } // STATIC Jet FastJets::mkJet(const PseudoJet& pj, const Particles& fsparticles, const Particles& tagparticles) { const PseudoJets pjconstituents = pj.constituents(); Particles constituents, tags; constituents.reserve(pjconstituents.size()); for (const fastjet::PseudoJet& pjc : pjconstituents) { // Pure ghosts don't have corresponding particles if (pjc.has_area() && pjc.is_pure_ghost()) continue; // Default user index = 0 doesn't give valid particle lookup if (pjc.user_index() == 0) continue; // Split by index sign into constituent & tag lookup if (pjc.user_index() > 0) { // Find constituents if index > 0 const size_t i = pjc.user_index() - 1; if (i >= fsparticles.size()) throw RangeError("FS particle lookup failed in jet construction"); constituents.push_back(fsparticles.at(i)); } else if (!tagparticles.empty()) { // Find tags if index < 0 const size_t i = abs(pjc.user_index()) - 1; if (i >= tagparticles.size()) throw RangeError("Tag particle lookup failed in jet construction"); tags.push_back(tagparticles.at(i)); } } return Jet(pj, constituents, tags); } // STATIC Jets FastJets::mkJets(const PseudoJets& pjs, const Particles& fsparticles, const Particles& tagparticles) { Jets rtn; rtn.reserve(pjs.size()); for (const PseudoJet pj : pjs) { rtn.push_back(FastJets::mkJet(pj, fsparticles, tagparticles)); } return rtn; } void FastJets::project(const Event& e) { // Assemble final state particles const string fskey = (_useInvisibles == JetAlg::Invisibles::NONE) ? "VFS" : "FS"; Particles fsparticles = applyProjection(e, fskey).particles(); // Remove prompt invisibles if needed (already done by VFS if using NO_INVISIBLES) if (_useInvisibles == JetAlg::Invisibles::DECAY) { ifilter_discard(fsparticles, [](const Particle& p) { return !p.isVisible() && p.isPrompt(); }); } // Remove prompt/all muons if needed if (_useMuons == JetAlg::Muons::DECAY) { ifilter_discard(fsparticles, [](const Particle& p) { return isMuon(p) && p.isPrompt(); }); } else if (_useMuons == JetAlg::Muons::NONE) { ifilter_discard(fsparticles, isMuon); } // Tagging particles const Particles chadrons = applyProjection(e, "HFHadrons").cHadrons(); const Particles bhadrons = applyProjection(e, "HFHadrons").bHadrons(); const Particles taus = applyProjection(e, "Taus").particles(); calc(fsparticles, chadrons+bhadrons+taus); } void FastJets::calc(const Particles& fsparticles, const Particles& tagparticles) { MSG_DEBUG("Finding jets from " << fsparticles.size() << " input particles + " << tagparticles.size() << " tagging particles"); _fsparticles = fsparticles; _tagparticles = tagparticles; // Make pseudojets, with mapping info to Rivet FS and tag particles PseudoJets pjs = mkClusterInputs(_fsparticles, _tagparticles); // Run either basic or area-calculating cluster sequence as reqd. if (_adef) { _cseq.reset(new fastjet::ClusterSequenceArea(pjs, _jdef, *_adef)); } else { _cseq.reset(new fastjet::ClusterSequence(pjs, _jdef)); } MSG_DEBUG("ClusterSequence constructed; Njets_tot = " << _cseq->inclusive_jets().size() << ", Njets(pT > 10 GeV) = " << _cseq->inclusive_jets(10*GeV).size()); } void FastJets::reset() { _yscales.clear(); _fsparticles.clear(); _tagparticles.clear(); /// @todo _cseq = fastjet::ClusterSequence(); } Jets FastJets::_jets() const { /// @todo Cache? return mkJets(pseudojets(), _fsparticles, _tagparticles); } /// @todo "Automate" trimming as part of project() with pre-registered Filters Jet FastJets::trimJet(const Jet& input, const fastjet::Filter& trimmer) const { if (input.pseudojet().associated_cluster_sequence() != clusterSeq().get()) throw Error("To trim a Rivet::Jet, its associated PseudoJet must have come from this FastJets' ClusterSequence"); PseudoJet pj = trimmer(input); return mkJet(pj, _fsparticles, _tagparticles); } PseudoJets FastJets::pseudoJets(double ptmin) const { return clusterSeq() ? clusterSeq()->inclusive_jets(ptmin) : PseudoJets(); } }