diff --git a/analyses/pluginCMS/CMS_2015_I1310737.cc b/analyses/pluginCMS/CMS_2015_I1310737.cc --- a/analyses/pluginCMS/CMS_2015_I1310737.cc +++ b/analyses/pluginCMS/CMS_2015_I1310737.cc @@ -1,193 +1,193 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/InvMassFinalState.hh" #include "Rivet/Projections/VisibleFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/ZFinder.hh" namespace Rivet { class CMS_2015_I1310737 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2015_I1310737); /// Book histograms and initialise projections before the run void init() { FinalState fs; ///< @todo No cuts? VisibleFinalState visfs(fs); // Prompt leptons only PromptFinalState pfs(fs); ZFinder zeeFinder(pfs, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV, PID::ELECTRON, 71.0*GeV, 111.0*GeV); declare(zeeFinder, "ZeeFinder"); ZFinder zmumuFinder(pfs, Cuts::abseta < 2.4 && Cuts::pT > 20*GeV, PID::MUON, 71.0*GeV, 111.0*GeV); declare(zmumuFinder, "ZmumuFinder"); VetoedFinalState jetConstits(visfs); jetConstits.addVetoOnThisFinalState(zeeFinder); jetConstits.addVetoOnThisFinalState(zmumuFinder); FastJets akt05Jets(jetConstits, FastJets::ANTIKT, 0.5); declare(akt05Jets, "AntiKt05Jets"); book(_h_excmult_jets_tot ,1, 1, 1); book(_h_incmult_jets_tot ,2, 1, 1); book(_h_leading_jet_pt_tot ,3, 1, 1); book(_h_second_jet_pt_tot ,4, 1, 1); book(_h_third_jet_pt_tot ,5, 1, 1); book(_h_fourth_jet_pt_tot ,6, 1, 1); book(_h_leading_jet_eta_tot ,7, 1, 1); book(_h_second_jet_eta_tot ,8, 1, 1); book(_h_third_jet_eta_tot ,9, 1, 1); book(_h_fourth_jet_eta_tot ,10, 1, 1); book(_h_ht1_tot ,11, 1, 1); book(_h_ht2_tot ,12, 1, 1); book(_h_ht3_tot ,13, 1, 1); book(_h_ht4_tot ,14, 1, 1); } /// Perform the per-event analysis void analyze(const Event& event) {; const ZFinder& zeeFS = apply(event, "ZeeFinder"); const ZFinder& zmumuFS = apply(event, "ZmumuFinder"); const Particles& zees = zeeFS.bosons(); const Particles& zmumus = zmumuFS.bosons(); // We did not find exactly one Z. No good. if (zees.size() + zmumus.size() != 1) { MSG_DEBUG("Did not find exactly one good Z candidate"); vetoEvent; } // Find the (dressed!) leptons const Particles& dressedLeptons = zees.size() ? zeeFS.constituents() : zmumuFS.constituents(); // Cluster jets // NB. Veto has already been applied on leptons and photons used for dressing const FastJets& fj = apply(event, "AntiKt05Jets"); const Jets& jets = fj.jetsByPt(Cuts::abseta < 2.4 && Cuts::pT > 30*GeV); // Perform lepton-jet overlap and HT calculation double ht = 0; Jets goodjets; for (const Jet& j : jets) { // Decide if this jet is "good", i.e. isolated from the leptons /// @todo Nice use-case for any() and a C++11 lambda bool overlap = false; for (const Particle& l : dressedLeptons) { - if (Rivet::deltaR(j, l) < 0.5) { + if (deltaR(j, l) < 0.5) { overlap = true; break; } } // Fill HT and good-jets collection if (overlap) continue; goodjets.push_back(j); ht += j.pT(); } // We don't care about events with no isolated jets if (goodjets.empty()) { MSG_DEBUG("No jets in event"); vetoEvent; } ///////////////// // Weight to be used for histo filling const double w = 0.5 * 1.0; // Fill jet number integral histograms _h_excmult_jets_tot->fill(goodjets.size(), w); /// @todo Could be better computed by toIntegral transform on exclusive histo for (size_t iJet = 1; iJet <= goodjets.size(); iJet++ ) _h_incmult_jets_tot->fill(iJet, w); // Fill leading jet histograms const Jet& j1 = goodjets[0]; _h_leading_jet_pt_tot->fill(j1.pT()/GeV, w); _h_leading_jet_eta_tot->fill(j1.abseta(), w); _h_ht1_tot->fill(ht/GeV, w); // Fill 2nd jet histograms if (goodjets.size() < 2) return; const Jet& j2 = goodjets[1]; _h_second_jet_pt_tot->fill(j2.pT()/GeV, w); _h_second_jet_eta_tot->fill(j2.abseta(), w); _h_ht2_tot->fill(ht/GeV, w); // Fill 3rd jet histograms if (goodjets.size() < 3) return; const Jet& j3 = goodjets[2]; _h_third_jet_pt_tot->fill(j3.pT()/GeV, w); _h_third_jet_eta_tot->fill(j3.abseta(), w); _h_ht3_tot->fill(ht/GeV, w); // Fill 4th jet histograms if (goodjets.size() < 4) return; const Jet& j4 = goodjets[3]; _h_fourth_jet_pt_tot->fill(j4.pT()/GeV, w); _h_fourth_jet_eta_tot->fill(j4.abseta(), w); _h_ht4_tot->fill(ht/GeV, w); } /// Normalise histograms etc., after the run void finalize() { const double norm = (sumOfWeights() != 0) ? crossSection()/sumOfWeights() : 1.0; MSG_INFO("Cross section = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << crossSection() << " pb"); MSG_INFO("# Events = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << numEvents() ); MSG_INFO("SumW = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(3) << sumOfWeights()); MSG_INFO("Norm factor = " << std::setfill(' ') << std::setw(14) << std::fixed << std::setprecision(6) << norm); scale(_h_excmult_jets_tot, norm ); scale(_h_incmult_jets_tot, norm ); scale(_h_leading_jet_pt_tot, norm ); scale(_h_second_jet_pt_tot, norm ); scale(_h_third_jet_pt_tot, norm ); scale(_h_fourth_jet_pt_tot, norm ); scale(_h_leading_jet_eta_tot, norm ); scale(_h_second_jet_eta_tot, norm ); scale(_h_third_jet_eta_tot, norm ); scale(_h_fourth_jet_eta_tot, norm ); scale(_h_ht1_tot, norm ); scale(_h_ht2_tot, norm ); scale(_h_ht3_tot, norm ); scale(_h_ht4_tot, norm ); } private: /// @name Histograms Histo1DPtr _h_excmult_jets_tot, _h_incmult_jets_tot; Histo1DPtr _h_leading_jet_pt_tot, _h_second_jet_pt_tot, _h_third_jet_pt_tot, _h_fourth_jet_pt_tot; Histo1DPtr _h_leading_jet_eta_tot, _h_second_jet_eta_tot, _h_third_jet_eta_tot, _h_fourth_jet_eta_tot; Histo1DPtr _h_ht1_tot, _h_ht2_tot, _h_ht3_tot, _h_ht4_tot; }; DECLARE_RIVET_PLUGIN(CMS_2015_I1310737); } 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,493 @@ // -*- 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 +namespace Rivet { - /// @brief Special dressed lepton finder + + /// @brief Special CMS dressed lepton finder /// /// Find dressed leptons by clustering all leptons and photons - class SpecialDressedLeptons : public FinalState { - public: + /// + /// @todo Merge functionality into main DressedLeptons, or avoid projection + class SpecialDressedLeptons_CMS_2016_I1491950 : public FinalState { + public: + /// The default constructor. May specify cuts - SpecialDressedLeptons(const FinalState& fs, const Cut& cut) + SpecialDressedLeptons_CMS_2016_I1491950(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)); + return unique_ptr(new SpecialDressedLeptons_CMS_2016_I1491950(*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.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.pid() == 0)) continue; - + DressedLepton lepton = DressedLepton(lepCand); - + for (const Particle& cand : jet.particles()) { + //if (isSame(cand, lepCand)) continue; if (cand == lepCand) continue; - if (cand.pid() != PID::PHOTON) 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: + + 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); + SpecialDressedLeptons_CMS_2016_I1491950 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) - { + void analyze(const Event& event) { // leptons - const SpecialDressedLeptons& dressedleptons_proj = applyProjection(event, "DressedLeptons"); + const auto& dressedleptons_proj = applyProjection(event, "DressedLeptons"); std::vector dressedLeptons = dressedleptons_proj.dressedLeptons(); - - if(dressedLeptons.size() != 1) return; - + if (dressedLeptons.size() != 1) vetoEvent; + // neutrinos const Particles neutrinos = applyProjection(event, "Neutrinos").particlesByPt(); _nusum = FourMomentum(0., 0., 0., 0.); - for(const Particle& neutrino : neutrinos) - { + 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_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_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; - + 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,179 @@ +// -*- 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/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 + /// @brief Special CMS dressed lepton finder /// /// Find dressed leptons by clustering all leptons and photons - class SpecialDressedLeptons : public FinalState { + class SpecialDressedLeptons_CMS_2016_PAS_TOP_15_006 : public FinalState { public: /// Constructor - SpecialDressedLeptons(const FinalState& fs, const Cut& cut) + SpecialDressedLeptons_CMS_2016_PAS_TOP_15_006(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)); + return unique_ptr(new SpecialDressedLeptons_CMS_2016_PAS_TOP_15_006(*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.pid() == 0)) continue; DressedLepton lepton = DressedLepton(lepCand); for (const Particle& cand : jet.particles()) { + //if (isSame(cand, lepCand)) continue; 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(); } } } 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); + SpecialDressedLeptons_CMS_2016_PAS_TOP_15_006 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"); + const auto& 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_I1610623.cc b/analyses/pluginCMS/CMS_2017_I1610623.cc --- a/analyses/pluginCMS/CMS_2017_I1610623.cc +++ b/analyses/pluginCMS/CMS_2017_I1610623.cc @@ -1,261 +1,253 @@ +// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" - -#include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/WFinder.hh" -#include "Rivet/AnalysisLoader.hh" -#include "Rivet/AnalysisInfo.hh" -#include "Rivet/Tools/RivetYODA.hh" - -#include - namespace Rivet { - /// @brief Add a short analysis description here class CMS_2017_I1610623 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CMS_2017_I1610623); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise and register projections FinalState fs; WFinder wfinder_mu(fs, Cuts::abseta < 2.4 && Cuts::pT > 0*GeV, PID::MUON, 0*GeV, 1000000*GeV, 0*GeV, 0.1, WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::YES, WFinder::MassWindow::MT); //WFinder wfinder_mu(fs, Cuts::abseta < 2.4 && Cuts::pT > 0*GeV, PID::MUON, 0*GeV, 1000000*GeV, 0*GeV, 0.1, WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT); declare(wfinder_mu, "WFinder_mu"); // Define veto FS VetoedFinalState vfs; vfs.addVetoOnThisFinalState(wfinder_mu); vfs.addVetoPairId(PID::MUON); vfs.vetoNeutrinos(); FastJets fastjets(vfs, FastJets::ANTIKT, 0.4); declare(fastjets, "Jets"); //------------- book(_hist_Mult_exc ,"d01-x01-y01"); book(_hist_inc_WJetMult ,"d02-x01-y01"); //------------- book(_hist_JetPt1j ,"d03-x01-y01"); book(_hist_JetPt2j ,"d04-x01-y01"); book(_hist_JetPt3j ,"d05-x01-y01"); book(_hist_JetPt4j ,"d06-x01-y01"); //------------- book(_hist_JetRap1j ,"d07-x01-y01"); book(_hist_JetRap2j ,"d08-x01-y01"); book(_hist_JetRap3j ,"d09-x01-y01"); book(_hist_JetRap4j ,"d10-x01-y01"); //------------- book(_hist_Ht_1j ,"d11-x01-y01"); book(_hist_Ht_2j ,"d12-x01-y01"); book(_hist_Ht_3j ,"d13-x01-y01"); book(_hist_Ht_4j ,"d14-x01-y01"); //------------- book(_hist_dphij1mu_1j , "d15-x01-y01"); book(_hist_dphij2mu_2j , "d16-x01-y01"); book(_hist_dphij3mu_3j , "d17-x01-y01"); book(_hist_dphij4mu_4j , "d18-x01-y01"); //------------- book(_hist_dRmuj_1j , "d19-x01-y01"); } // define function used for filiing inc Njets histo - void Fill(Histo1DPtr& _histJetMult, std::vector& finaljet_list){ + void _fill(Histo1DPtr& _histJetMult, std::vector& finaljet_list){ _histJetMult->fill(0); for (size_t i=0 ; ifill(i+1); // inclusive multiplicity } } /// Perform the per-event analysis void analyze(const Event& event) { /// @todo Do the event by event analysis here const WFinder& wfinder_mu = applyProjection(event, "WFinder_mu"); if (wfinder_mu.bosons().size() != 1) { vetoEvent; } if (wfinder_mu.bosons().size() == 1) { const FourMomentum lepton0 = wfinder_mu.constituentLepton().momentum(); const FourMomentum neutrino = wfinder_mu.constituentNeutrino().momentum(); double WmT = wfinder_mu.mT(); if (WmT < 50.0*GeV) vetoEvent; double pt0 = lepton0.pT(); double eta0 = lepton0.eta(); if ( (fabs(eta0) > 2.4) || (pt0 < 25.0*GeV) ) vetoEvent; // Obtain the jets. vector finaljet_list; vector jet100_list; double HT = 0.0; // loop over jets in an event, pushback in finaljet_list collection for (const Jet& j : applyProjection(event, "Jets").jetsByPt(30.0*GeV)) { const double jrap = j.momentum().rap(); const double jpt = j.momentum().pT(); if ( (fabs(jrap) < 2.4) && (deltaR(lepton0, j.momentum()) > 0.4) ) { if(jpt > 30.0*GeV) { finaljet_list.push_back(j.momentum()); HT += j.momentum().pT(); } if(jpt > 100.0*GeV) { jet100_list.push_back(j.momentum()); } } } // end looping over jets //---------------------- FILL HISTOGRAMS ------------------ // Multiplicity exc plot. _hist_Mult_exc->fill(finaljet_list.size()); // Multiplicity inc plot. - Fill(_hist_inc_WJetMult, finaljet_list); + _fill(_hist_inc_WJetMult, finaljet_list); // dRmuj plot. double mindR(99999); if(jet100_list.size()>=1) { for (unsigned ji = 0; ji < jet100_list.size(); ji++){ double dr_(9999); dr_ = fabs(deltaR(lepton0, jet100_list[ji])); if (dr_ < mindR){ mindR = dr_; } } if(jet100_list[0].pT() > 300.0*GeV){ _hist_dRmuj_1j->fill(mindR); } } if(finaljet_list.size()>=1) { _hist_JetPt1j->fill(finaljet_list[0].pT()); _hist_JetRap1j->fill(fabs(finaljet_list[0].rap())); _hist_Ht_1j->fill(HT); _hist_dphij1mu_1j->fill(deltaPhi(finaljet_list[0].phi(), lepton0.phi())); } if(finaljet_list.size()>=2) { _hist_JetPt2j->fill(finaljet_list[1].pT()); _hist_JetRap2j->fill(fabs(finaljet_list[1].rap())); _hist_Ht_2j->fill(HT); _hist_dphij2mu_2j->fill(deltaPhi(finaljet_list[1].phi(), lepton0.phi())); } if(finaljet_list.size()>=3) { _hist_JetPt3j->fill(finaljet_list[2].pT()); _hist_JetRap3j->fill(fabs(finaljet_list[2].rap())); _hist_Ht_3j->fill(HT); _hist_dphij3mu_3j->fill(deltaPhi(finaljet_list[2].phi(), lepton0.phi())); } if(finaljet_list.size()>=4) { _hist_JetPt4j->fill(finaljet_list[3].pT()); _hist_JetRap4j->fill(fabs(finaljet_list[3].rap())); _hist_Ht_4j->fill(HT); _hist_dphij4mu_4j->fill(deltaPhi(finaljet_list[3].phi(), lepton0.phi())); } } // close the Wboson loop } //void loop /// Normalise histograms etc., after the run void finalize() { const double crossec = !std::isnan(crossSectionPerEvent()) ? crossSection() : 61526.7*picobarn; if (std::isnan(crossSectionPerEvent())){ MSG_INFO("No valid cross-section given, using NNLO xsec calculated by FEWZ " << crossec/picobarn << " pb"); } scale(_hist_Mult_exc, crossec/picobarn/sumOfWeights()); scale(_hist_inc_WJetMult, crossec/picobarn/sumOfWeights()); scale(_hist_JetPt1j, crossec/picobarn/sumOfWeights()); scale(_hist_JetPt2j, crossec/picobarn/sumOfWeights()); scale(_hist_JetPt3j, crossec/picobarn/sumOfWeights()); scale(_hist_JetPt4j, crossec/picobarn/sumOfWeights()); scale(_hist_JetRap1j, crossec/picobarn/sumOfWeights()); scale(_hist_JetRap2j, crossec/picobarn/sumOfWeights()); scale(_hist_JetRap3j, crossec/picobarn/sumOfWeights()); scale(_hist_JetRap4j, crossec/picobarn/sumOfWeights()); scale(_hist_Ht_1j, crossec/picobarn/sumOfWeights()); scale(_hist_Ht_2j, crossec/picobarn/sumOfWeights()); scale(_hist_Ht_3j, crossec/picobarn/sumOfWeights()); scale(_hist_Ht_4j, crossec/picobarn/sumOfWeights()); scale(_hist_dphij1mu_1j, crossec/picobarn/sumOfWeights()); scale(_hist_dphij2mu_2j, crossec/picobarn/sumOfWeights()); scale(_hist_dphij3mu_3j, crossec/picobarn/sumOfWeights()); scale(_hist_dphij4mu_4j, crossec/picobarn/sumOfWeights()); scale(_hist_dRmuj_1j, crossec/picobarn/sumOfWeights()); } //@} private: /// @name Histograms //@{ Histo1DPtr _hist_Mult_exc; Histo1DPtr _hist_inc_WJetMult; Histo1DPtr _hist_JetPt1j; Histo1DPtr _hist_JetPt2j; Histo1DPtr _hist_JetPt3j; Histo1DPtr _hist_JetPt4j; Histo1DPtr _hist_JetRap1j; Histo1DPtr _hist_JetRap2j; Histo1DPtr _hist_JetRap3j; Histo1DPtr _hist_JetRap4j; Histo1DPtr _hist_Ht_1j; Histo1DPtr _hist_Ht_2j; Histo1DPtr _hist_Ht_3j; Histo1DPtr _hist_Ht_4j; Histo1DPtr _hist_dphij1mu_1j; Histo1DPtr _hist_dphij2mu_2j; Histo1DPtr _hist_dphij3mu_3j; Histo1DPtr _hist_dphij4mu_4j; Histo1DPtr _hist_dRmuj_1j; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CMS_2017_I1610623); } diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh --- a/include/Rivet/Analysis.hh +++ b/include/Rivet/Analysis.hh @@ -1,1335 +1,1335 @@ // -*- C++ -*- #ifndef RIVET_Analysis_HH #define RIVET_Analysis_HH #include "Rivet/Config/RivetCommon.hh" #include "Rivet/AnalysisInfo.hh" #include "Rivet/Event.hh" #include "Rivet/Projection.hh" #include "Rivet/ProjectionApplier.hh" #include "Rivet/ProjectionHandler.hh" #include "Rivet/AnalysisLoader.hh" #include "Rivet/Tools/Cuts.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleUtils.hh" #include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/Tools/RivetMT2.hh" #include "Rivet/Tools/RivetYODA.hh" #include "Rivet/Tools/Percentile.hh" #include "Rivet/Projections/CentralityProjection.hh" #include /// @def vetoEvent /// Preprocessor define for vetoing events, including the log message and return. #define vetoEvent \ do { MSG_DEBUG("Vetoing event on line " << __LINE__ << " of " << __FILE__); return; } while(0) namespace Rivet { // Convenience for analysis writers using std::cout; using std::cerr; using std::endl; using std::tuple; using std::stringstream; using std::swap; using std::numeric_limits; // Forward declaration class AnalysisHandler; /// @brief This is the base class of all analysis classes in Rivet. /// /// There are /// three virtual functions which should be implemented in base classes: /// /// void init() is called by Rivet before a run is started. Here the /// analysis class should book necessary histograms. The needed /// projections should probably rather be constructed in the /// constructor. /// /// void analyze(const Event&) is called once for each event. Here the /// analysis class should apply the necessary Projections and fill the /// histograms. /// /// void finalize() is called after a run is finished. Here the analysis /// class should do whatever manipulations are necessary on the /// histograms. Writing the histograms to a file is, however, done by /// the Rivet class. class Analysis : public ProjectionApplier { /// The AnalysisHandler is a friend. friend class AnalysisHandler; public: /// @name Standard constructors and destructors. //@{ // /// The default constructor. // Analysis(); /// Constructor Analysis(const std::string& name); /// The destructor. virtual ~Analysis() {} //@} public: /// @name Main analysis methods //@{ /// Initialize this analysis object. A concrete class should here /// book all necessary histograms. An overridden function must make /// sure it first calls the base class function. virtual void init() { } /// Analyze one event. A concrete class should here apply the /// necessary projections on the \a event and fill the relevant /// histograms. An overridden function must make sure it first calls /// the base class function. virtual void analyze(const Event& event) = 0; /// Finalize this analysis object. A concrete class should here make /// all necessary operations on the histograms. Writing the /// histograms to a file is, however, done by the Rivet class. An /// overridden function must make sure it first calls the base class /// function. virtual void finalize() { } //@} public: /// @name Metadata /// Metadata is used for querying from the command line and also for /// building web pages and the analysis pages in the Rivet manual. //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored. const AnalysisInfo& info() const { assert(_info && "No AnalysisInfo object :O"); return *_info; } /// @brief Get the name of the analysis. /// /// By default this is computed by combining the results of the /// experiment, year and Spires ID metadata methods and you should /// only override it if there's a good reason why those won't /// work. If options has been set for this instance, a /// corresponding string is appended at the end. virtual std::string name() const { return ( (info().name().empty()) ? _defaultname : info().name() ) + _optstring; } /// Get name of reference data file, which could be different from plugin name virtual std::string getRefDataName() const { return (info().getRefDataName().empty()) ? _defaultname : info().getRefDataName(); } /// Set name of reference data file, which could be different from plugin name virtual void setRefDataName(const std::string& ref_data="") { info().setRefDataName(!ref_data.empty() ? ref_data : name()); } /// Get the Inspire ID code for this analysis. virtual std::string inspireId() const { return info().inspireId(); } /// Get the SPIRES ID code for this analysis (~deprecated). virtual std::string spiresId() const { return info().spiresId(); } /// @brief Names & emails of paper/analysis authors. /// /// Names and email of authors in 'NAME \' format. The first /// name in the list should be the primary contact person. virtual std::vector authors() const { return info().authors(); } /// @brief Get a short description of the analysis. /// /// Short (one sentence) description used as an index entry. /// Use @a description() to provide full descriptive paragraphs /// of analysis details. virtual std::string summary() const { return info().summary(); } /// @brief Get a full description of the analysis. /// /// Full textual description of this analysis, what it is useful for, /// what experimental techniques are applied, etc. Should be treated /// as a chunk of restructuredText (http://docutils.sourceforge.net/rst.html), /// with equations to be rendered as LaTeX with amsmath operators. virtual std::string description() const { return info().description(); } /// @brief Information about the events needed as input for this analysis. /// /// Event types, energies, kinematic cuts, particles to be considered /// stable, etc. etc. Should be treated as a restructuredText bullet list /// (http://docutils.sourceforge.net/rst.html) virtual std::string runInfo() const { return info().runInfo(); } /// Experiment which performed and published this analysis. virtual std::string experiment() const { return info().experiment(); } /// Collider on which the experiment ran. virtual std::string collider() const { return info().collider(); } /// When the original experimental analysis was published. virtual std::string year() const { return info().year(); } /// The luminosity in inverse femtobarn virtual std::string luminosityfb() const { return info().luminosityfb(); } /// Journal, and preprint references. virtual std::vector references() const { return info().references(); } /// BibTeX citation key for this article. virtual std::string bibKey() const { return info().bibKey(); } /// BibTeX citation entry for this article. virtual std::string bibTeX() const { return info().bibTeX(); } /// Whether this analysis is trusted (in any way!) virtual std::string status() const { return (info().status().empty()) ? "UNVALIDATED" : info().status(); } /// Any work to be done on this analysis. virtual std::vector todos() const { return info().todos(); } /// Return the allowed pairs of incoming beams required by this analysis. virtual const std::vector& requiredBeams() const { return info().beams(); } /// Declare the allowed pairs of incoming beams required by this analysis. virtual Analysis& setRequiredBeams(const std::vector& requiredBeams) { info().setBeams(requiredBeams); return *this; } /// Sets of valid beam energy pairs, in GeV virtual const std::vector >& requiredEnergies() const { return info().energies(); } /// Get vector of analysis keywords virtual const std::vector & keywords() const { return info().keywords(); } /// Declare the list of valid beam energy pairs, in GeV virtual Analysis& setRequiredEnergies(const std::vector >& requiredEnergies) { info().setEnergies(requiredEnergies); return *this; } //@} /// @name Internal metadata modifying methods //@{ /// Get the actual AnalysisInfo object in which all this metadata is stored (non-const). AnalysisInfo& info() { assert(_info && "No AnalysisInfo object :O"); return *_info; } //@} /// @name Run conditions //@{ /// Incoming beams for this run const ParticlePair& beams() const; /// Incoming beam IDs for this run const PdgIdPair beamIds() const; /// Centre of mass energy for this run double sqrtS() const; /// Check if we are running rivet-merge. bool merging() const { return sqrtS() <= 0.0; } //@} /// @name Analysis / beam compatibility testing /// @todo Replace with beamsCompatible() with no args (calling beams() function internally) /// @todo Add beamsMatch() methods with same (shared-code?) tolerance as in beamsCompatible() //@{ /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const ParticlePair& beams) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(PdgId beam1, PdgId beam2, double e1, double e2) const; /// Check if analysis is compatible with the provided beam particle IDs and energies bool isCompatible(const PdgIdPair& beams, const std::pair& energies) const; //@} /// Access the controlling AnalysisHandler object. AnalysisHandler& handler() const { return *_analysishandler; } protected: /// Get a Log object based on the name() property of the calling analysis object. Log& getLog() const; /// Get the process cross-section in pb. Throws if this hasn't been set. double crossSection() const; /// Get the process cross-section per generated event in pb. Throws if this /// hasn't been set. double crossSectionPerEvent() const; /// @brief Get the number of events seen (via the analysis handler). /// /// @note Use in the finalize phase only. size_t numEvents() const; /// @brief Get the sum of event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW() const; /// Alias double sumOfWeights() const { return sumW(); } /// @brief Get the sum of squared event weights seen (via the analysis handler). /// /// @note Use in the finalize phase only. double sumW2() const; protected: /// @name Histogram paths //@{ /// Get the canonical histogram "directory" path for this analysis. const std::string histoDir() const; /// Get the canonical histogram path for the named histogram in this analysis. const std::string histoPath(const std::string& hname) const; /// Get the canonical histogram path for the numbered histogram in this analysis. const std::string histoPath(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; /// Get the internal histogram name for given d, x and y (cf. HepData) const std::string mkAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const; //@} /// @name Histogram reference data //@{ /// Get reference data for a named histo /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(const string& hname) const { _cacheRefData(); MSG_TRACE("Using histo bin edges for " << name() << ":" << hname); if (!_refdata[hname]) { MSG_ERROR("Can't find reference histogram " << hname); throw Exception("Reference data " + hname + " not found."); } return dynamic_cast(*_refdata[hname]); } /// Get reference data for a numbered histo /// @todo SFINAE to ensure that the type inherits from YODA::AnalysisObject? template const T& refData(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { const string hname = mkAxisCode(datasetId, xAxisId, yAxisId); return refData(hname); } //@} /// @name Counter booking //@{ /// Book a counter. CounterPtr & book(CounterPtr &, const std::string& name, const std::string& title=""); // const std::string& valtitle="" /// Book a counter, using a path generated from the dataset and axis ID codes /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. CounterPtr & book(CounterPtr &, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title=""); // const std::string& valtitle="" //@} /// @name 1D histogram booking //@{ /// Book a 1D histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Histo1DPtr & book(Histo1DPtr &,const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr & book(Histo1DPtr &,const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with non-uniform bins defined by the vector of bin edges @a binedges . Histo1DPtr & book(Histo1DPtr &,const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram with binning from a reference scatter. Histo1DPtr & book(Histo1DPtr &,const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. Histo1DPtr & book(Histo1DPtr &,const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo1DPtr & book(Histo1DPtr &,unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D histogram booking //@{ /// Book a 2D histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a /// ylower - @a yupper respectively along the x- and y-axis. Histo2DPtr & book(Histo2DPtr &,const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr & book(Histo2DPtr &,const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with non-uniform bins defined by the /// vectors of bin edges @a xbinedges and @a ybinedges. Histo2DPtr & book(Histo2DPtr &,const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram with binning from a reference scatter. Histo2DPtr & book(Histo2DPtr &,const std::string& name, const Scatter3D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram, using the binnings in the reference data histogram. Histo2DPtr & book(Histo2DPtr &,const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Histo2DPtr & book(Histo2DPtr &,unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); //@} /// @name 1D profile histogram booking //@{ /// Book a 1D profile histogram with @a nbins uniformly distributed across the range @a lower - @a upper . Profile1DPtr & book(Profile1DPtr &, const std::string& name, size_t nbins, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr & book(Profile1DPtr &, const std::string& name, const std::vector& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with non-uniform bins defined by the vector of bin edges @a binedges . Profile1DPtr & book(Profile1DPtr &, const std::string& name, const std::initializer_list& binedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram with binning from a reference scatter. Profile1DPtr & book(Profile1DPtr &, const std::string& name, const Scatter2D& refscatter, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. Profile1DPtr & book(Profile1DPtr &, const std::string& name, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// Book a 1D profile histogram, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. Profile1DPtr & book(Profile1DPtr &, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} /// @name 2D profile histogram booking //@{ /// Book a 2D profile histogram with @a nxbins and @a nybins uniformly /// distributed across the ranges @a xlower - @a xupper and @a ylower - @a /// yupper respectively along the x- and y-axis. Profile2DPtr & book(Profile2DPtr &, const std::string& name, size_t nxbins, double xlower, double xupper, size_t nybins, double ylower, double yupper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr & book(Profile2DPtr &, const std::string& name, const std::vector& xbinedges, const std::vector& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// Book a 2D profile histogram with non-uniform bins defined by the vectorx /// of bin edges @a xbinedges and @a ybinedges. Profile2DPtr & book(Profile2DPtr &, const std::string& name, const std::initializer_list& xbinedges, const std::initializer_list& ybinedges, const std::string& title="", const std::string& xtitle="", const std::string& ytitle="", const std::string& ztitle=""); /// @todo REINSTATE // /// Book a 2D profile histogram with binning from a reference scatter. // Profile2DPtr & book(const Profile2DPtr &, const std::string& name, // const Scatter3D& refscatter, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); // /// Book a 2D profile histogram, using the binnings in the reference data histogram. // Profile2DPtr & book(const Profile2DPtr &, const std::string& name, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); // /// Book a 2D profile histogram, using the binnings in the reference data histogram. // /// // /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. // Profile2DPtr & book(const Profile2DPtr &, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, // const std::string& title="", // const std::string& xtitle="", // const std::string& ytitle="", // const std::string& ztitle=""); //@} /// @name 2D scatter booking //@{ /// @brief Book a 2-dimensional data point set with the given name. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed... /// assuming that there is a reference histo with the same name: if there /// isn't, an exception will be thrown. Scatter2DPtr & book(Scatter2DPtr & s2d, const string& hname, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set, using the binnings in the reference data histogram. /// /// The paper, dataset and x/y-axis IDs will be used to build the histo name in the HepData standard way. /// /// @note Unlike histogram booking, scatter booking by default makes no /// attempt to use reference data to pre-fill the data object. If you want /// this, which is sometimes useful e.g. when the x-position is not really /// meaningful and can't be extracted from the data, then set the @a /// copy_pts parameter to true. This creates points to match the reference /// data's x values and errors, but with the y values and errors zeroed. Scatter2DPtr & book(Scatter2DPtr & s2d, unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId, bool copy_pts=false, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set with equally spaced x-points in a range. /// /// The y values and errors will be set to 0. Scatter2DPtr & book(Scatter2DPtr & s2d, const string& hname, size_t npts, double lower, double upper, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); /// @brief Book a 2-dimensional data point set based on provided contiguous "bin edges". /// /// The y values and errors will be set to 0. Scatter2DPtr & book(Scatter2DPtr & s2d, const string& hname, const std::vector& binedges, const std::string& title, const std::string& xtitle, const std::string& ytitle); /// Book a 2-dimensional data point set with x-points from an existing scatter and a new path. Scatter2DPtr book(Scatter2DPtr & s2d, const Scatter2DPtr & scPtr, const std::string& path, const std::string& title="", const std::string& xtitle="", const std::string& ytitle=""); //@} public: /// @name Accessing options for this Analysis instance. //@{ /// Get an option for this analysis instance as a string. std::string getOption(std::string optname) { if ( _options.find(optname) != _options.end() ) return _options.find(optname)->second; return ""; } /// Get an option for this analysis instance converted to a /// specific type (given by the specified @a def value). template T getOption(std::string optname, T def) { if (_options.find(optname) == _options.end()) return def; std::stringstream ss; ss << _options.find(optname)->second; T ret; ss >> ret; return ret; } //@} /// @brief Book a CentralityProjection /// /// Using a SingleValueProjection, @a proj, giving the value of an /// experimental observable to be used as a centrality estimator, /// book a CentralityProjection based on the experimentally /// measured pecentiles of this observable (as given by the /// reference data for the @a calHistName histogram in the @a /// calAnaName analysis. If a preloaded file with the output of a /// run using the @a calAnaName analysis contains a valid /// generated @a calHistName histogram, it will be used as an /// optional percentile binning. Also if this preloaded file /// contains a histogram with the name @a calHistName with an /// appended "_IMP" This histogram will be used to add an optional /// centrality percentile based on the generated impact /// parameter. If @increasing is true, a low (high) value of @proj /// is assumed to correspond to a more peripheral (central) event. const CentralityProjection& declareCentrality(const SingleValueProjection &proj, string calAnaName, string calHistName, const string projName, bool increasing = false); /// @brief Book a Percentile wrapper around AnalysisObjects. /// /// Based on a previously registered CentralityProjection named @a /// projName book one AnalysisObject for each @a centralityBin and /// name them according to the corresponding code in the @a ref /// vector. /// /// @todo Convert to just be called book() cf. others template Percentile bookPercentile(string projName, vector > centralityBins, vector > ref) { typedef typename ReferenceTraits::RefT RefT; typedef rivet_shared_ptr> WrapT; Percentile pctl(this, projName); const int nCent = centralityBins.size(); for (int iCent = 0; iCent < nCent; ++iCent) { const string axisCode = mkAxisCode(std::get<0>(ref[iCent]), std::get<1>(ref[iCent]), std::get<2>(ref[iCent])); const RefT & refscatter = refData(axisCode); WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode))); wtf = addAnalysisObject(wtf); CounterPtr cnt(_weightNames(), Counter(histoPath("TMP/COUNTER/" + axisCode))); cnt = addAnalysisObject(cnt); pctl.add(wtf, cnt, centralityBins[iCent]); } return pctl; } // /// @brief Book Percentile wrappers around AnalysisObjects. // /// // /// Based on a previously registered CentralityProjection named @a // /// projName book one (or several) AnalysisObject(s) named // /// according to @a ref where the x-axis will be filled according // /// to the percentile output(s) of the @projName. // /// // /// @todo Convert to just be called book() cf. others // template // PercentileXaxis bookPercentileXaxis(string projName, // tuple ref) { // typedef typename ReferenceTraits::RefT RefT; // typedef rivet_shared_ptr> WrapT; // PercentileXaxis pctl(this, projName); // const string axisCode = mkAxisCode(std::get<0>(ref), // std::get<1>(ref), // std::get<2>(ref)); // const RefT & refscatter = refData(axisCode); // WrapT wtf(_weightNames(), T(refscatter, histoPath(axisCode))); // wtf = addAnalysisObject(wtf); // CounterPtr cnt(_weightNames(), Counter()); // cnt = addAnalysisObject(cnt); // pctl.add(wtf, cnt); // return pctl; // } private: // Functions that have to be defined in the .cc file to avoid circular #includes /// Get the list of weight names from the handler vector _weightNames() const; /// Get the default/nominal weight index size_t _defaultWeightIndex() const; /// Get an AO from another analysis MultiweightAOPtr _getOtherAnalysisObject(const std::string & ananame, const std::string& name); /// Check that analysis objects aren't being booked/registered outside the init stage void _checkBookInit() const; private: /// To be used in finalize context only: class CounterAdapter { public: CounterAdapter(double x) : x_(x ) {} CounterAdapter(const YODA::Counter & c) : x_(c.val() ) {} // CounterAdapter(CounterPtr cp) : x_(cp->val() ) {} CounterAdapter(const YODA::Scatter1D & s) : x_(s.points()[0].x()) { assert( s.numPoints() == 1 || "Can only scale by a single value."); } // CounterAdapter(Scatter1DPtr sp) : x_(sp->points()[0].x()) { // assert( sp->numPoints() == 1 || "Can only scale by a single value."); // } operator double() const { return x_; } private: double x_; }; public: double dbl(double x) { return x; } double dbl(const YODA::Counter & c) { return c.val(); } double dbl(const YODA::Scatter1D & s) { assert( s.numPoints() == 1 ); return s.points()[0].x(); } /// @name Analysis object manipulation /// @todo Should really be protected: only public to keep BinnedHistogram happy for now... //@{ /// Multiplicatively scale the given counter, @a cnt, by factor @s factor. void scale(CounterPtr cnt, CounterAdapter factor); /// Multiplicatively scale the given counters, @a cnts, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of CounterPtrs void scale(const std::vector& cnts, CounterAdapter factor) { for (auto& c : cnts) scale(c, factor); } /// @todo YUCK! template void scale(const CounterPtr (&cnts)[array_size], CounterAdapter factor) { // for (size_t i = 0; i < std::extent::value; ++i) scale(cnts[i], factor); for (auto& c : cnts) scale(c, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo1DPtr histo, CounterAdapter norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void normalize(const std::vector& histos, CounterAdapter norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo1DPtr (&histos)[array_size], CounterAdapter norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo1DPtr histo, CounterAdapter factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo1DPtrs void scale(const std::vector& histos, CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo1DPtr (&histos)[array_size], CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// Normalize the given histogram, @a histo, to area = @a norm. void normalize(Histo2DPtr histo, CounterAdapter norm=1.0, bool includeoverflows=true); /// Normalize the given histograms, @a histos, to area = @a norm. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void normalize(const std::vector& histos, CounterAdapter norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// @todo YUCK! template void normalize(const Histo2DPtr (&histos)[array_size], CounterAdapter norm=1.0, bool includeoverflows=true) { for (auto& h : histos) normalize(h, norm, includeoverflows); } /// Multiplicatively scale the given histogram, @a histo, by factor @s factor. void scale(Histo2DPtr histo, CounterAdapter factor); /// Multiplicatively scale the given histograms, @a histos, by factor @s factor. /// @note Constness intentional, if weird, to allow passing rvalue refs of smart ptrs (argh) /// @todo Use SFINAE for a generic iterable of Histo2DPtrs void scale(const std::vector& histos, CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// @todo YUCK! template void scale(const Histo2DPtr (&histos)[array_size], CounterAdapter factor) { for (auto& h : histos) scale(h, factor); } /// Helper for counter division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(CounterPtr c1, CounterPtr c2, Scatter1DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Counter& c1, const YODA::Counter& c2, Scatter1DPtr s) const; /// Helper for histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile1DPtr p1, Profile1DPtr p2, Scatter2DPtr s) const; /// Helper for profile histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile1D& p1, const YODA::Profile1D& p2, Scatter2DPtr s) const; /// Helper for 2D histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Histo2DPtr h1, Histo2DPtr h2, Scatter3DPtr s) const; /// Helper for 2D histogram division with raw YODA objects. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Histo2D& h1, const YODA::Histo2D& h2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(Profile2DPtr p1, Profile2DPtr p2, Scatter3DPtr s) const; /// Helper for 2D profile histogram division with raw YODA objects /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void divide(const YODA::Profile2D& p1, const YODA::Profile2D& p2, Scatter3DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram efficiency calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void efficiency(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(Histo1DPtr h1, Histo1DPtr h2, Scatter2DPtr s) const; /// Helper for histogram asymmetry calculation. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void asymm(const YODA::Histo1D& h1, const YODA::Histo1D& h2, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(Histo1DPtr h, Scatter2DPtr s) const; /// Helper for converting a differential histo to an integral one. /// /// @note Assigns to the (already registered) output scatter, @a s. Preserves the path information of the target. void integrate(const Histo1D& h, Scatter2DPtr s) const; //@} public: /// List of registered analysis data objects const vector& analysisObjects() const { return _analysisobjects; } protected: /// @name Data object registration, retrieval, and removal //@{ /// Register a data object in the histogram system template AO addAnalysisObject(const AO & aonew) { _checkBookInit(); for (const MultiweightAOPtr& ao : analysisObjects()) { // Check AO base-name first ao.get()->setActiveWeightIdx(_defaultWeightIndex()); aonew.get()->setActiveWeightIdx(_defaultWeightIndex()); if (ao->path() != aonew->path()) continue; // If base-name matches, check compatibility // NB. This evil is because dynamic_ptr_cast can't work on rivet_shared_ptr directly AO aoold = AO(dynamic_pointer_cast(ao.get())); //< OMG if ( !aoold || !bookingCompatible(aonew, aoold) ) { MSG_WARNING("Found incompatible pre-existing data object with same base path " << aonew->path() << " for " << name()); throw LookupError("Found incompatible pre-existing data object with same base path during AO booking"); } // Finally, check all weight variations for (size_t weightIdx = 0; weightIdx < _weightNames().size(); ++weightIdx) { aoold.get()->setActiveWeightIdx(weightIdx); aonew.get()->setActiveWeightIdx(weightIdx); if (aoold->path() != aonew->path()) { MSG_WARNING("Found incompatible pre-existing data object with different weight-path " << aonew->path() << " for " << name()); throw LookupError("Found incompatible pre-existing data object with same weight-path during AO booking"); } } // They're fully compatible: bind and return aoold.get()->unsetActiveWeight(); MSG_TRACE("Bound pre-existing data object " << aoold->path() << " for " << name()); return aoold; } // No equivalent found MSG_TRACE("Registered " << aonew->annotation("Type") << " " << aonew->path() << " for " << name()); aonew.get()->unsetActiveWeight(); _analysisobjects.push_back(aonew); return aonew; } /// Unregister a data object from the histogram system (by name) void removeAnalysisObject(const std::string& path); /// Unregister a data object from the histogram system (by pointer) void removeAnalysisObject(const MultiweightAOPtr & ao); // /// Get all data objects, for all analyses, from the AnalysisHandler // /// @todo Can we remove this? Why not call handler().getData()? // vector getAllData(bool includeorphans) const; /// Get a Rivet data object from the histogram system template const AO getAnalysisObject(const std::string& aoname) const { for (const MultiweightAOPtr& ao : analysisObjects()) { ao.get()->setActiveWeightIdx(_defaultWeightIndex()); if (ao->path() == histoPath(aoname)) { // return dynamic_pointer_cast(ao); return AO(dynamic_pointer_cast(ao.get())); } } throw LookupError("Data object " + histoPath(aoname) + " not found"); } // /// Get a data object from the histogram system // template // const std::shared_ptr getAnalysisObject(const std::string& name) const { // foreach (const AnalysisObjectPtr& ao, analysisObjects()) { // if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); // } // throw LookupError("Data object " + histoPath(name) + " not found"); // } // /// Get a data object from the histogram system (non-const) // template // std::shared_ptr getAnalysisObject(const std::string& name) { // foreach (const AnalysisObjectPtr& ao, analysisObjects()) { // if (ao->path() == histoPath(name)) return dynamic_pointer_cast(ao); // } // throw LookupError("Data object " + histoPath(name) + " not found"); // } /// Get a data object from another analysis (e.g. preloaded /// calibration histogram). template AO getAnalysisObject(const std::string& ananame, const std::string& aoname) { MultiweightAOPtr ao = _getOtherAnalysisObject(ananame, aoname); // return dynamic_pointer_cast(ao); return AO(dynamic_pointer_cast(ao.get())); } // /// Get a named Histo1D object from the histogram system // const Histo1DPtr getHisto1D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Histo1D object from the histogram system (non-const) // Histo1DPtr getHisto1D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Histo1D object from the histogram system by axis ID codes (non-const) // const Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Histo1D object from the histogram system by axis ID codes (non-const) // Histo1DPtr getHisto1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a named Histo2D object from the histogram system // const Histo2DPtr getHisto2D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Histo2D object from the histogram system (non-const) // Histo2DPtr getHisto2D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Histo2D object from the histogram system by axis ID codes (non-const) // const Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Histo2D object from the histogram system by axis ID codes (non-const) // Histo2DPtr getHisto2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a named Profile1D object from the histogram system // const Profile1DPtr getProfile1D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Profile1D object from the histogram system (non-const) // Profile1DPtr getProfile1D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Profile1D object from the histogram system by axis ID codes (non-const) // const Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Profile1D object from the histogram system by axis ID codes (non-const) // Profile1DPtr getProfile1D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a named Profile2D object from the histogram system // const Profile2DPtr getProfile2D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Profile2D object from the histogram system (non-const) // Profile2DPtr getProfile2D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Profile2D object from the histogram system by axis ID codes (non-const) // const Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Profile2D object from the histogram system by axis ID codes (non-const) // Profile2DPtr getProfile2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a named Scatter2D object from the histogram system // const Scatter2DPtr getScatter2D(const std::string& name) const { // return getAnalysisObject(name); // } // /// Get a named Scatter2D object from the histogram system (non-const) // Scatter2DPtr getScatter2D(const std::string& name) { // return getAnalysisObject(name); // } // /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) // const Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } // /// Get a Scatter2D object from the histogram system by axis ID codes (non-const) // Scatter2DPtr getScatter2D(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) { // return getAnalysisObject(makeAxisCode(datasetId, xAxisId, yAxisId)); // } //@} private: /// Name passed to constructor (used to find .info analysis data file, and as a fallback) string _defaultname; /// Pointer to analysis metadata object unique_ptr _info; /// Storage of all plot objects /// @todo Make this a map for fast lookup by path? vector _analysisobjects; /// @name Cross-section variables //@{ double _crossSection; bool _gotCrossSection; //@} /// The controlling AnalysisHandler object. AnalysisHandler* _analysishandler; /// Collection of cached refdata to speed up many autobookings: the /// reference data file should only be read once. mutable std::map _refdata; /// Options the (this instance of) the analysis map _options; /// The string of options. string _optstring; private: /// @name Utility functions //@{ /// Get the reference data for this paper and cache it. void _cacheRefData() const; //@} /// The assignment operator is private and must never be called. /// In fact, it should not even be implemented. Analysis& operator=(const Analysis&); }; } // Include definition of analysis plugin system so that analyses automatically see it when including Analysis.hh #include "Rivet/AnalysisBuilder.hh" /// @def DECLARE_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism. -#define DECLARE_RIVET_PLUGIN(clsname) Rivet::AnalysisBuilder plugin_ ## clsname +#define DECLARE_RIVET_PLUGIN(clsname) ::Rivet::AnalysisBuilder plugin_ ## clsname /// @def DECLARE_ALIASED_RIVET_PLUGIN /// Preprocessor define to prettify the global-object plugin hook mechanism, with an extra alias name for this analysis. // #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) Rivet::AnalysisBuilder plugin_ ## clsname ## ( ## #alias ## ) #define DECLARE_ALIASED_RIVET_PLUGIN(clsname, alias) DECLARE_RIVET_PLUGIN(clsname)( #alias ) /// @def DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR /// Preprocessor define to prettify the manky constructor with name string argument #define DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) clsname() : Analysis(# clsname) {} /// @def DEFAULT_RIVET_ANALYSIS_CTOR /// Slight abbreviation for DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR #define DEFAULT_RIVET_ANALYSIS_CTOR(clsname) DEFAULT_RIVET_ANALYSIS_CONSTRUCTOR(clsname) #endif