diff --git a/analyses/pluginALICE/ALICE_2016_I1507157.cc b/analyses/pluginALICE/ALICE_2016_I1507157.cc --- a/analyses/pluginALICE/ALICE_2016_I1507157.cc +++ b/analyses/pluginALICE/ALICE_2016_I1507157.cc @@ -1,201 +1,186 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/AliceCommon.hh" #include "Rivet/Projections/PrimaryParticles.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/EventMixingFinalState.hh" namespace Rivet { /// @brief Correlations of identified particles in pp. /// Also showcasing use of EventMixingFinalState. + class ALICE_2016_I1507157 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ALICE_2016_I1507157); /// @name Analysis methods //@{ - + /// @brief Calculate angular distance between particles. double phaseDif(double a1, double a2){ double dif = a1 - a2; while (dif < -M_PI/2) dif += 2*M_PI; while (dif > 3*M_PI/2) dif -= 2*M_PI; return dif; } + /// Book histograms and initialise projections before the run void init() { double etamax = 0.8; - double pTmin = 0.5; // GeV - + double pTmin = 0.5; // GeV + // Trigger declare(ALICE::V0AndTrigger(), "V0-AND"); // Charged tracks used to manage the mixing observable. ChargedFinalState cfsMult(Cuts::abseta < etamax); declare(cfsMult, "CFSMult"); - + // Primary particles. - PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS, - Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON, + PrimaryParticles pp({Rivet::PID::PIPLUS, Rivet::PID::KPLUS, + Rivet::PID::K0S, Rivet::PID::K0L, Rivet::PID::PROTON, Rivet::PID::NEUTRON, Rivet::PID::LAMBDA, Rivet::PID::SIGMAMINUS, - Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0, + Rivet::PID::SIGMAPLUS, Rivet::PID::XIMINUS, Rivet::PID::XI0, Rivet::PID::OMEGAMINUS},Cuts::abseta < etamax && Cuts::pT > pTmin*GeV); declare(pp,"APRIM"); // The event mixing projection - declare(EventMixingFinalState(cfsMult, pp, 5, 0, 100, 10),"EVM"); + declare(EventMixingFinalState(&cfsMult, pp, 5, 0, 100, 10),"EVM"); // The particle pairs. pid = {{211, -211}, {321, -321}, {2212, -2212}, {3122, -3122}, {211, 211}, {321, 321}, {2212, 2212}, {3122, 3122}, {2212, 3122}, {2212, -3122}}; // The associated histograms in the data file. vector refdata = {"d04-x01-y01","d04-x01-y02","d04-x01-y03", "d06-x01-y02","d05-x01-y01","d05-x01-y02","d05-x01-y03","d06-x01-y01", "d01-x01-y02","d02-x01-y02"}; ratio.resize(refdata.size()); signal.resize(refdata.size()); background.resize(refdata.size()); for (int i = 0, N = refdata.size(); i < N; ++i) { // The ratio plots. - book(ratio[i], refdata[i], true); - // Signal and mixed background. - book(signal[i], "/TMP/" + refdata[i] + "-s", *ratio[i], refdata[i] + "-s"); - book(background[i], "/TMP/" + refdata[i] + "-b", *ratio[i], refdata[i] + "-b"); + book(ratio[i],refdata[i], true); + // Signal and mixed background. + book(signal[i],"/TMP/" + refdata[i] + "-s", *ratio[i], refdata[i] + "-s"); + book(background[i],"/TMP/" + refdata[i] + "-b", *ratio[i], refdata[i] + "-b"); // Number of signal and mixed pairs. - nsp.push_back(0.); + nsp.push_back(0.); nmp.push_back(0.); - } + } } /// Perform the per-event analysis void analyze(const Event& event) { - // Triggering if (!apply(event, "V0-AND")()) return; // The projections - const PrimaryParticles& pp = applyProjection(event,"APRIM"); - const EventMixingFinalState& evm = + const PrimaryParticles& pp = + applyProjection(event,"APRIM"); + const EventMixingFinalState& evm = applyProjection(event, "EVM"); - // Get all mixing events - vector mixEvents = evm.getMixingEvents(); - - // If there are not enough events to mix, don't fill any histograms - if(mixEvents.size() == 0) - return; - - // Make a vector of mixed event particles - vector mixParticles; - size_t pSize = 0; - for(size_t i = 0; i < mixEvents.size(); ++i) - pSize+=mixEvents[i].size(); - mixParticles.reserve(pSize); - for(size_t i = 0; i < mixEvents.size(); ++i) - mixParticles.insert(mixParticles.end(), mixEvents[i].begin(), - mixEvents[i].end()); - - // Shuffle the particles in the mixing events - random_shuffle(mixParticles.begin(), mixParticles.end()); + // Test if we have enough mixing events available to continue. + if (!evm.hasMixingEvents()) return; for(const Particle& p1 : pp.particles()) { // Start by doing the signal distributions for(const Particle& p2 : pp.particles()) { if(p1 == p2) continue; double dEta = abs(p1.eta() - p2.eta()); double dPhi = phaseDif(p1.phi(), p2.phi()); if(dEta < 1.3) { for (int i = 0, N = pid.size(); i < N; ++i) { int pid1 = pid[i].first; int pid2 = pid[i].second; bool samesign = (pid1 * pid2 > 0); - if (samesign && ((pid1 == p1.pid() && pid2 == p2.pid()) || + if (samesign && ((pid1 == p1.pid() && pid2 == p2.pid()) || (pid1 == -p1.pid() && pid2 == -p2.pid()))) { signal[i]->fill(dPhi); nsp[i] += 1.0; } if (!samesign && abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == p2.pid()) { signal[i]->fill(dPhi); nsp[i] += 1.0; } - if (!samesign && abs(pid1) != abs(pid2) && + if (!samesign && abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == p2.pid()) || (pid2 == p1.pid() && pid1 == p2.pid()) ) ) { signal[i]->fill(dPhi); nsp[i] += 1.0; } } } } // Then do the background distribution - for(const Particle& pMix : mixParticles){ + for(const Particle& pMix : evm.particles()){ double dEta = abs(p1.eta() - pMix.eta()); double dPhi = phaseDif(p1.phi(), pMix.phi()); if(dEta < 1.3) { for (int i = 0, N = pid.size(); i < N; ++i) { int pid1 = pid[i].first; int pid2 = pid[i].second; bool samesign = (pid1 * pid2 > 0); - if (samesign && ((pid1 == p1.pid() && pid2 == pMix.pid()) || + if (samesign && ((pid1 == p1.pid() && pid2 == pMix.pid()) || (pid1 == -p1.pid() && pid2 == -pMix.pid()))) { background[i]->fill(dPhi); nmp[i] += 1.0; } if (!samesign && abs(pid1) == abs(pid2) && pid1 == p1.pid() && pid2 == pMix.pid()) { background[i]->fill(dPhi); nmp[i] += 1.0; } - if (!samesign && abs(pid1) != abs(pid2) && + if (!samesign && abs(pid1) != abs(pid2) && ( (pid1 == p1.pid() && pid2 == pMix.pid()) || (pid2 == p1.pid() && pid1 == pMix.pid()) ) ) { background[i]->fill(dPhi); nmp[i] += 1.0; } } } } } } /// Normalise histograms etc., after the run void finalize() { for (int i = 0, N = pid.size(); i < N; ++i) { double sc = nmp[i] / nsp[i]; - signal[i]->scaleW(sc); - divide(signal[i],background[i],ratio[i]); + signal[i]->scaleW(sc); + divide(signal[i],background[i],ratio[i]); } } //@} /// @name Histograms //@{ vector > pid; vector signal; vector background; vector ratio; vector nsp; vector nmp; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ALICE_2016_I1507157); } 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/analyses/pluginMC/MC_Cent_pPb_Eta.cc.needspercentile b/analyses/pluginMC/MC_Cent_pPb_Eta.cc.needspercentile new file mode 100644 --- /dev/null +++ b/analyses/pluginMC/MC_Cent_pPb_Eta.cc.needspercentile @@ -0,0 +1,78 @@ +// -*- C++ -*- +#include "Rivet/Analyses/MC_Cent_pPb.hh" +#include "Rivet/Tools/Percentile.hh" + +namespace Rivet { + + +class MC_Cent_pPb_Eta : public Analysis { + +public: + + DEFAULT_RIVET_ANALYSIS_CTOR(MC_Cent_pPb_Eta); + + /// Book histograms and initialise projections before the run + void init() { + + MSG_INFO("CENT parameter set to " << getOption("cent","REF")); + + // The centrality projection. + declareCentrality(MC_SumETFwdPbCentrality(), + "MC_Cent_pPb_Calib", "SumETPb", "CENT"); + + // The trigger projection. + declare(MC_pPbMinBiasTrigger(), "Trigger"); + + // The particles to be analysed. + declare(ChargedFinalState(Cuts::eta > -2.7 && Cuts::eta < 2.7 && + Cuts::pT > 0.1*GeV), "CFS"); + + // The centrality bins and the corresponding histograms. + std::vector< std::pair > centralityBins = + { {0, 1}, {1, 5}, {5, 10}, {10, 20}, + {20, 30}, {30, 40}, {40, 60}, {60, 90} }; + // std::vector< std::tuple > refData = + // { {2, 1, 8}, {2, 1, 7}, {2, 1, 6}, {2, 1, 5}, + // {2, 1, 4}, {2, 1, 3}, {2, 1, 2}, {2, 1, 1} }; + std::vector< std::tuple > refData; + for ( int i = 8; i > 0; --i ) + refData.push_back(std::tuple(2, 1, i)); + + // The centrality-binned histograms. + _hEta = bookPercentile("CENT", centralityBins, refData); + + } + + /// Perform the per-event analysis + void analyze(const Event& event) { + const double weight = event.weight(); + + if ( !apply(event, "Trigger")() ) vetoEvent; + + _hEta->init(event); + for ( const auto &p : apply(event,"CFS").particles() ) + _hEta->fill(p.eta(), weight); + + } + + /// Finalize + void finalize() { + + // Scale by the inverse sum of event weights in each centrality + // bin. + _hEta->normalizePerEvent(); + + } + +private: + + /// The histograms binned in centrality. + Percentile _hEta; + +}; + + +// The hook for the plugin system +DECLARE_RIVET_PLUGIN(MC_Cent_pPb_Eta); + +} diff --git a/analyses/pluginMC/TEST.cc.needspercentile b/analyses/pluginMC/TEST.cc.needspercentile new file mode 100644 --- /dev/null +++ b/analyses/pluginMC/TEST.cc.needspercentile @@ -0,0 +1,83 @@ +// -*- C++ -*- +#include "Rivet/Analysis.hh" +#include "Rivet/Projections/ChargedFinalState.hh" +#include "Rivet/Projections/PrimaryParticles.hh" +#include "Rivet/Tools/Correlators.hh" + + +namespace Rivet { + + + class TEST : public CumulantAnalysis { + public: + + /// @name Constructors etc. + //@{ + + /// Constructor + TEST() : CumulantAnalysis("TEST") { + } + //@} + + public: + + /// @name Analysis methods + //@{ + /// Book histograms and initialise projections before the run + void init() { + + ChargedFinalState cfs(Cuts::abseta < 1.0); + declare(cfs, "CFS"); + ChargedFinalState pp(Cuts::abseta < 2.0); + declare(pp, "PP"); + book(h_c22, "c22",120,0,120); + book(h_c23, "c23",120,0,120); + book(h_v22pT, "v22pT",10,0,10); + ec22 = bookECorrelator<2,2>("ec22",h_c22); + ec23 = bookECorrelator<3,2>("ec32",h_c22); + ec22pT = bookECorrelator<2,2>("ec22pT",h_v22pT); + pair max = getMaxValues(); + // Declare correlator projections. + declare(Correlators(pp, max.first, max.second, h_v22pT),"CRS"); + } + /// Perform the per-event analysis + void analyze(const Event& event) { + const Correlators& c = apply(event,"CRS"); + ec22->fill(apply(event,"CFS").particles().size(), + c, event.weight()); + ec23->fill(apply(event,"CFS").particles().size(), + c, event.weight()); + ec22pT->fill(c, event.weight()); + } + /// Normalise histograms etc., after the run + void finalize() { + stream(); + cnTwoInt(h_c22,ec22); + cnTwoInt(h_c23,ec23); + vnTwoDiff(h_v22pT,ec22pT); + + } + + + //@} + private: + + + /// @name Histograms + //@{ + Scatter2DPtr h_c22; + Scatter2DPtr h_v22pT; + ECorrPtr ec22; + ECorrPtr ec22pT; + Scatter2DPtr h_c23; + ECorrPtr ec23; + //@} + + }; + + + + // The hook for the plugin system + DECLARE_RIVET_PLUGIN(TEST); + +} 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 diff --git a/include/Rivet/Projections/EventMixingFinalState.hh b/include/Rivet/Projections/EventMixingFinalState.hh --- a/include/Rivet/Projections/EventMixingFinalState.hh +++ b/include/Rivet/Projections/EventMixingFinalState.hh @@ -1,104 +1,225 @@ // -*- C++ -*- #ifndef RIVET_EventMixingFinalState_HH #define RIVET_EventMixingFinalState_HH #include "Rivet/Projection.hh" #include "Rivet/Projections/ParticleFinder.hh" +#include "Rivet/Tools/Random.hh" #include #include + namespace Rivet { - using std::deque; + // @brief Projects out an event mixed of several events, given + // a mixing observable (eg. number of final state particles), + // defining what should qualify as a mixable event. + // Binning in the mixing observable is defined in the constructor, + // as is the number of events one wants to mix with. + // The method calculateMixingObs() must can be overloaded + // in derived classes, to provide the definition of the mixing observable, + // on the provided projection, eg. centrality or something more elaborate. + // + // The implementation can correcly handle mixing of weighted events, but + // not multi-weighted events. Nor does the implementation attemt to handle + // calculation of one (or more) event weights for the mixed events. For most + // common use cases, like calculating a background sample, this is sufficient. + // If a more elaborate use case ever turns up, this must be reevaluated. + // + // + // @author Christian Bierlich + // Weighted random shuffle, similar to std::random_shuffle, which + // allows the passing of a weight for each element to be shuffled. + template + void weighted_shuffle(RandomAccessIterator first, RandomAccessIterator last, + WeightIterator fw, WeightIterator lw, RandomNumberGenerator& g) { + while(first != last && fw != lw) { + std::discrete_distribution weightDist(fw, lw); + int i = weightDist(g); + if(i){ + std::iter_swap(first, next(first, i)); + std::iter_swap(fw, next(fw, i)); + } + ++first; + ++fw; + } + } + // A MixEvent is a vector of particles with and associated weight. + typedef pair MixEvent; + typedef map > MixMap; + // EventMixingBase is the base class for event mixing projections. + // Most methods are defined in this base class as they should. + // In order to use it, a derived class should be implemented where: + // - The constructor is reimplmented, giving the derived projection type + // from which the mixing observable is calculated. The constructor must also + // be declared public in the derived class. + // - The calculateMixingObs is implemented. + // To examples of such derived classes are given below, + // 1) EventMixingFinalState, where the mixing observable are calculated + // on a multiplicity of a charged final state, and: + // 2) EventMixingCentrality, where the mixing observable is centrality. - /// Projects out an event mixed of several events, given - /// a mixing observable (eg. number of final state particles), - /// defining what should qualify as a mixable event. - /// Binning in the mixing observable is defined in the constructor, - /// as is the number of events one wants to mix with. - /// The protected method calculateMixingObs() can be overloaded - /// in derived classes, to define other mixing observables, eg. - /// centrality or something even more elaborate. - /// - /// @warning !!!DISCLAIMER!!! The projection makes no attempt at correct handling - /// of event weights - ie. what one should use as event weight for several - /// mixed events. Proceed with caution if you do not use an MC with - /// unit weights. - /// - /// @author Christian Bierlich - /// - typedef map > MixMap; - class EventMixingFinalState : public Projection { - public: + class EventMixingBase : public Projection { + protected: // Constructor - EventMixingFinalState(const ParticleFinder& fsp, const ParticleFinder& mix, size_t nMixIn, - double oMin, double oMax, double deltao ) : nMix(nMixIn){ - setName("EventMixingFinalState"); - declare(fsp, "FS"); - declare(mix, "MIX"); - MSG_WARNING("EventMixingFinalState is a naive implementation, not currently " << - "validated. Use with caution."); + EventMixingBase(const Projection* mixObsProjPtr, const ParticleFinder& mix, + size_t nMixIn, double oMin, double oMax, double deltao) : nMix(nMixIn), + unitWeights(true) { + // The base class contructor should be called explicitly in derived classes + // to add projections below. + setName("EventMixingBase"); + declare(*mixObsProjPtr,"OBS"); + declare(mix,"MIX"); + MSG_WARNING("EventMixing is not fully validated. Use with caution."); - // Set up the map for mixing events - for (double o = oMin; o < oMax; o+=deltao) - mixEvents[o] = deque(); + // Set up the map for mixing events. + for(double o = oMin; o < oMax; o+=deltao ) + mixEvents[o] = std::deque(); } - // Clone on the heap - DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); - - + public: + + // Test if we have enough mixing events available for projected, + // current mixing observable. + bool hasMixingEvents() const { + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + if(mixItr == mixEvents.end() || mixItr->second.size() < nMix + 1) + return false; + return true; + } + // Return a vector of mixing events. - vector getMixingEvents() const { - MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); - if(mixItr == mixEvents.end() || mixItr->second.size() < nMix + 1) - return vector(); - return vector(mixItr->second.begin(), mixItr->second.end() - 1); + vector getMixingEvents() const { + if (!hasMixingEvents()) + return vector(); + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + return vector(mixItr->second.begin(), mixItr->second.end() - 1); } + // Return a vector of particles from the mixing events. Can + // be overloaded in derived classes, though normally not neccesary. + virtual const Particles particles() const { + // Test if we have enough mixing events. + if(!hasMixingEvents()) + return vector(); + // Get mixing events for the current, projected mixing observable. + MixMap::const_iterator mixItr = mixEvents.lower_bound(mObs); + vector mixEvents = + vector(mixItr->second.begin(), mixItr->second.end() - 1); + // Make the vector of mixed particles. + Particles mixParticles; + vector weights; + size_t pSize = 0; + for(size_t i = 0; i < mixEvents.size(); ++i) + pSize+=mixEvents[i].first.size(); + mixParticles.reserve(pSize); + weights.reserve(pSize); + // Put the particles in the vector. + for(size_t i = 0; i < mixEvents.size(); ++i) { + mixParticles.insert(mixParticles.end(), mixEvents[i].first.begin(), + mixEvents[i].first.end()); + vector tmp = vector(mixEvents[i].first.size(), + mixEvents[i].second); + weights.insert(weights.end(), tmp.begin(), tmp.end()); + } + + // Shuffle the particles. + if (unitWeights) { + // Use the thread safe random number generator. + auto rnd = [&] (int i) {return rng()()%i;}; + random_shuffle(mixParticles.begin(), mixParticles.end(), rnd); + return mixParticles; + } + else { + weighted_shuffle(mixParticles.begin(), mixParticles.end(), + weights.begin(), weights.end(), rng()); + Particles tmp = vector(mixParticles.begin(), + mixParticles.begin() + size_t(ceil(mixParticles.size() / 2))); + return tmp; + } + } + protected: // Calulate mixing observable. - // Can be overloaded to define more elaborate observables. - virtual void calculateMixingObs(const Particles& parts){ - mObs = parts.size(); + // Must be overloaded in derived classes. + virtual void calculateMixingObs(const Projection* mProj) = 0; + + /// Perform the projection on the Event. + void project(const Event& e){ + const Projection* mixObsProjPtr = &applyProjection(e, "OBS"); + calculateMixingObs(mixObsProjPtr); + MixMap::iterator mixItr = mixEvents.lower_bound(mObs); + if(mixItr == mixEvents.end()){ + // We are out of bounds. + MSG_DEBUG("Mixing observable out of bounds."); + return; + } + const Particles mix = applyProjection(e, "MIX").particles(); + mixItr->second.push_back(make_pair(mix,e.weights()[0])); + // Assume unit weights until we see otherwise. + if (unitWeights && e.weights()[0] != 1.0 ) { + unitWeights = false; + nMix *= 2; + } + if(mixItr->second.size() > nMix + 1) + mixItr->second.pop_front(); } - /// Perform the projection on the Event - void project(const Event& e){ - const Particles parts = applyProjection(e, "FS").particles(); - - calculateMixingObs(parts); - - MixMap::iterator mixItr = mixEvents.lower_bound(mObs); - if(mixItr == mixEvents.end()){ - // We are out of bounds. - MSG_DEBUG("Mixing observable out of bounds."); - return; - } - const Particles mix = applyProjection(e, "MIX").particles(); - - mixItr->second.push_back(mix); - if(mixItr->second.size() > nMix + 1) - mixItr->second.pop_front(); + /// Compare with other projections. + CmpState compare(const Projection& p) const { + return mkNamedPCmp(p,"OBS"); } - - /// Compare with other projections - CmpState compare(const Projection& p) const { - return mkNamedPCmp(p, "FS"); - } - + + // The mixing observable of the current event. + double mObs; + private: - // The number of event to mix with + // The number of event to mix with. size_t nMix; - // The mixing observable of the current event - double mObs; - // The event map; + // The event map. MixMap mixEvents; + // Using unit weights or not. + bool unitWeights; }; + // EventMixingFinalState has multiplicity in the mixing projection. + // as the mixing observable. + class EventMixingFinalState : public EventMixingBase { + public: + EventMixingFinalState(const ParticleFinder* mixObsProjPtr, + const ParticleFinder& mix, size_t nMixIn, double oMin, double oMax, + double deltao) : + EventMixingBase(mixObsProjPtr, mix, nMixIn, oMin, oMax, deltao) { + setName("EventMixingFinalState"); + } + DEFAULT_RIVET_PROJ_CLONE(EventMixingFinalState); + + protected: + // Calulate mixing observable. + virtual void calculateMixingObs(const Projection* mProj) { + mObs = ((ParticleFinder*) mProj)->particles().size(); + } + }; + + // EventMixingCentrality has centrality as the mixing observable. + class EventMixingCentrality : public EventMixingBase { + public: + EventMixingCentrality(const CentralityProjection* mixObsProjPtr, + const ParticleFinder& mix, size_t nMixIn, double oMin, double oMax, + double deltao) : + EventMixingBase(mixObsProjPtr, mix, nMixIn, oMin, oMax, deltao) { + setName("EventMixingCentrality"); + } + + DEFAULT_RIVET_PROJ_CLONE(EventMixingCentrality); + protected: + virtual void calculateMixingObs(const Projection* mProj) { + mObs = ((CentralityProjection*) mProj)->operator()(); + } + }; } - #endif