diff --git a/analyses/pluginATLAS/ATLAS_2016_I1426523.cc b/analyses/pluginATLAS/ATLAS_2016_I1426523.cc --- a/analyses/pluginATLAS/ATLAS_2016_I1426523.cc +++ b/analyses/pluginATLAS/ATLAS_2016_I1426523.cc @@ -1,426 +1,426 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { /// @brief Measurement of the WZ production cross section at 8 TeV class ATLAS_2016_I1426523 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2016_I1426523); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Lepton cuts Cut FS_Zlept = Cuts::abseta < 2.5 && Cuts::pT > 15*GeV; const FinalState fs; Cut fs_z = Cuts::abseta < 2.5 && Cuts::pT > 15*GeV; Cut fs_j = Cuts::abseta < 4.5 && Cuts::pT > 25*GeV; // Get photons to dress leptons PromptFinalState photons(Cuts::abspid == PID::PHOTON); // Electrons and muons in Fiducial PS PromptFinalState leptons(fs_z && (Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON)); leptons.acceptTauDecays(false); DressedLeptons dressedleptons(photons, leptons, 0.1, FS_Zlept, true); declare(dressedleptons, "DressedLeptons"); // Electrons and muons in Total PS PromptFinalState leptons_total(Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON); leptons_total.acceptTauDecays(false); DressedLeptons dressedleptonsTotal(photons, leptons_total, 0.1, Cuts::open(), true); declare(dressedleptonsTotal, "DressedLeptonsTotal"); // Neutrinos IdentifiedFinalState nu_id; nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(false); declare(neutrinos, "Neutrinos"); MSG_WARNING("\033[91;1mLIMITED VALIDITY - check info file for details!\033[m"); // Jets VetoedFinalState veto; veto.addVetoOnThisFinalState(dressedleptons); FastJets jets(veto, FastJets::ANTIKT, 0.4); declare(jets, "Jets"); // Book histograms book(_h["eee"] , 1, 1, 1); book(_h["mee"] , 1, 1, 2); book(_h["emm"] , 1, 1, 3); book(_h["mmm"] , 1, 1, 4); book(_h["fid"] , 1, 1, 5); book(_h["eee_Plus"] , 2, 1, 1); book(_h["mee_Plus"] , 2, 1, 2); book(_h["emm_Plus"] , 2, 1, 3); book(_h["mmm_Plus"] , 2, 1, 4); book(_h["fid_Plus"] , 2, 1, 5); book(_h["eee_Minus"] , 3, 1, 1); book(_h["mee_Minus"] , 3, 1, 2); book(_h["emm_Minus"] , 3, 1, 3); book(_h["mmm_Minus"] , 3, 1, 4); book(_h["fid_Minus"] , 3, 1, 5); book(_h["total"] , 5, 1, 1); book(_h["Njets"] , 27, 1, 1); book(_h["Njets_norm"], 41, 1, 1); bookHandler("ZpT", 12); bookHandler("ZpT_Plus", 13); bookHandler("ZpT_Minus", 14); bookHandler("WpT", 15); bookHandler("WpT_Plus", 16); bookHandler("WpT_Minus", 17); bookHandler("mTWZ", 18); bookHandler("mTWZ_Plus", 19); bookHandler("mTWZ_Minus", 20); bookHandler("pTv", 21); bookHandler("pTv_Plus", 22); bookHandler("pTv_Minus", 23); bookHandler("Deltay", 24); bookHandler("Deltay_Plus", 25); bookHandler("Deltay_Minus", 26); bookHandler("mjj", 28); bookHandler("Deltayjj", 29); bookHandler("ZpT_norm", 30); bookHandler("ZpT_Plus_norm", 31); bookHandler("ZpT_Minus_norm", 32); bookHandler("WpT_norm", 33); bookHandler("mTWZ_norm", 34); bookHandler("pTv_norm", 35); bookHandler("pTv_Plus_norm", 36); bookHandler("pTv_Minus_norm", 37); bookHandler("Deltay_norm", 38); bookHandler("Deltay_Minus_norm", 39); bookHandler("Deltay_Plus_norm", 40); bookHandler("mjj_norm", 42); bookHandler("Deltayjj_norm", 43); } void bookHandler(const string& tag, size_t ID) { book(_s[tag], ID, 1, 1); - const string code1 = makeAxisCode(ID, 1, 1); - const string code2 = makeAxisCode(ID, 1, 2); + const string code1 = mkAxisCode(ID, 1, 1); + const string code2 = mkAxisCode(ID, 1, 2); book(_h[tag], code2, refData(code1)); } /// Perform the per-event analysis void analyze(const Event& event) { const vector& dressedleptons = apply(event, "DressedLeptons").dressedLeptons(); const vector& dressedleptonsTotal = apply(event, "DressedLeptonsTotal").dressedLeptons(); const Particles& neutrinos = apply(event, "Neutrinos").particlesByPt(); Jets jets = apply(event, "Jets").jetsByPt( (Cuts::abseta < 4.5) && (Cuts::pT > 25*GeV) ); if ((dressedleptonsTotal.size()<3) || (neutrinos.size()<1)) vetoEvent; //---Total PS: assign leptons to W and Z bosons using Resonant shape algorithm // NB: This resonant shape algorithm assumes the Standard Model and can therefore // NOT be used for reinterpretation in terms of new-physics models. int i, j, k; double MassZ01 = 0., MassZ02 = 0., MassZ12 = 0.; double MassW0 = 0., MassW1 = 0., MassW2 = 0.; double WeightZ1, WeightZ2, WeightZ3; double WeightW1, WeightW2, WeightW3; double M1, M2, M3; double WeightTotal1, WeightTotal2, WeightTotal3; //try Z pair of leptons 01 if ( (dressedleptonsTotal[0].pid()==-(dressedleptonsTotal[1].pid())) && (dressedleptonsTotal[2].abspid()==neutrinos[0].abspid()-1)){ MassZ01 = (dressedleptonsTotal[0].momentum()+dressedleptonsTotal[1].momentum()).mass(); MassW2 = (dressedleptonsTotal[2].momentum()+neutrinos[0].momentum()).mass(); } //try Z pair of leptons 02 if ( (dressedleptonsTotal[0].pid()==-(dressedleptonsTotal[2].pid())) && (dressedleptonsTotal[1].abspid()==neutrinos[0].abspid()-1)){ MassZ02 = (dressedleptonsTotal[0].momentum()+dressedleptonsTotal[2].momentum()).mass(); MassW1 = (dressedleptonsTotal[1].momentum()+neutrinos[0].momentum()).mass(); } //try Z pair of leptons 12 if ( (dressedleptonsTotal[1].pid()==-(dressedleptonsTotal[2].pid())) && (dressedleptonsTotal[0].abspid()==neutrinos[0].abspid()-1)){ MassZ12 = (dressedleptonsTotal[1].momentum()+dressedleptonsTotal[2].momentum()).mass(); MassW0 = (dressedleptonsTotal[0].momentum()+neutrinos[0].momentum()).mass(); } WeightZ1 = 1/(pow(MassZ01*MassZ01 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2)); WeightW1 = 1/(pow(MassW2*MassW2 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2)); WeightTotal1 = WeightZ1*WeightW1; M1 = -1*WeightTotal1; WeightZ2 = 1/(pow(MassZ02*MassZ02- MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2)); WeightW2 = 1/(pow(MassW1*MassW1- MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2)); WeightTotal2 = WeightZ2*WeightW2; M2 = -1*WeightTotal2; WeightZ3 = 1/(pow(MassZ12*MassZ12 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2)); WeightW3 = 1/(pow(MassW0*MassW0 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2)); WeightTotal3 = WeightZ3*WeightW3; M3 = -1*WeightTotal3; if( (M1 < M2 && M1 < M3) || (MassZ01 != 0 && MassW2 != 0 && MassZ02 == 0 && MassZ12 == 0) ){ i = 0; j = 1; k = 2; } if( (M2 < M1 && M2 < M3) || (MassZ02 != 0 && MassW1 != 0 && MassZ01 == 0 && MassZ12 == 0) ){ i = 0; j = 2; k = 1; } if( (M3 < M1 && M3 < M2) || (MassZ12 != 0 && MassW0 != 0 && MassZ01 == 0 && MassZ02 == 0) ){ i = 1; j = 2; k = 0; } FourMomentum ZbosonTotal = dressedleptonsTotal[i].momentum()+dressedleptonsTotal[j].momentum(); if ( (ZbosonTotal.mass() >= 66*GeV) && (ZbosonTotal.mass() <= 116*GeV) ) _h["total"]->fill(8000); //---end Total PS //---Fiducial PS: assign leptons to W and Z bosons using Resonant shape algorithm if (dressedleptons.size() < 3 || neutrinos.size() < 1) vetoEvent; int EventType = -1; int Nel = 0, Nmu = 0; for (const DressedLepton& l : dressedleptons) { if (l.abspid() == 11) ++Nel; if (l.abspid() == 13) ++Nmu; } if ( Nel == 3 && Nmu==0 ) EventType = 3; if ( Nel == 2 && Nmu==1 ) EventType = 2; if ( Nel == 1 && Nmu==2 ) EventType = 1; if ( Nel == 0 && Nmu==3 ) EventType = 0; int EventCharge = -dressedleptons[0].charge()*dressedleptons[1].charge()*dressedleptons[2].charge(); MassZ01 = 0; MassZ02 = 0; MassZ12 = 0; MassW0 = 0; MassW1 = 0; MassW2 = 0; //try Z pair of leptons 01 if ( (dressedleptons[0].pid()==-(dressedleptons[1].pid())) && (dressedleptons[2].abspid()==neutrinos[0].abspid()-1)){ MassZ01 = (dressedleptons[0].momentum()+dressedleptons[1].momentum()).mass(); MassW2 = (dressedleptons[2].momentum()+neutrinos[0].momentum()).mass(); } //try Z pair of leptons 02 if ( (dressedleptons[0].pid()==-(dressedleptons[2].pid())) && (dressedleptons[1].abspid()==neutrinos[0].abspid()-1)){ MassZ02 = (dressedleptons[0].momentum()+dressedleptons[2].momentum()).mass(); MassW1 = (dressedleptons[1].momentum()+neutrinos[0].momentum()).mass(); } //try Z pair of leptons 12 if ( (dressedleptons[1].pid()==-(dressedleptons[2].pid())) && (dressedleptons[0].abspid()==neutrinos[0].abspid()-1)){ MassZ12 = (dressedleptons[1].momentum()+dressedleptons[2].momentum()).mass(); MassW0 = (dressedleptons[0].momentum()+neutrinos[0].momentum()).mass(); } WeightZ1 = 1/(pow(MassZ01*MassZ01 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2)); WeightW1 = 1/(pow(MassW2*MassW2 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2)); WeightTotal1 = WeightZ1*WeightW1; M1 = -1*WeightTotal1; WeightZ2 = 1/(pow(MassZ02*MassZ02- MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2)); WeightW2 = 1/(pow(MassW1*MassW1- MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2)); WeightTotal2 = WeightZ2*WeightW2; M2 = -1*WeightTotal2; WeightZ3 = 1/(pow(MassZ12*MassZ12 - MZ_PDG*MZ_PDG,2) + pow(MZ_PDG*GammaZ_PDG,2)); WeightW3 = 1/(pow(MassW0*MassW0 - MW_PDG*MW_PDG,2) + pow(MW_PDG*GammaW_PDG,2)); WeightTotal3 = WeightZ3*WeightW3; M3 = -1*WeightTotal3; if( (M1 < M2 && M1 < M3) || (MassZ01 != 0 && MassW2 != 0 && MassZ02 == 0 && MassZ12 == 0) ){ i = 0; j = 1; k = 2; } if( (M2 < M1 && M2 < M3) || (MassZ02 != 0 && MassW1 != 0 && MassZ01 == 0 && MassZ12 == 0) ){ i = 0; j = 2; k = 1; } if( (M3 < M1 && M3 < M2) || (MassZ12 != 0 && MassW0 != 0 && MassZ01 == 0 && MassZ02 == 0) ){ i = 1; j = 2; k = 0; } FourMomentum Zlepton1 = dressedleptons[i].momentum(); FourMomentum Zlepton2 = dressedleptons[j].momentum(); FourMomentum Wlepton = dressedleptons[k].momentum(); FourMomentum Zboson = dressedleptons[i].momentum()+dressedleptons[j].momentum(); FourMomentum Wboson = dressedleptons[k].momentum()+neutrinos[0].momentum(); double Wboson_mT = sqrt( 2 * Wlepton.pT() * neutrinos[0].pt() * (1 - cos(deltaPhi(Wlepton, neutrinos[0]))) )/GeV; if (fabs(Zboson.mass()-MZ_PDG)>=10.) vetoEvent; if (Wboson_mT<=30.) vetoEvent; if (Wlepton.pT()<=20.) vetoEvent; if (deltaR(Zlepton1,Zlepton2) < 0.2) vetoEvent; if (deltaR(Zlepton1,Wlepton) < 0.3) vetoEvent; if (deltaR(Zlepton2,Wlepton) < 0.3) vetoEvent; double WZ_pt = Zlepton1.pt() + Zlepton2.pt() + Wlepton.pt() + neutrinos[0].pt(); double WZ_px = Zlepton1.px() + Zlepton2.px() + Wlepton.px() + neutrinos[0].px(); double WZ_py = Zlepton1.py() + Zlepton2.py() + Wlepton.py() + neutrinos[0].py(); double mTWZ = sqrt( pow(WZ_pt, 2) - ( pow(WZ_px, 2) + pow(WZ_py,2) ) )/GeV; double AbsDeltay = fabs(Zboson.rapidity()-Wlepton.rapidity()); if (EventType == 3) _h["eee"]->fill(8000.); if (EventType == 2) _h["mee"]->fill(8000.); if (EventType == 1) _h["emm"]->fill(8000.); if (EventType == 0) _h["mmm"]->fill(8000.); _h["fid"]->fill(8000.); if (EventCharge == 1) { if (EventType == 3) _h["eee_Plus"]->fill(8000.); if (EventType == 2) _h["mee_Plus"]->fill(8000.); if (EventType == 1) _h["emm_Plus"]->fill(8000.); if (EventType == 0) _h["mmm_Plus"]->fill(8000.); _h["fid_Plus"]->fill(8000.); _h["Deltay_Plus"]->fill(AbsDeltay); _h["Deltay_Plus_norm"]->fill(AbsDeltay); fillWithOverflow("ZpT_Plus", Zboson.pT()/GeV, 220); fillWithOverflow("WpT_Plus", Wboson.pT()/GeV, 220); fillWithOverflow("mTWZ_Plus", mTWZ, 600); fillWithOverflow("pTv_Plus", neutrinos[0].pt(), 90); fillWithOverflow("ZpT_Plus_norm", Zboson.pT()/GeV, 220); fillWithOverflow("pTv_Plus_norm", neutrinos[0].pt()/GeV, 90); } else { if (EventType == 3) _h["eee_Minus"]->fill(8000.); if (EventType == 2) _h["mee_Minus"]->fill(8000.); if (EventType == 1) _h["emm_Minus"]->fill(8000.); if (EventType == 0) _h["mmm_Minus"]->fill(8000.); _h["fid_Minus"]->fill(8000.); _h["Deltay_Minus"]->fill(AbsDeltay); _h["Deltay_Minus_norm"]->fill(AbsDeltay); fillWithOverflow("ZpT_Minus", Zboson.pT()/GeV, 220); fillWithOverflow("WpT_Minus", Wboson.pT()/GeV, 220); fillWithOverflow("mTWZ_Minus", mTWZ, 600); fillWithOverflow("pTv_Minus", neutrinos[0].pt()/GeV, 90); fillWithOverflow("ZpT_Minus_norm", Zboson.pT()/GeV, 220); fillWithOverflow("pTv_Minus_norm", neutrinos[0].pt()/GeV, 90); } fillWithOverflow("ZpT", Zboson.pT()/GeV, 220); fillWithOverflow("WpT", Wboson.pT()/GeV, 220); fillWithOverflow("mTWZ", mTWZ, 600); fillWithOverflow("pTv", neutrinos[0].pt()/GeV, 90); _h["Deltay"]->fill(AbsDeltay); fillWithOverflow("Njets", jets.size(), 5); fillWithOverflow("Njets_norm", jets.size(), 5); fillWithOverflow("ZpT_norm", Zboson.pT()/GeV, 220); fillWithOverflow("WpT_norm", Wboson.pT()/GeV, 220); fillWithOverflow("mTWZ_norm", mTWZ, 600); fillWithOverflow("pTv_norm", neutrinos[0].pt()/GeV, 90); _h["Deltay_norm"]->fill(AbsDeltay); if (jets.size()>1) { double mjj = (jets[0].momentum()+jets[1].momentum()).mass()/GeV; fillWithOverflow("mjj", mjj, 800); fillWithOverflow("mjj_norm", mjj, 800); double DeltaYjj = fabs(jets[0].rapidity()-jets[1].rapidity()); fillWithOverflow("Deltayjj", DeltaYjj, 5); fillWithOverflow("Deltayjj_norm", DeltaYjj, 5); } } void fillWithOverflow(const string& tag, const double value, const double overflow){ if (value < overflow) _h[tag]->fill(value); else _h[tag]->fill(overflow - 0.45); } /// Normalise histograms etc., after the run void finalize() { const double xs_pb(crossSection() / picobarn); const double xs_fb(crossSection() / femtobarn); const double sumw(sumOfWeights()); MSG_INFO("Cross-Section/pb: " << xs_pb ); MSG_INFO("Cross-Section/fb: " << xs_fb ); MSG_INFO("Sum of weights : " << sumw ); MSG_INFO("nEvents : " << numEvents()); const double sf_pb(xs_pb / sumw); const double sf_fb(xs_fb / sumw); MSG_INFO("sf_pb : " << sf_pb); MSG_INFO("sf_fb : " << sf_fb); float totalBR= 4*0.1086*0.033658; // W and Z leptonic branching fractions for (map::iterator it = _h.begin(); it != _h.end(); ++it) { if (it->first.find("total") != string::npos) scale(it->second, sf_pb/totalBR); else if (it->first.find("norm") != string::npos) normalize(it->second); else if (it->first.find("fid") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("Njets") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("ZpT") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("WpT") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("mTWZ") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("pTv") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("Deltay") != string::npos) scale(it->second, sf_fb/4.); else if (it->first.find("mjj") != string::npos) scale(it->second, sf_fb/4.); else scale(it->second, sf_fb); } for (map::iterator it = _s.begin(); it != _s.end(); ++it) { makeScatterWithoutDividingByBinwidth(it->first); removeAnalysisObject(_h[it->first]); } } void makeScatterWithoutDividingByBinwidth(const string& tag) { vector points; //size_t nBins = _dummy->numBins(); for (const HistoBin1D &bin : _h[tag]->bins()) { double x = bin.midpoint(); double y = bin.sumW(); double ex = bin.xWidth()/2; double ey = sqrt(bin.sumW2()); points.push_back(Point2D(x, y, ex, ey)); } _s[tag]->addPoints(points); } //@} private: /// @name Histograms //@{ map _h; map _s; //@} double MZ_PDG = 91.1876; double MW_PDG = 83.385; double GammaZ_PDG = 2.4952; double GammaW_PDG = 2.085; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2016_I1426523); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1609448.cc b/analyses/pluginATLAS/ATLAS_2017_I1609448.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1609448.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1609448.cc @@ -1,285 +1,285 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" namespace Rivet { /// ATLAS pTmiss+jets cross-section ratios at 13 TeV class ATLAS_2017_I1609448 : public Analysis { public: /// Constructor ATLAS_2017_I1609448(string name="ATLAS_2017_I1609448") : Analysis(name) { _mode = 0; // using Z -> nunu channel by default } struct HistoHandler { Histo1DPtr histo; Scatter2DPtr scatter; unsigned int d, x, y; void fill(double value) { histo->fill(value); } }; /// Initialize void init() { // Prompt photons PromptFinalState photon_fs(Cuts::abspid == PID::PHOTON && Cuts::abseta < 4.9); // Prompt electrons PromptFinalState el_fs(Cuts::abseta < 4.9 && Cuts::abspid == PID::ELECTRON); // Prompt muons PromptFinalState mu_fs(Cuts::abseta < 4.9 && Cuts::abspid == PID::MUON); // Dressed leptons Cut lep_cuts = Cuts::pT > 7*GeV && Cuts::abseta < 2.5; DressedLeptons dressed_leps(photon_fs, (_mode == 2 ? el_fs : mu_fs), 0.1, lep_cuts); declare(dressed_leps, "DressedLeptons"); // In-acceptance leptons for lepton veto PromptFinalState veto_lep_fs(Cuts::abseta < 4.9 && (Cuts::abspid == PID::ELECTRON || Cuts::abspid == PID::MUON)); veto_lep_fs.acceptTauDecays(); veto_lep_fs.acceptMuonDecays(); DressedLeptons veto_lep(photon_fs, veto_lep_fs, 0.1, lep_cuts); declare(veto_lep, "VetoLeptons"); // MET VetoedFinalState met_fs(!(Cuts::abseta > 2.5 && Cuts::abspid == PID::MUON)); if (_mode) met_fs.addVetoOnThisFinalState(dressed_leps); declare(MissingMomentum(met_fs), "MET"); // Jet collection FastJets jets(FinalState(Cuts::abseta < 4.9), FastJets::ANTIKT, 0.4, JetAlg::Muons::NONE, JetAlg::Invisibles::NONE); declare(jets, "Jets"); _h["met_mono"] = bookHandler(1, 1, 2); _h["met_vbf" ] = bookHandler(2, 1, 2); _h["mjj_vbf" ] = bookHandler(3, 1, 2); _h["dphijj_vbf"] = bookHandler(4, 1, 2); } HistoHandler bookHandler(unsigned int id_d, unsigned int id_x, unsigned int id_y) { HistoHandler dummy; if (_mode < 2) { // numerator mode - const string histName = "_" + makeAxisCode(id_d, id_x, id_y); + const string histName = "_" + mkAxisCode(id_d, id_x, id_y); book(dummy.histo, histName, refData(id_d, id_x, id_y)); // hidden auxiliary output book(dummy.scatter, id_d, id_x, id_y - 1, true); // ratio dummy.d = id_d; dummy.x = id_x; dummy.y = id_y; } else { book(dummy.histo, id_d, id_x, 4); // denominator mode } return dummy; } bool isBetweenJets(const Jet& probe, const Jet& boundary1, const Jet& boundary2) { const double y_p = probe.rapidity(); const double y_b1 = boundary1.rapidity(); const double y_b2 = boundary2.rapidity(); const double y_min = std::min(y_b1, y_b2); const double y_max = std::max(y_b1, y_b2); return (y_p > y_min && y_p < y_max); } int centralJetVeto(Jets& jets) { if (jets.size() < 2) return 0; const Jet bj1 = jets.at(0); const Jet bj2 = jets.at(1); // Start loop at the 3rd hardest pT jet int n_between = 0; for (size_t i = 2; i < jets.size(); ++i) { const Jet j = jets.at(i); if (isBetweenJets(j, bj1, bj2) && j.pT() > 25*GeV) ++n_between; } return n_between; } /// Perform the per-event analysis void analyze(const Event& event) { // Require 0 (Znunu) or 2 (Zll) dressed leptons bool isZll = bool(_mode); const vector &vetoLeptons = applyProjection(event, "VetoLeptons").dressedLeptons(); const vector &all_leps = applyProjection(event, "DressedLeptons").dressedLeptons(); if (!isZll && vetoLeptons.size()) vetoEvent; if ( isZll && all_leps.size() != 2) vetoEvent; vector leptons; bool pass_Zll = true; if (isZll) { // Sort dressed leptons by pT if (all_leps[0].pt() > all_leps[1].pt()) { leptons.push_back(all_leps[0]); leptons.push_back(all_leps[1]); } else { leptons.push_back(all_leps[1]); leptons.push_back(all_leps[0]); } // Leading lepton pT cut pass_Zll &= leptons[0].pT() > 80*GeV; // Opposite-charge requirement pass_Zll &= charge3(leptons[0]) + charge3(leptons[1]) == 0; // Z-mass requirement const double Zmass = (leptons[0].mom() + leptons[1].mom()).mass(); pass_Zll &= (Zmass >= 66*GeV && Zmass <= 116*GeV); } if (!pass_Zll) vetoEvent; // Get jets and remove those within dR = 0.5 of a dressed lepton Jets jets = applyProjection(event, "Jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::absrap < 4.4); for (const DressedLepton& lep : leptons) ifilter_discard(jets, deltaRLess(lep, 0.5)); const size_t njets = jets.size(); if (!njets) vetoEvent; const int njets_gap = centralJetVeto(jets); double jpt1 = jets[0].pT(); double jeta1 = jets[0].eta(); double mjj = 0., jpt2 = 0., dphijj = 0.; if (njets >= 2) { mjj = (jets[0].momentum() + jets[1].momentum()).mass(); jpt2 = jets[1].pT(); dphijj = deltaPhi(jets[0], jets[1]); } // MET Vector3 met_vec = apply(event, "MET").vectorMPT(); double met = met_vec.mod(); // Cut on deltaPhi between MET and first 4 jets, but only if jet pT > 30 GeV bool dphi_fail = false; for (size_t i = 0; i < jets.size() && i < 4; ++i) { dphi_fail |= (deltaPhi(jets[i], met_vec) < 0.4 && jets[i].pT() > 30*GeV); } const bool pass_met_dphi = met > 200*GeV && !dphi_fail; const bool pass_vbf = pass_met_dphi && mjj > 200*GeV && jpt1 > 80*GeV && jpt2 > 50*GeV && njets >= 2 && !njets_gap; const bool pass_mono = pass_met_dphi && jpt1 > 120*GeV && fabs(jeta1) < 2.4; if (pass_mono) _h["met_mono"].fill(met); if (pass_vbf) { _h["met_vbf"].fill(met/GeV); _h["mjj_vbf"].fill(mjj/GeV); _h["dphijj_vbf"].fill(dphijj); } } /// Normalise, scale and otherwise manipulate histograms here void finalize() { const double sf(crossSection() / femtobarn / sumOfWeights()); for (map::iterator hit = _h.begin(); hit != _h.end(); ++hit) { scale(hit->second.histo, sf); if (_mode < 2) constructRmiss(hit->second); } } void constructRmiss(HistoHandler& handler) { // Load transfer function from reference data file const YODA::Scatter2D& rmiss = refData(handler.d, handler.x, handler.y); const YODA::Scatter2D& numer = refData(handler.d, handler.x, handler.y + 1); const YODA::Scatter2D& denom = refData(handler.d, handler.x, handler.y + 2); for (size_t i = 0; i < handler.scatter->numPoints(); ++i) { const Point2D& r = rmiss.point(i); // SM Rmiss const Point2D& n = numer.point(i); // SM numerator const Point2D& d = denom.point(i); // SM denominator const HistoBin1D& b = handler.histo->bin(i); // BSM double bsmy; try { bsmy = b.height(); } catch (const Exception&) { // LowStatsError or WeightError bsmy = 0; } double bsmey; try { bsmey = b.heightErr(); } catch (const Exception&) { // LowStatsError or WeightError bsmey = 0; } // Combined numerator double sm_plus_bsm = n.y() + bsmy; // Rmiss central value double rmiss_y = safediv(sm_plus_bsm, d.y()); // Ratio error (Rmiss = SM_num/SM_denom + BSM/SM_denom ~ Rmiss_SM + BSM/SM_denom double rmiss_p = sqrt(r.yErrPlus()*r.yErrPlus() + safediv(bsmey*bsmey, d.y()*d.y())); double rmiss_m = sqrt(r.yErrMinus()*r.yErrMinus() + safediv(bsmey*bsmey, d.y()*d.y())); // Set new values Point2D& p = handler.scatter->point(i); // (SM + BSM) Rmiss p.setY(rmiss_y); p.setYErrMinus(rmiss_m); p.setYErrPlus(rmiss_p); } } protected: // Analysis-mode switch size_t _mode; /// Histograms map _h; }; /// ATLAS pTmiss+jets specialisation for Znunu channel class ATLAS_2017_I1609448_Znunu : public ATLAS_2017_I1609448 { public: ATLAS_2017_I1609448_Znunu() : ATLAS_2017_I1609448("ATLAS_2017_I1609448_Znunu") { _mode = 0; } }; /// ATLAS pTmiss+jets specialisation for Zmumu channel class ATLAS_2017_I1609448_Zmumu : public ATLAS_2017_I1609448 { public: ATLAS_2017_I1609448_Zmumu() : ATLAS_2017_I1609448("ATLAS_2017_I1609448_Zmumu") { _mode = 1; } }; /// ATLAS pTmiss+jets specialisation for Zee channel class ATLAS_2017_I1609448_Zee : public ATLAS_2017_I1609448 { public: ATLAS_2017_I1609448_Zee() : ATLAS_2017_I1609448("ATLAS_2017_I1609448_Zee") { _mode = 2; } }; // Hooks for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1609448); DECLARE_RIVET_PLUGIN(ATLAS_2017_I1609448_Znunu); DECLARE_RIVET_PLUGIN(ATLAS_2017_I1609448_Zmumu); DECLARE_RIVET_PLUGIN(ATLAS_2017_I1609448_Zee); } diff --git a/include/Rivet/Analysis.hh b/include/Rivet/Analysis.hh --- a/include/Rivet/Analysis.hh +++ b/include/Rivet/Analysis.hh @@ -1,961 +1,956 @@ // -*- 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" /// @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::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. virtual std::string name() const { return (info().name().empty()) ? _defaultname : info().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; //@} /// @name Analysis / beam compatibility testing //@{ /// 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 sumOfWeights() 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; - /// Alias - /// @deprecated Prefer the "mk" form, consistent with other "making function" names - const std::string makeAxisCode(unsigned int datasetId, unsigned int xAxisId, unsigned int yAxisId) const { - return mkAxisCode(datasetId, xAxisId, yAxisId); - } //@} /// @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 = makeAxisCode(datasetId, xAxisId, yAxisId); + 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=""); /// Book a 2D profile histogram with binning from a reference scatter. // Profile2DPtr bookProfile2D(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 bookProfile2D(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 bookProfile2D(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); //@} 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, and removal //@{ /// Register a data object in the histogram system void addAnalysisObject(const MultiweightAOPtr & ao); /// 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); //@} 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; //@} /// 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; 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 /// @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