diff --git a/analyses/pluginATLAS/ATLAS_2012_I1083318.cc b/analyses/pluginATLAS/ATLAS_2012_I1083318.cc --- a/analyses/pluginATLAS/ATLAS_2012_I1083318.cc +++ b/analyses/pluginATLAS/ATLAS_2012_I1083318.cc @@ -1,252 +1,245 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" namespace Rivet { /// ATLAS W + jets production at 7 TeV class ATLAS_2012_I1083318 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2012_I1083318); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs; IdentifiedFinalState allleptons; allleptons.acceptIdPair(PID::ELECTRON); allleptons.acceptIdPair(PID::MUON); Cut cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV; DressedLeptons leptons(fs, allleptons, 0.1, cuts); declare(leptons, "leptons"); // Leading neutrinos for Etmiss LeadingParticlesFinalState neutrinos(fs); neutrinos.addParticleIdPair(PID::NU_E); neutrinos.addParticleIdPair(PID::NU_MU); neutrinos.setLeadingOnly(true); declare(neutrinos, "neutrinos"); // Input for the jets: "Neutrinos, electrons, and muons from decays of the // massive W boson were not used" VetoedFinalState veto; veto.addVetoOnThisFinalState(leptons); veto.addVetoOnThisFinalState(neutrinos); - FastJets jets(veto, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(veto, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); for (size_t i = 0; i < 2; ++i) { book(_h_NjetIncl[i] ,1, 1, i+1); book(_h_RatioNjetIncl[i], 2, 1, i+1); book(_h_FirstJetPt_1jet[i] ,3, 1, i+1); book(_h_FirstJetPt_2jet[i] ,4, 1, i+1); book(_h_FirstJetPt_3jet[i] ,5, 1, i+1); book(_h_FirstJetPt_4jet[i] ,6, 1, i+1); book(_h_SecondJetPt_2jet[i] ,7, 1, i+1); book(_h_SecondJetPt_3jet[i] ,8, 1, i+1); book(_h_SecondJetPt_4jet[i] ,9, 1, i+1); book(_h_ThirdJetPt_3jet[i] ,10, 1, i+1); book(_h_ThirdJetPt_4jet[i] ,11, 1, i+1); book(_h_FourthJetPt_4jet[i] ,12, 1, i+1); book(_h_Ht_1jet[i] ,13, 1, i+1); book(_h_Ht_2jet[i] ,14, 1, i+1); book(_h_Ht_3jet[i] ,15, 1, i+1); book(_h_Ht_4jet[i] ,16, 1, i+1); book(_h_Minv_2jet[i] ,17, 1, i+1); book(_h_Minv_3jet[i] ,18, 1, i+1); book(_h_Minv_4jet[i] ,19, 1, i+1); book(_h_JetRapidity[i] ,20, 1, i+1); book(_h_DeltaYElecJet[i] ,21, 1, i+1); book(_h_SumYElecJet[i] ,22, 1, i+1); book(_h_DeltaR_2jet[i] ,23, 1, i+1); book(_h_DeltaY_2jet[i] ,24, 1, i+1); book(_h_DeltaPhi_2jet[i] ,25, 1, i+1); } } /// Perform the per-event analysis void analyze(const Event& event) { const vector& leptons = apply(event, "leptons").dressedLeptons(); Particles neutrinos = apply(event, "neutrinos").particlesByPt(); - if (leptons.size() != 1 || (neutrinos.size() == 0)) { - vetoEvent; - } + if (leptons.size() != 1 || (neutrinos.size() == 0)) vetoEvent; FourMomentum lepton = leptons[0].momentum(); FourMomentum p_miss = neutrinos[0].momentum(); - if (p_miss.Et() < 25.0*GeV) { - vetoEvent; - } + if (p_miss.Et() < 25.0*GeV) vetoEvent; double mT = sqrt(2.0 * lepton.pT() * p_miss.Et() * (1.0 - cos( lepton.phi()-p_miss.phi()) ) ); - if (mT < 40.0*GeV) { - vetoEvent; - } + if (mT < 40.0*GeV) vetoEvent; double jetcuts[] = { 30.0*GeV, 20.0*GeV }; const FastJets& jetpro = apply(event, "jets"); for (size_t i = 0; i < 2; ++i) { vector jets; double HT = lepton.pT() + p_miss.pT(); for (const Jet& jet : jetpro.jetsByPt(jetcuts[i])) { if (jet.absrap() < 4.4 && deltaR(lepton, jet.momentum()) > 0.5) { jets.push_back(jet.momentum()); HT += jet.pT(); } } _h_NjetIncl[i]->fill(0.0); // Njet>=1 observables if (jets.size() < 1) continue; _h_NjetIncl[i]->fill(1.0); _h_FirstJetPt_1jet[i]->fill(jets[0].pT()); _h_JetRapidity[i]->fill(jets[0].rapidity()); _h_Ht_1jet[i]->fill(HT); _h_DeltaYElecJet[i]->fill(lepton.rapidity()-jets[0].rapidity()); _h_SumYElecJet[i]->fill(lepton.rapidity()+jets[0].rapidity()); // Njet>=2 observables if (jets.size() < 2) continue; _h_NjetIncl[i]->fill(2.0); _h_FirstJetPt_2jet[i]->fill(jets[0].pT()); _h_SecondJetPt_2jet[i]->fill(jets[1].pT()); _h_Ht_2jet[i]->fill(HT); double m2_2jet = FourMomentum(jets[0]+jets[1]).mass2(); _h_Minv_2jet[i]->fill(m2_2jet>0.0 ? sqrt(m2_2jet) : 0.0); _h_DeltaR_2jet[i]->fill(deltaR(jets[0], jets[1])); _h_DeltaY_2jet[i]->fill(jets[0].rapidity()-jets[1].rapidity()); _h_DeltaPhi_2jet[i]->fill(deltaPhi(jets[0], jets[1])); // Njet>=3 observables if (jets.size() < 3) continue; _h_NjetIncl[i]->fill(3.0); _h_FirstJetPt_3jet[i]->fill(jets[0].pT()); _h_SecondJetPt_3jet[i]->fill(jets[1].pT()); _h_ThirdJetPt_3jet[i]->fill(jets[2].pT()); _h_Ht_3jet[i]->fill(HT); double m2_3jet = FourMomentum(jets[0]+jets[1]+jets[2]).mass2(); _h_Minv_3jet[i]->fill(m2_3jet>0.0 ? sqrt(m2_3jet) : 0.0); // Njet>=4 observables if (jets.size() < 4) continue; _h_NjetIncl[i]->fill(4.0); _h_FirstJetPt_4jet[i]->fill(jets[0].pT()); _h_SecondJetPt_4jet[i]->fill(jets[1].pT()); _h_ThirdJetPt_4jet[i]->fill(jets[2].pT()); _h_FourthJetPt_4jet[i]->fill(jets[3].pT()); _h_Ht_4jet[i]->fill(HT); double m2_4jet = FourMomentum(jets[0]+jets[1]+jets[2]+jets[3]).mass2(); _h_Minv_4jet[i]->fill(m2_4jet>0.0 ? sqrt(m2_4jet) : 0.0); // Njet>=5 observables if (jets.size() < 5) continue; _h_NjetIncl[i]->fill(5.0); } } /// Normalise histograms etc., after the run void finalize() { for (size_t i = 0; i < 2; ++i) { // Construct jet multiplicity ratio for (size_t n = 1; n < _h_NjetIncl[i]->numBins(); ++n) { YODA::HistoBin1D& b0 = _h_NjetIncl[i]->bin(n-1); YODA::HistoBin1D& b1 = _h_NjetIncl[i]->bin(n); double val = 0.0, err= 0.0; if (b0.height() && b1.height()) { val = b1.height() / b0.height(); err = b1.height() / b0.height() * (b0.relErr() + b1.relErr()); } _h_RatioNjetIncl[i]->addPoint(n, val, 0.5, err); } // Scale all histos to the cross section const double factor = crossSection()/sumOfWeights(); scale(_h_DeltaPhi_2jet[i], factor); scale(_h_DeltaR_2jet[i], factor); scale(_h_DeltaY_2jet[i], factor); scale(_h_DeltaYElecJet[i], factor); scale(_h_FirstJetPt_1jet[i], factor); scale(_h_FirstJetPt_2jet[i], factor); scale(_h_FirstJetPt_3jet[i], factor); scale(_h_FirstJetPt_4jet[i], factor); scale(_h_FourthJetPt_4jet[i], factor); scale(_h_Ht_1jet[i], factor); scale(_h_Ht_2jet[i], factor); scale(_h_Ht_3jet[i], factor); scale(_h_Ht_4jet[i], factor); scale(_h_JetRapidity[i], factor); scale(_h_Minv_2jet[i], factor); scale(_h_Minv_3jet[i], factor); scale(_h_Minv_4jet[i], factor); scale(_h_NjetIncl[i], factor); scale(_h_SecondJetPt_2jet[i], factor); scale(_h_SecondJetPt_3jet[i], factor); scale(_h_SecondJetPt_4jet[i], factor); scale(_h_SumYElecJet[i], factor); scale(_h_ThirdJetPt_3jet[i], factor); scale(_h_ThirdJetPt_4jet[i], factor); } } //@} private: /// @name Histograms //@{ Histo1DPtr _h_DeltaPhi_2jet[2]; Histo1DPtr _h_DeltaR_2jet[2]; Histo1DPtr _h_DeltaY_2jet[2]; Histo1DPtr _h_DeltaYElecJet[2]; Histo1DPtr _h_FirstJetPt_1jet[2]; Histo1DPtr _h_FirstJetPt_2jet[2]; Histo1DPtr _h_FirstJetPt_3jet[2]; Histo1DPtr _h_FirstJetPt_4jet[2]; Histo1DPtr _h_FourthJetPt_4jet[2]; Histo1DPtr _h_Ht_1jet[2]; Histo1DPtr _h_Ht_2jet[2]; Histo1DPtr _h_Ht_3jet[2]; Histo1DPtr _h_Ht_4jet[2]; Histo1DPtr _h_JetRapidity[2]; Histo1DPtr _h_Minv_2jet[2]; Histo1DPtr _h_Minv_3jet[2]; Histo1DPtr _h_Minv_4jet[2]; Histo1DPtr _h_NjetIncl[2]; Scatter2DPtr _h_RatioNjetIncl[2]; Histo1DPtr _h_SecondJetPt_2jet[2]; Histo1DPtr _h_SecondJetPt_3jet[2]; Histo1DPtr _h_SecondJetPt_4jet[2]; Histo1DPtr _h_SumYElecJet[2]; Histo1DPtr _h_ThirdJetPt_3jet[2]; Histo1DPtr _h_ThirdJetPt_4jet[2]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2012_I1083318); } diff --git a/analyses/pluginATLAS/ATLAS_2013_I1217863_W.cc b/analyses/pluginATLAS/ATLAS_2013_I1217863_W.cc --- a/analyses/pluginATLAS/ATLAS_2013_I1217863_W.cc +++ b/analyses/pluginATLAS/ATLAS_2013_I1217863_W.cc @@ -1,198 +1,197 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" namespace Rivet { class ATLAS_2013_I1217863_W : public Analysis { public: /// Constructor ATLAS_2013_I1217863_W(string name="ATLAS_2013_I1217863_W") : Analysis(name) { // the electron mode is used by default _mode = 1; } /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs; declare(fs, "FS"); Cut cuts = Cuts::abseta < 2.47 && Cuts::pT > 25*GeV; // W finder for electrons and muons WFinder wf(fs, cuts, _mode==3? PID::MUON : PID::ELECTRON, 0.0*GeV, 1000.0*GeV, 35.0*GeV, 0.1, WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT); declare(wf, "WF"); // leading photon LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < 2.37 && Cuts::pT > 15*GeV)); photonfs.addParticleId(PID::PHOTON); declare(photonfs, "LeadingPhoton"); // jets VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("WF")); jet_fs.addVetoOnThisFinalState(getProjection("LeadingPhoton")); - FastJets jets(jet_fs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(jet_fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "Jets"); // FS excluding the leading photon VetoedFinalState vfs(fs); vfs.addVetoOnThisFinalState(photonfs); declare(vfs, "isolatedFS"); // Book histograms book(_hist_EgammaT_incl , 7, 1, _mode); // dSigma / dE^gamma_T for Njet >= 0 book(_hist_EgammaT_excl , 8, 1, _mode); // dSigma / dE^gamma_T for Njet = 0 book(_hist_Njet_EgammaT15 ,15, 1, _mode); // dSigma / dNjet for E^gamma_T > 15 book(_hist_Njet_EgammaT60 ,16, 1, _mode); // dSigma / dNjet for E^gamma_T > 60 book(_hist_mWgammaT ,19, 1, _mode); // dSigma / dm^{Wgamma}_T } /// Perform the per-event analysis void analyze(const Event& event) { // retrieve leading photon Particles photons = apply(event, "LeadingPhoton").particles(); if (photons.size() != 1) vetoEvent; const Particle& leadingPhoton = photons[0]; if (leadingPhoton.Et() < 15.0*GeV) vetoEvent; if (leadingPhoton.abseta() > 2.37) vetoEvent; // check photon isolation double coneEnergy(0.0); Particles fs = apply(event, "isolatedFS").particles(); for (const Particle& p : fs) { if ( deltaR(leadingPhoton, p) < 0.4 ) coneEnergy += p.E(); } if ( coneEnergy / leadingPhoton.E() >= 0.5 ) vetoEvent; // retrieve W boson candidate const WFinder& wf = apply(event, "WF"); if ( wf.bosons().size() != 1 ) vetoEvent; // only one W boson candidate //const Particle& Wboson = wf.boson(); // retrieve constituent neutrino const Particle& neutrino = wf.constituentNeutrino(); if ( !(neutrino.pT() > 35.0*GeV) ) vetoEvent; // retrieve constituent lepton const Particle& lepton = wf.constituentLepton(); if ( !(lepton.pT() > 25.0*GeV && lepton.abseta() < 2.47) ) vetoEvent; // check photon-lepton overlap if ( !(deltaR(leadingPhoton, lepton) > 0.7) ) vetoEvent; // count jets const FastJets& jetfs = apply(event, "Jets"); Jets jets = jetfs.jets(cmpMomByEt); int goodJets = 0; for (const Jet& j : jets) { if ( !(j.Et() > 30.0*GeV) ) break; if ( (j.abseta() < 4.4) && \ (deltaR(leadingPhoton, j) > 0.3) && \ (deltaR(lepton, j) > 0.3) ) ++goodJets; } double Njets = double(goodJets) + 0.5; double photonEt = leadingPhoton.Et()*GeV; const FourMomentum& lep_gamma = lepton.momentum() + leadingPhoton.momentum(); double term1 = sqrt(lep_gamma.mass2() + lep_gamma.pT2()) + neutrino.Et(); double term2 = (lep_gamma + neutrino.momentum()).pT2(); double mWgammaT = sqrt(term1 * term1 - term2) * GeV; _hist_EgammaT_incl->fill(photonEt); _hist_Njet_EgammaT15->fill(Njets); if ( !goodJets ) _hist_EgammaT_excl->fill(photonEt); if (photonEt > 40.0*GeV) { _hist_mWgammaT->fill(mWgammaT); if (photonEt > 60.0*GeV) _hist_Njet_EgammaT60->fill(Njets); } } /// Normalise histograms etc., after the run void finalize() { const double xs_fb = crossSection()/femtobarn; const double sumw = sumOfWeights(); const double sf = xs_fb / sumw; scale(_hist_EgammaT_excl, sf); scale(_hist_EgammaT_incl, sf); normalize(_hist_Njet_EgammaT15); normalize(_hist_Njet_EgammaT60); normalize(_hist_mWgammaT); } //@} protected: size_t _mode; private: /// @name Histograms //@{ Histo1DPtr _hist_EgammaT_incl; Histo1DPtr _hist_EgammaT_excl; Histo1DPtr _hist_Njet_EgammaT15; Histo1DPtr _hist_Njet_EgammaT60; Histo1DPtr _hist_mWgammaT; //@} }; class ATLAS_2013_I1217863_W_EL : public ATLAS_2013_I1217863_W { public: ATLAS_2013_I1217863_W_EL() : ATLAS_2013_I1217863_W("ATLAS_2013_I1217863_W_EL") { _mode = 2; } }; class ATLAS_2013_I1217863_W_MU : public ATLAS_2013_I1217863_W { public: ATLAS_2013_I1217863_W_MU() : ATLAS_2013_I1217863_W("ATLAS_2013_I1217863_W_MU") { _mode = 3; } }; DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217863_W); DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217863_W_EL); DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217863_W_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2013_I1217863_Z.cc b/analyses/pluginATLAS/ATLAS_2013_I1217863_Z.cc --- a/analyses/pluginATLAS/ATLAS_2013_I1217863_Z.cc +++ b/analyses/pluginATLAS/ATLAS_2013_I1217863_Z.cc @@ -1,193 +1,192 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" namespace Rivet { class ATLAS_2013_I1217863_Z : public Analysis { public: /// Constructor ATLAS_2013_I1217863_Z(string name="ATLAS_2013_I1217863_Z") : Analysis(name) { // the electron mode is used by default _mode = 1; } /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs; declare(fs, "FS"); Cut cuts = Cuts::abseta < 2.47 && Cuts::pT > 25*GeV; // Z finder ZFinder zf(fs, cuts, _mode==3? PID::MUON : PID::ELECTRON, 40.0*GeV, 1000.0*GeV, 0.1, ZFinder::ClusterPhotons::NODECAY, ZFinder::AddPhotons::NO); declare(zf, "ZF"); // leading photon LeadingParticlesFinalState photonfs(FinalState(Cuts::abseta < 2.37 && Cuts::pT > 15*GeV)); photonfs.addParticleId(PID::PHOTON); declare(photonfs, "LeadingPhoton"); // jets VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("ZF")); jet_fs.addVetoOnThisFinalState(getProjection("LeadingPhoton")); - FastJets jets(jet_fs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(jet_fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "Jets"); // FS excluding the leading photon VetoedFinalState vfs(fs); vfs.addVetoOnThisFinalState(photonfs); declare(vfs, "isolatedFS"); // Book histograms book(_hist_EgammaT_incl ,11, 1, _mode); // dSigma / dE^gamma_T for Njet >= 0 book(_hist_EgammaT_excl ,12, 1, _mode); // dSigma / dE^gamma_T for Njet = 0 book(_hist_Njet_EgammaT15 ,17, 1, _mode); // dSigma / dNjet for E^gamma_T >= 15 book(_hist_Njet_EgammaT60 ,18, 1, _mode); // dSigma / dNjet for E^gamma_T >= 60 book(_hist_mZgamma ,20, 1, _mode); // dSigma / dm^{Zgamma} } /// Perform the per-event analysis void analyze(const Event& event) { // retrieve leading photon Particles photons = apply(event, "LeadingPhoton").particles(); if (photons.size() != 1) vetoEvent; const Particle& leadingPhoton = photons[0]; if (leadingPhoton.Et() < 15.0*GeV) vetoEvent; if (leadingPhoton.abseta() > 2.37) vetoEvent; // check photon isolation double coneEnergy(0.0); Particles fs = apply(event, "isolatedFS").particles(); for (const Particle& p : fs) { if ( deltaR(leadingPhoton, p) < 0.4 ) coneEnergy += p.E(); } if (coneEnergy / leadingPhoton.E() >= 0.5 ) vetoEvent; // retrieve Z boson candidate const ZFinder& zf = apply(event, "ZF"); if ( zf.bosons().size() != 1 ) vetoEvent; // only one Z boson candidate const Particle& Zboson = zf.boson(); if ( !(Zboson.mass() > 40.0*GeV) ) vetoEvent; // check charge of constituent leptons const Particles& leptons = zf.constituents(); if (leptons.size() != 2 || leptons[0].charge() * leptons[1].charge() > 0.) vetoEvent; // check photon-lepton overlap for (const Particle& p : leptons) { if ( !(p.pT() > 25.0*GeV && p.abseta() < 2.47 && deltaR(leadingPhoton, p) > 0.7) ) vetoEvent; } // count jets const FastJets& jetfs = apply(event, "Jets"); Jets jets = jetfs.jets(cmpMomByEt); int goodJets = 0; for (const Jet& j : jets) { if ( !(j.Et() > 30.0*GeV) ) break; if ( (j.abseta() < 4.4) && \ (deltaR(leadingPhoton, j) > 0.3) && \ (deltaR(leptons[0], j) > 0.3) && \ (deltaR(leptons[1], j) > 0.3) ) ++goodJets; } double Njets = double(goodJets) + 0.5; double photonEt = leadingPhoton.Et()*GeV; double mZgamma = (Zboson.momentum() + leadingPhoton.momentum()).mass() * GeV; _hist_EgammaT_incl->fill(photonEt); _hist_Njet_EgammaT15->fill(Njets); if ( !goodJets ) _hist_EgammaT_excl->fill(photonEt); if (photonEt >= 40.0*GeV) { _hist_mZgamma->fill(mZgamma); if (photonEt >= 60.0*GeV) _hist_Njet_EgammaT60->fill(Njets); } } /// Normalise histograms etc., after the run void finalize() { const double xs_fb = crossSection()/femtobarn; const double sumw = sumOfWeights(); const double sf = xs_fb / sumw; scale(_hist_EgammaT_excl, sf); scale(_hist_EgammaT_incl, sf); normalize(_hist_Njet_EgammaT15); normalize(_hist_Njet_EgammaT60); normalize(_hist_mZgamma); } //@} protected: size_t _mode; private: /// @name Histograms //@{ Histo1DPtr _hist_EgammaT_incl; Histo1DPtr _hist_EgammaT_excl; Histo1DPtr _hist_Njet_EgammaT15; Histo1DPtr _hist_Njet_EgammaT60; Histo1DPtr _hist_mZgamma; //@} }; class ATLAS_2013_I1217863_Z_EL : public ATLAS_2013_I1217863_Z { public: ATLAS_2013_I1217863_Z_EL() : ATLAS_2013_I1217863_Z("ATLAS_2013_I1217863_Z_EL") { _mode = 2; } }; class ATLAS_2013_I1217863_Z_MU : public ATLAS_2013_I1217863_Z { public: ATLAS_2013_I1217863_Z_MU() : ATLAS_2013_I1217863_Z("ATLAS_2013_I1217863_Z_MU") { _mode = 3; } }; DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217863_Z); DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217863_Z_EL); DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217863_Z_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2013_I1217867.cc b/analyses/pluginATLAS/ATLAS_2013_I1217867.cc --- a/analyses/pluginATLAS/ATLAS_2013_I1217867.cc +++ b/analyses/pluginATLAS/ATLAS_2013_I1217867.cc @@ -1,153 +1,152 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { class ATLAS_2013_I1217867 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor ATLAS_2013_I1217867() : Analysis("ATLAS_2013_I1217867") { m_njet = 4; _h_dI.resize(2, std::vector(m_njet)); _h_dI_ratio.resize(2, std::vector(m_njet-1)); } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { // Initialise projections FinalState fs(Cuts::abseta < 5.0); IdentifiedFinalState bareElectrons(fs); bareElectrons.acceptIdPair(PID::ELECTRON); Cut cuts = (Cuts::absetaIn(0, 1.37) || Cuts::absetaIn(1.52, 2.47)) && Cuts::pT > 20*GeV; DressedLeptons electronClusters(fs, bareElectrons, 0.1, cuts); declare(electronClusters, "electronClusters"); IdentifiedFinalState bareMuons(fs); bareMuons.acceptIdPair(PID::MUON); Cut mucuts = Cuts::abseta < 2.4 && Cuts::pT > 20*GeV; DressedLeptons muonClusters(fs, bareMuons, 0.1, mucuts); declare(muonClusters, "muonClusters"); IdentifiedFinalState neutrinos(Cuts::pT > 25*GeV); neutrinos.acceptNeutrinos(); declare(neutrinos, "neutrinos"); VetoedFinalState jetFS(fs); jetFS.addVetoOnThisFinalState(electronClusters); jetFS.addVetoOnThisFinalState(muonClusters); jetFS.addVetoOnThisFinalState(neutrinos); - FastJets jetpro(jetFS, FastJets::KT, 0.6); - jetpro.useInvisibles(true); + FastJets jetpro(jetFS, FastJets::KT, 0.6, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jetpro, "jets"); // Book histograms for (size_t flav = 0; flav < 2; ++flav) { for (size_t i = 0; i < m_njet; ++i) book(_h_dI[flav][i], i+1, 1, flav+1); for (size_t i = 0; i < m_njet-1; ++i) book(_h_dI_ratio[flav][i], 4+i+1, 1, flav+1); } } /// Perform the per-event analysis void analyze(const Event& e) { const DressedLeptons& electronClusters = apply(e, "electronClusters"); const DressedLeptons& muonClusters = apply(e, "muonClusters"); int ne = electronClusters.dressedLeptons().size(); int nmu = muonClusters.dressedLeptons().size(); FourMomentum lepton; size_t flav = 2; if (ne==1) { lepton=electronClusters.dressedLeptons()[0].momentum(); flav = 0; if (nmu > 0) vetoEvent; } else if (nmu == 1) { lepton=muonClusters.dressedLeptons()[0].momentum(); flav = 1; if (ne > 0) vetoEvent; } else { vetoEvent; } const Particles& neutrinos = apply(e, "neutrinos").particlesByPt(); if (neutrinos.size() < 1) vetoEvent; FourMomentum neutrino = neutrinos[0].momentum(); double mtW=sqrt(2.0*lepton.pT()*neutrino.pT()*(1-cos(lepton.phi()-neutrino.phi()))); if (mtW<40.0*GeV) vetoEvent; const shared_ptr seq = apply(e, "jets").clusterSeq(); if (seq) { for (size_t i = 0; i < min(m_njet,(size_t)seq->n_particles()); ++i) { double d_ij = sqrt(seq->exclusive_dmerge_max(i)); _h_dI[flav][i]->fill(d_ij); if (i20.0*GeV) { double d_ijplus1 = sqrt(seq->exclusive_dmerge_max(i+1)); _h_dI_ratio[flav][i]->fill(d_ijplus1/d_ij); } } } } } /// Normalise histograms etc., after the run void finalize() { for (size_t flav = 0; flav < 2; ++flav) { for (size_t i = 0; i < m_njet; ++i) { normalize(_h_dI[flav][i], 1.0, false); if (i < m_njet-1) normalize(_h_dI_ratio[flav][i], 1.0, false); } } } //@} private: /// @name Histograms //@{ vector< vector > _h_dI; vector< vector > _h_dI_ratio; //@} size_t m_njet; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2013_I1217867); } diff --git a/analyses/pluginATLAS/ATLAS_2013_I1230812.cc b/analyses/pluginATLAS/ATLAS_2013_I1230812.cc --- a/analyses/pluginATLAS/ATLAS_2013_I1230812.cc +++ b/analyses/pluginATLAS/ATLAS_2013_I1230812.cc @@ -1,310 +1,309 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" namespace Rivet { /// Z + jets in pp at 7 TeV (combined channel / base class) /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes class ATLAS_2013_I1230812 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor ATLAS_2013_I1230812(string name="ATLAS_2013_I1230812") : Analysis(name) { // This class uses the combined e+mu mode _mode = 1; } //@} /// Book histograms and initialise projections before the run void init() { // Determine the e/mu decay channels used (NB Prompt leptons only). /// @todo Note that Zs are accepted with any rapidity: the cuts are on the e/mu: is this correct? Cut pt20 = Cuts::pT >= 20.0*GeV; if (_mode == 1) { // Combined ZFinder zfinder(FinalState(Cuts::abseta < 2.5), pt20, PID::ELECTRON, 66*GeV, 116*GeV); declare(zfinder, "zfinder"); } else if (_mode == 2) { // Electron const Cut eta_e = Cuts::abseta < 1.37 || Cuts::absetaIn(1.52, 2.47); ZFinder zfinder(FinalState(eta_e), pt20, PID::ELECTRON, 66*GeV, 116*GeV); declare(zfinder, "zfinder"); } else if (_mode == 3) { // Muon ZFinder zfinder(FinalState(Cuts::abseta < 2.4), pt20, PID::MUON, 66*GeV, 116*GeV); declare(zfinder, "zfinder"); } else { MSG_ERROR("Unknown decay channel mode!!!"); } // Define veto FS in order to prevent Z-decay products entering the jet algorithm VetoedFinalState had_fs; had_fs.addVetoOnThisFinalState(getProjection("zfinder")); - FastJets jets(had_fs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(had_fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); book(_h_njet_incl , 1, 1, _mode); book(_h_njet_incl_ratio , 2, 1, _mode, true); book(_h_njet_excl , 3, 1, _mode); book(_h_njet_excl_ratio , 4, 1, _mode, true); book(_h_njet_excl_pt150 , 5, 1, _mode); book(_h_njet_excl_pt150_ratio , 6, 1, _mode, true); book(_h_njet_excl_vbf , 7, 1, _mode); book(_h_njet_excl_vbf_ratio , 8, 1, _mode, true); book(_h_ptlead , 9, 1, _mode); book(_h_ptseclead , 10, 1, _mode); book(_h_ptthirdlead , 11, 1, _mode); book(_h_ptfourthlead , 12, 1, _mode); book(_h_ptlead_excl , 13, 1, _mode); book(_h_pt_ratio , 14, 1, _mode); book(_h_pt_z , 15, 1, _mode); book(_h_pt_z_excl , 16, 1, _mode); book(_h_ylead , 17, 1, _mode); book(_h_yseclead , 18, 1, _mode); book(_h_ythirdlead , 19, 1, _mode); book(_h_yfourthlead , 20, 1, _mode); book(_h_deltay , 21, 1, _mode); book(_h_mass , 22, 1, _mode); book(_h_deltaphi , 23, 1, _mode); book(_h_deltaR , 24, 1, _mode); book(_h_ptthirdlead_vbf , 25, 1, _mode); book(_h_ythirdlead_vbf , 26, 1, _mode); book(_h_ht , 27, 1, _mode); book(_h_st , 28, 1, _mode); } /// Perform the per-event analysis void analyze(const Event& event) { const ZFinder& zfinder = apply(event, "zfinder"); MSG_DEBUG("# Z constituents = " << zfinder.constituents().size()); if (zfinder.constituents().size() != 2) vetoEvent; FourMomentum z = zfinder.boson().momentum(); FourMomentum lp = zfinder.constituents()[0].momentum(); FourMomentum lm = zfinder.constituents()[1].momentum(); if (deltaR(lp, lm) < 0.2) vetoEvent; Jets jets = apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 4.4); ifilter_discard(jets, deltaRLess(lp, 0.5)); ifilter_discard(jets, deltaRLess(lm, 0.5)); // Fill jet multiplicities for (size_t ijet = 0; ijet <= jets.size(); ++ijet) { _h_njet_incl->fill(ijet); } _h_njet_excl->fill(jets.size()); // Require at least one jet if (jets.size() >= 1) { // Leading jet histos const double ptlead = jets[0].pT()/GeV; const double yabslead = fabs(jets[0].rapidity()); const double ptz = z.pT()/GeV; _h_ptlead->fill(ptlead); _h_ylead ->fill(yabslead); _h_pt_z ->fill(ptz); // Fill jet multiplicities if (ptlead > 150) _h_njet_excl_pt150->fill(jets.size()); // Loop over selected jets, fill inclusive distributions double st = 0; double ht = lp.pT()/GeV + lm.pT()/GeV; for (size_t ijet = 0; ijet < jets.size(); ++ijet) { ht += jets[ijet].pT()/GeV; st += jets[ijet].pT()/GeV; } _h_ht->fill(ht); _h_st->fill(st); // Require exactly one jet if (jets.size() == 1) { _h_ptlead_excl->fill(ptlead); _h_pt_z_excl ->fill(ptz); } } // Require at least two jets if (jets.size() >= 2) { // Second jet histos const double ptlead = jets[0].pT()/GeV; const double pt2ndlead = jets[1].pT()/GeV; const double ptratio = pt2ndlead/ptlead; const double yabs2ndlead = fabs(jets[1].rapidity()); _h_ptseclead->fill(pt2ndlead); _h_yseclead->fill( yabs2ndlead); _h_pt_ratio->fill( ptratio); // Dijet histos const double deltaphi = fabs(deltaPhi(jets[1], jets[0])); const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ; const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY)); const double mass = (jets[0].momentum() + jets[1].momentum()).mass()/GeV; _h_mass->fill( mass); _h_deltay->fill( deltarap); _h_deltaphi->fill(deltaphi); _h_deltaR->fill( deltar); if (mass > 350 && deltarap > 3) _h_njet_excl_vbf->fill(jets.size()); } // Require at least three jets if (jets.size() >= 3) { // Third jet histos const double pt3rdlead = jets[2].pT()/GeV; const double yabs3rdlead = fabs(jets[2].rapidity()); _h_ptthirdlead->fill(pt3rdlead); _h_ythirdlead->fill( yabs3rdlead); //Histos after VBF preselection const double deltarap = fabs(jets[0].rapidity() - jets[1].rapidity()) ; const double mass = (jets[0].momentum() + jets[1].momentum()).mass(); if (mass > 350 && deltarap > 3) { _h_ptthirdlead_vbf->fill(pt3rdlead); _h_ythirdlead_vbf->fill( yabs3rdlead); } } // Require at least four jets if (jets.size() >= 4) { // Fourth jet histos const double pt4thlead = jets[3].pT()/GeV; const double yabs4thlead = fabs(jets[3].rapidity()); _h_ptfourthlead->fill(pt4thlead); _h_yfourthlead->fill( yabs4thlead); } } /// @name Ratio calculator util functions //@{ /// Calculate the efficiency error, being careful about div-by-zero double err_incl(const HistoBin1D &M, const HistoBin1D &N, bool hasWeights) { double r = safediv(M.sumW(), N.sumW()); if (hasWeights) { // use F. James's approximation for weighted events return sqrt( safediv((1 - 2 * r) * M.sumW2() + r * r * N.sumW2(), N.sumW() * N.sumW()) ); } return sqrt( safediv(r * (1 - r), N.sumW()) ); } /// Calculate the ratio error, being careful about div-by-zero double err_excl(const HistoBin1D &A, const HistoBin1D &B) { double r = safediv(A.sumW(), B.sumW()); double dAsquared = safediv(A.sumW2(), A.sumW() * A.sumW()); // squared relative error of A double dBsquared = safediv(B.sumW2(), B.sumW() * B.sumW()); // squared relative error of B return r * sqrt(dAsquared + dBsquared); } //@} void finalize() { bool hasWeights = _h_njet_incl->effNumEntries() != _h_njet_incl->numEntries(); for (size_t i = 0; i < 6; ++i) { _h_njet_incl_ratio->point(i).setY(safediv(_h_njet_incl->bin(i + 1).sumW(), _h_njet_incl->bin(i).sumW()), err_incl(_h_njet_incl->bin(i + 1), _h_njet_incl->bin(i), hasWeights)); _h_njet_excl_ratio->point(i).setY(safediv(_h_njet_excl->bin(i + 1).sumW(), _h_njet_excl->bin(i).sumW()), err_excl(_h_njet_excl->bin(i + 1), _h_njet_excl->bin(i))); if (i >= 1) { _h_njet_excl_pt150_ratio->point(i - 1).setY(safediv(_h_njet_excl_pt150->bin(i).sumW(), _h_njet_excl_pt150->bin(i - 1).sumW()), err_excl(_h_njet_excl_pt150->bin(i), _h_njet_excl_pt150->bin(i - 1))); if (i >= 2) { _h_njet_excl_vbf_ratio->point(i - 2).setY(safediv(_h_njet_excl_vbf->bin(i).sumW(), _h_njet_excl_vbf->bin(i - 1).sumW()), err_excl(_h_njet_excl_vbf->bin(i), _h_njet_excl_vbf->bin(i - 1))); } } } const double xs = crossSectionPerEvent()/picobarn; scale({_h_njet_incl, _h_njet_excl, _h_njet_excl_pt150, _h_njet_excl_vbf}, xs); scale({_h_ptlead, _h_ptseclead, _h_ptthirdlead, _h_ptfourthlead, _h_ptlead_excl}, xs); scale({_h_pt_ratio, _h_pt_z, _h_pt_z_excl}, xs); scale({_h_ylead, _h_yseclead, _h_ythirdlead, _h_yfourthlead}, xs); scale({_h_deltay, _h_mass, _h_deltaphi, _h_deltaR}, xs); scale({_h_ptthirdlead_vbf, _h_ythirdlead_vbf}, xs); scale({_h_ht, _h_st}, xs); } //@} protected: size_t _mode; private: Scatter2DPtr _h_njet_incl_ratio; Scatter2DPtr _h_njet_excl_ratio; Scatter2DPtr _h_njet_excl_pt150_ratio; Scatter2DPtr _h_njet_excl_vbf_ratio; Histo1DPtr _h_njet_incl; Histo1DPtr _h_njet_excl; Histo1DPtr _h_njet_excl_pt150; Histo1DPtr _h_njet_excl_vbf; Histo1DPtr _h_ptlead; Histo1DPtr _h_ptseclead; Histo1DPtr _h_ptthirdlead; Histo1DPtr _h_ptfourthlead; Histo1DPtr _h_ptlead_excl; Histo1DPtr _h_pt_ratio; Histo1DPtr _h_pt_z; Histo1DPtr _h_pt_z_excl; Histo1DPtr _h_ylead; Histo1DPtr _h_yseclead; Histo1DPtr _h_ythirdlead; Histo1DPtr _h_yfourthlead; Histo1DPtr _h_deltay; Histo1DPtr _h_mass; Histo1DPtr _h_deltaphi; Histo1DPtr _h_deltaR; Histo1DPtr _h_ptthirdlead_vbf; Histo1DPtr _h_ythirdlead_vbf; Histo1DPtr _h_ht; Histo1DPtr _h_st; }; class ATLAS_2013_I1230812_EL : public ATLAS_2013_I1230812 { public: ATLAS_2013_I1230812_EL() : ATLAS_2013_I1230812("ATLAS_2013_I1230812_EL") { _mode = 2; } }; class ATLAS_2013_I1230812_MU : public ATLAS_2013_I1230812 { public: ATLAS_2013_I1230812_MU() : ATLAS_2013_I1230812("ATLAS_2013_I1230812_MU") { _mode = 3; } }; DECLARE_RIVET_PLUGIN(ATLAS_2013_I1230812); DECLARE_RIVET_PLUGIN(ATLAS_2013_I1230812_EL); DECLARE_RIVET_PLUGIN(ATLAS_2013_I1230812_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2014_I1319490.cc b/analyses/pluginATLAS/ATLAS_2014_I1319490.cc --- a/analyses/pluginATLAS/ATLAS_2014_I1319490.cc +++ b/analyses/pluginATLAS/ATLAS_2014_I1319490.cc @@ -1,232 +1,231 @@ #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { class ATLAS_2014_I1319490 : public Analysis { public: ATLAS_2014_I1319490(string name = "ATLAS_2014_I1319490") : Analysis(name) { _mode = 0; // using electron channel for combined data by default } // Book histograms and initialise projections before the run void init() { FinalState fs; Cut cuts; if (_mode == 2) { // muon channel cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.4, 2.4); } else if (_mode) { // electron channel cuts = (Cuts::pT > 25.0*GeV) & ( Cuts::etaIn(-2.47, -1.52) | Cuts::etaIn(-1.37, 1.37) | Cuts::etaIn(1.52, 2.47) ); } else { // combined data extrapolated to common phase space cuts = (Cuts::pT > 25.0*GeV) & Cuts::etaIn(-2.5, 2.5); } // bosons WFinder wfinder(fs, cuts, _mode > 1? PID::MUON : PID::ELECTRON, 40.0*GeV, DBL_MAX, 0.0*GeV, 0.1, WFinder::ChargedLeptons::PROMPT, WFinder::ClusterPhotons::NODECAY, WFinder::AddPhotons::NO, WFinder::MassWindow::MT); declare(wfinder, "WF"); // jets VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("WF")); - FastJets jets(jet_fs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(jet_fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "Jets"); // book histograms book(histos["h_N_incl"] ,1,1,_mode+1); book(histos["h_N"] ,4,1,_mode+1); book(histos["h_pt_jet1_1jet"] ,5,1,_mode+1); book(histos["h_pt_jet1_1jet_excl"] ,6,1,_mode+1); book(histos["h_pt_jet1_2jet"] ,7,1,_mode+1); book(histos["h_pt_jet1_3jet"] ,8,1,_mode+1); book(histos["h_pt_jet2_2jet"] ,9,1,_mode+1); book(histos["h_pt_jet3_3jet"] ,10,1,_mode+1); book(histos["h_pt_jet4_4jet"] ,11,1,_mode+1); book(histos["h_pt_jet5_5jet"] ,12,1,_mode+1); book(histos["h_y_jet1_1jet"] ,13,1,_mode+1); book(histos["h_y_jet2_2jet"] ,14,1,_mode+1); book(histos["h_HT_1jet"] ,15,1,_mode+1); book(histos["h_HT_1jet_excl"] ,16,1,_mode+1); book(histos["h_HT_2jet"] ,17,1,_mode+1); book(histos["h_HT_2jet_excl"] ,18,1,_mode+1); book(histos["h_HT_3jet"] ,19,1,_mode+1); book(histos["h_HT_3jet_excl"] ,20,1,_mode+1); book(histos["h_HT_4jet"] ,21,1,_mode+1); book(histos["h_HT_5jet"] ,22,1,_mode+1); book(histos["h_deltaPhi_jet12"] ,23,1,_mode+1); book(histos["h_deltaRap_jet12"] ,24,1,_mode+1); book(histos["h_deltaR_jet12"] ,25,1,_mode+1); book(histos["h_M_Jet12_2jet"] ,26,1,_mode+1); book(histos["h_y_jet3_3jet"] ,27,1,_mode+1); book(histos["h_y_jet4_4jet"] ,28,1,_mode+1); book(histos["h_y_jet5_5jet"] ,29,1,_mode+1); book(histos["h_ST_1jet"] ,30,1,_mode+1); book(histos["h_ST_2jet"] ,31,1,_mode+1); book(histos["h_ST_2jet_excl"] ,32,1,_mode+1); book(histos["h_ST_3jet"] ,33,1,_mode+1); book(histos["h_ST_3jet_excl"] ,34,1,_mode+1); book(histos["h_ST_4jet"] ,35,1,_mode+1); book(histos["h_ST_5jet"] ,36,1,_mode+1); } void fillPlots(const Particle& lepton, const double& missET, Jets& all_jets) { // do jet-lepton overlap removal Jets jets; double ST = 0.0; // scalar pT sum of all selected jets for (const Jet &j : all_jets) { if (deltaR(j, lepton) > 0.5) { jets += j; ST += j.pT() / GeV; } } const size_t njets = jets.size(); const double HT = ST + lepton.pT() / GeV + missET; histos["h_N"]->fill(njets + 0.5); for (size_t i = 0; i <= njets; ++i) { histos["h_N_incl"]->fill(i + 0.5); } if (njets) { const double pT1 = jets[0].pT() / GeV; const double rap1 = jets[0].absrap(); histos["h_pt_jet1_1jet" ]->fill(pT1); histos["h_y_jet1_1jet"]->fill(rap1); histos["h_HT_1jet"]->fill(HT); histos["h_ST_1jet"]->fill(ST); if (njets == 1) { histos["h_pt_jet1_1jet_excl"]->fill(pT1); histos["h_HT_1jet_excl"]->fill(HT); } else { const double pT2 = jets[1].pT() / GeV; const double rap2 = jets[1].absrap(); const double dR = deltaR(jets[0], jets[1]); const double dRap = deltaRap(jets[0], jets[1]); const double dPhi = deltaPhi(jets[0], jets[1]); const double mjj = (jets[0].momentum() + jets[1].momentum()).mass() / GeV; histos["h_pt_jet1_2jet"]->fill(pT1); histos["h_pt_jet2_2jet"]->fill(pT2); histos["h_y_jet2_2jet"]->fill(rap2); histos["h_M_Jet12_2jet"]->fill(mjj); histos["h_HT_2jet"]->fill(HT); histos["h_ST_2jet"]->fill(ST); histos["h_deltaPhi_jet12"]->fill(dPhi); histos["h_deltaRap_jet12"]->fill(dRap); histos["h_deltaR_jet12"]->fill(dR); if (njets == 2) { histos["h_ST_2jet_excl"]->fill(ST); histos["h_HT_2jet_excl"]->fill(HT); } else { const double pT3 = jets[2].pT() / GeV; const double rap3 = jets[2].absrap(); histos["h_pt_jet1_3jet"]->fill(pT1); histos["h_pt_jet3_3jet"]->fill(pT3); histos["h_y_jet3_3jet"]->fill(rap3); histos["h_HT_3jet"]->fill(HT); histos["h_ST_3jet"]->fill(ST); if(njets == 3) { histos["h_ST_3jet_excl"]->fill(ST); histos["h_HT_3jet_excl"]->fill(HT); } else { const double pT4 = jets[3].pT() / GeV; const double rap4 = jets[3].absrap(); histos["h_pt_jet4_4jet"]->fill(pT4); histos["h_y_jet4_4jet"]->fill(rap4); histos["h_HT_4jet"]->fill(HT); histos["h_ST_4jet"]->fill(ST); if (njets > 4) { const double pT5 = jets[4].pT() / GeV; const double rap5 = jets[4].absrap(); histos["h_pt_jet5_5jet"]->fill(pT5); histos["h_y_jet5_5jet"]->fill(rap5); histos["h_HT_5jet"]->fill(HT); histos["h_ST_5jet"]->fill(ST); } } } } } } // Perform the per-event analysis void analyze(const Event& event) { // Retrieve boson candidate const WFinder& wf = apply(event, "WF"); if (wf.empty()) vetoEvent; // Retrieve jets const JetAlg& jetfs = apply(event, "Jets"); Jets all_jets = jetfs.jetsByPt(Cuts::pT > 30.0*GeV && Cuts::absrap < 4.4); const Particles& leptons = wf.constituentLeptons(); const double missET = wf.constituentNeutrino().pT() / GeV; if (leptons.size() == 1 && missET > 25.0 && wf.mT() > 40.0*GeV) { const Particle& lep = leptons[0]; fillPlots(lep, missET, all_jets); } } void finalize() { const double scalefactor(crossSection() / sumOfWeights()); /// @todo Update to use C++11 range-for for (map::iterator hit = histos.begin(); hit != histos.end(); ++hit) { scale(hit->second, scalefactor); } } protected: size_t _mode; private: map histos; }; class ATLAS_2014_I1319490_EL : public ATLAS_2014_I1319490 { public: ATLAS_2014_I1319490_EL() : ATLAS_2014_I1319490("ATLAS_2014_I1319490_EL") { _mode = 1; } }; class ATLAS_2014_I1319490_MU : public ATLAS_2014_I1319490 { public: ATLAS_2014_I1319490_MU() : ATLAS_2014_I1319490("ATLAS_2014_I1319490_MU") { _mode = 2; } }; // The hooks for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490_EL); DECLARE_RIVET_PLUGIN(ATLAS_2014_I1319490_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc b/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc --- a/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc +++ b/analyses/pluginATLAS/ATLAS_2015_CONF_2015_041.cc @@ -1,160 +1,159 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { /// Z + jets in pp at 13 TeV /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes class ATLAS_2015_CONF_2015_041 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor ATLAS_2015_CONF_2015_041(string name="ATLAS_2015_CONF_2015_041") : Analysis(name), _weights(5) { // This class uses the combined e+mu mode _mode = 0; } //@} /// Book histograms and initialise projections before the run void init() { const FinalState fs; Cut cuts = (Cuts::pT > 25*GeV) & (Cuts::abseta < 2.5); ZFinder zfinder(fs, cuts, _mode? PID::MUON : PID::ELECTRON, 66*GeV, 116*GeV); declare(zfinder, "zfinder"); // Define veto FS in order to prevent Z-decay products entering the jet algorithm VetoedFinalState had_fs; had_fs.addVetoOnThisFinalState(zfinder); - FastJets jets(had_fs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(had_fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); // individual channels book(_hNjets ,1, 1, _mode + 1); book(_hNjetsRatio ,2, 1, _mode + 1, true); // combination book(_hNjets_comb ,1, 2, _mode + 1); book(_hNjetsRatio_comb ,2, 2, _mode + 1, true); for (size_t i = 0; i < 5; i++) book(_weights[i], "_weights" + to_str(i)); } /// Perform the per-event analysis void analyze(const Event& event) { const ZFinder& zfinder = apply(event, "zfinder"); const Particles& leptons = zfinder.constituents(); if (leptons.size() != 2) vetoEvent; Jets jets; for (Jet j : apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 2.5)) { bool keep = true; for(const Particle& l : leptons) keep &= deltaR(j, l) > 0.4; if (keep) jets += j; } size_t njets = jets.size(); for(size_t i = 0; i <= njets; ++i) { _hNjets->fill(i + 0.5); _hNjets_comb->fill(i + 0.5); } for (size_t i = 0; i < 5; ++i) { if (njets >= i) _weights[i]->fill(); } } /// @name Ratio calculator util functions //@{ /// Calculate the ratio, being careful about div-by-zero double ratio(double a, double b) { return (b != 0) ? a/b : 0; } /// Calculate the ratio error, being careful about div-by-zero double ratio_err(double a, double b) { return (b != 0) ? sqrt(a/b*(1-a/b)/b) : 0; } //@} void finalize() { for (size_t i = 0; i < 4; ++i) { double n = _hNjets->bin(i + 1).sumW(); double dN = _hNjets->bin(i + 1).sumW2(); double d = _hNjets->bin(i).sumW(); double dD = _hNjets->bin(i).sumW2(); double r = safediv(n, d); double e = sqrt( safediv(r * (1 - r), d) ); if ( _hNjets->effNumEntries() != _hNjets->numEntries() ) { // use F. James's approximation for weighted events: e = sqrt( safediv((1 - 2 * r) * dN + r * r * dD, d * d) ); } _hNjetsRatio->point(i).setY(r, e); _hNjetsRatio_comb->point(i).setY(r, e); } scale(_hNjets, crossSectionPerEvent() ); scale(_hNjets_comb, crossSectionPerEvent() ); } //@} protected: size_t _mode; private: vector _weights; Scatter2DPtr _hNjetsRatio, _hNjetsRatio_comb; Histo1DPtr _hNjets, _hNjets_comb; }; class ATLAS_2015_CONF_2015_041_EL : public ATLAS_2015_CONF_2015_041 { public: ATLAS_2015_CONF_2015_041_EL() : ATLAS_2015_CONF_2015_041("ATLAS_2015_CONF_2015_041_EL") { _mode = 0; } }; class ATLAS_2015_CONF_2015_041_MU : public ATLAS_2015_CONF_2015_041 { public: ATLAS_2015_CONF_2015_041_MU() : ATLAS_2015_CONF_2015_041("ATLAS_2015_CONF_2015_041_MU") { _mode = 1; } }; DECLARE_RIVET_PLUGIN(ATLAS_2015_CONF_2015_041); DECLARE_RIVET_PLUGIN(ATLAS_2015_CONF_2015_041_EL); DECLARE_RIVET_PLUGIN(ATLAS_2015_CONF_2015_041_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2015_I1404878.cc b/analyses/pluginATLAS/ATLAS_2015_I1404878.cc --- a/analyses/pluginATLAS/ATLAS_2015_I1404878.cc +++ b/analyses/pluginATLAS/ATLAS_2015_I1404878.cc @@ -1,260 +1,259 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VisibleFinalState.hh" namespace Rivet { class ATLAS_2015_I1404878 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1404878); void init() { // Eta ranges Cut eta_full = (Cuts::abseta < 4.2) & (Cuts::pT >= 1.0*MeV); Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV); // All final state particles FinalState fs(eta_full); // Get photons to dress leptons IdentifiedFinalState photons(fs); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); electrons.acceptTauDecays(true); declare(electrons, "electrons"); DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true); declare(dressedelectrons, "dressedelectrons"); DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true); declare(ewdressedelectrons, "ewdressedelectrons"); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); declare(muons, "muons"); DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true); declare(dressedmuons, "dressedmuons"); DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true); declare(ewdressedmuons, "ewdressedmuons"); // Projection to find neutrinos IdentifiedFinalState nu_id; nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); // get MET from generic invisibles VetoedFinalState inv_fs(fs); inv_fs.addVetoOnThisFinalState(VisibleFinalState(fs)); declare(inv_fs, "InvisibleFS"); // Jet clustering. VetoedFinalState vfs; vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); - FastJets jets(vfs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); // Histogram booking book(_h["massttbar"] , 1, 1, 1); book(_h["massttbar_norm"] , 2, 1, 1); book(_h["ptttbar"] , 3, 1, 1); book(_h["ptttbar_norm"] , 4, 1, 1); book(_h["absrapttbar"] , 5, 1, 1); book(_h["absrapttbar_norm"] , 6, 1, 1); book(_h["ptpseudotophadron"] , 7, 1, 1); book(_h["ptpseudotophadron_norm"] , 8, 1, 1); book(_h["absrappseudotophadron"] , 9, 1, 1); book(_h["absrappseudotophadron_norm"] ,10, 1, 1); book(_h["absPout"] ,11, 1, 1); book(_h["absPout_norm"] ,12, 1, 1); book(_h["dPhittbar"] ,13, 1, 1); book(_h["dPhittbar_norm"] ,14, 1, 1); book(_h["HTttbar"] ,15, 1, 1); book(_h["HTttbar_norm"] ,16, 1, 1); book(_h["Yboost"] ,17, 1, 1); book(_h["Yboost_norm"] ,18, 1, 1); book(_h["chittbar"] ,19, 1, 1); book(_h["chittbar_norm"] ,20, 1, 1); book(_h["RWt"] ,21, 1, 1); book(_h["RWt_norm"] ,22, 1, 1); } void analyze(const Event& event) { // Get the selected objects, using the projections. vector electrons = applyProjection(event, "dressedelectrons").dressedLeptons(); vector muons = applyProjection(event, "dressedmuons").dressedLeptons(); const Jets& jets = applyProjection(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); const FinalState& ifs = applyProjection(event, "InvisibleFS"); // Calculate MET FourMomentum met; for (const Particle& p : ifs.particles()) met += p.momentum(); // Count the number of b-tags Jets bjets, lightjets; for (const Jet& jet : jets){ bool b_tagged = jet.bTags(Cuts::pT > 5*GeV).size(); if ( b_tagged && bjets.size() < 2 ) bjets += jet; else lightjets += jet; } bool single_electron = (electrons.size() == 1) && (muons.empty()); bool single_muon = (muons.size() == 1) && (electrons.empty()); DressedLepton* lepton = nullptr; if (single_electron) lepton = &electrons[0]; else if (single_muon) lepton = &muons[0]; if(!single_electron && !single_muon) vetoEvent; if (jets.size() < 4 || bjets.size() < 2) vetoEvent; FourMomentum pbjet1; // Momentum of bjet1 FourMomentum pbjet2; // Momentum of bjet2 if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) { pbjet1 = bjets[0].momentum(); pbjet2 = bjets[1].momentum(); } else { pbjet1 = bjets[1].momentum(); pbjet2 = bjets[0].momentum(); } double bestWmass = 1000.0*TeV; double mWPDG = 80.399*GeV; int Wj1index = -1, Wj2index = -1; for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) { for (unsigned int j = i + 1; j < lightjets.size(); ++j) { double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass(); if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) { bestWmass = wmass; Wj1index = i; Wj2index = j; } } } // Compute hadronic W boson FourMomentum pWhadron = lightjets[Wj1index].momentum() + lightjets[Wj2index].momentum(); double pz = _computeneutrinoz(lepton->momentum(), met); FourMomentum ppseudoneutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz); // Compute leptonic, hadronic, combined pseudo-top FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1; FourMomentum ppseudotophadron = pbjet2 + pWhadron; FourMomentum pttbar = ppseudotoplepton + ppseudotophadron; Vector3 z_versor(0,0,1); Vector3 vpseudotophadron = ppseudotophadron.vector3(); Vector3 vpseudotoplepton = ppseudotoplepton.vector3(); // Observables double ystar = 0.5 * deltaRap(ppseudotophadron, ppseudotoplepton); double chi_ttbar = exp(2 * fabs(ystar)); double deltaPhi_ttbar = deltaPhi(ppseudotoplepton,ppseudotophadron); double HT_ttbar = ppseudotophadron.pt() + ppseudotoplepton.pt(); double Yboost = 0.5 * fabs(ppseudotophadron.rapidity() + ppseudotoplepton.rapidity()); double R_Wt = pWhadron.pt() / ppseudotophadron.pt(); double absPout = fabs(vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod()))); // absolute cross sections _h["ptpseudotophadron"]->fill( ppseudotophadron.pt()); //pT of pseudo top hadron _h["ptttbar"]->fill( pttbar.pt()); //fill pT of ttbar in combined channel _h["absrappseudotophadron"]->fill(ppseudotophadron.absrap()); _h["absrapttbar"]->fill( pttbar.absrap()); _h["massttbar"]->fill( pttbar.mass()); _h["absPout"]->fill( absPout); _h["chittbar"]->fill( chi_ttbar); _h["dPhittbar"]->fill( deltaPhi_ttbar); _h["HTttbar"]->fill( HT_ttbar); _h["Yboost"]->fill( Yboost); _h["RWt"]->fill( R_Wt); // normalised cross sections _h["ptpseudotophadron_norm"]->fill( ppseudotophadron.pt()); //pT of pseudo top hadron _h["ptttbar_norm"]->fill( pttbar.pt()); //fill pT of ttbar in combined channel _h["absrappseudotophadron_norm"]->fill(ppseudotophadron.absrap()); _h["absrapttbar_norm"]->fill( pttbar.absrap()); _h["massttbar_norm"]->fill( pttbar.mass()); _h["absPout_norm"]->fill( absPout); _h["chittbar_norm"]->fill( chi_ttbar); _h["dPhittbar_norm"]->fill( deltaPhi_ttbar); _h["HTttbar_norm"]->fill( HT_ttbar); _h["Yboost_norm"]->fill( Yboost); _h["RWt_norm"]->fill( R_Wt); } void finalize() { // Normalize to cross-section const double sf = crossSection() / sumOfWeights(); for (auto& k_h : _h) { scale(k_h.second, sf); if (k_h.first.find("_norm") != string::npos) normalize(k_h.second); } } private: // Compute z component of neutrino momentum given lepton and met double _computeneutrinoz(const FourMomentum& lepton, FourMomentum& met) const { double m_W = 80.399; // in GeV, given in the paper double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py()); double a = sqr ( lepton.E() )- sqr ( lepton.pz() ); double b = -2*k*lepton.pz(); double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k ); double discriminant = sqr(b) - 4 * a * c; double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns double pzneutrino; if (discriminant < 0) { // if the discriminant is negative: pzneutrino = - b / (2 * a); } else { // if the discriminant is positive, take the soln with smallest absolute value pzneutrino = (fabs(quad[0]) < fabs(quad[1])) ? quad[0] : quad[1]; } return pzneutrino; } /// @todo Replace with central version double _mT(const FourMomentum &l, FourMomentum &nu) const { return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) ); } /// @name Objects that are used by the event selection decisions map _h; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2015_I1404878); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1495243.cc b/analyses/pluginATLAS/ATLAS_2017_I1495243.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1495243.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1495243.cc @@ -1,221 +1,220 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" namespace Rivet { /// @brief $t\bar{t}$ + jets at 13 TeV class ATLAS_2017_I1495243 : public Analysis { public: DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1495243); void init() { Cut eta_full = Cuts::abseta < 5.0 && Cuts::pT > 1.0*MeV; Cut eta_lep = Cuts::abseta < 2.5; // Collect final state particles FinalState FS(eta_full); // Get photons to dress leptons IdentifiedFinalState photons(FS); photons.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(FS); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState electrons(el_id); electrons.acceptTauDecays(false); DressedLeptons dressedelectrons(photons, electrons, 0.1, Cuts::abseta< 2.5 && Cuts::pT > 25.0*GeV, true); declare(dressedelectrons, "electrons"); DressedLeptons fulldressedelectrons(photons, electrons, 0.1, eta_full, true); // Projection to find the muons IdentifiedFinalState mu_id(FS); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(false); DressedLeptons dressedmuons(photons, muons, 0.1, Cuts::abseta < 2.5 && Cuts::pT > 25.0*GeV, true); declare(dressedmuons, "muons"); DressedLeptons fulldressedmuons(photons, muons, 0.1, eta_full, true); // Projection to find neutrinos to exclude from jets IdentifiedFinalState nu_id; nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(false); // Jet clustering VetoedFinalState vfs; vfs.addVetoOnThisFinalState(fulldressedelectrons); vfs.addVetoOnThisFinalState(fulldressedmuons); vfs.addVetoOnThisFinalState(neutrinos); - FastJets jets(vfs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); // Book Histograms book(_h["bjet_pt"] , 5,1,1); book(_h["2bjet_pt"], 6,1,1); book(_h["ljet_pt"] , 7,1,1); for (size_t i = 0; i < 4; ++i) { book(_h["njet" + to_str(i)], i+1, 1, 1); book(_h["Q0" + to_str(i)], "_Q0" + to_str(i+ 7), refData((i>1?"d":"d0") + to_str(i+ 8) + "-x01-y01")); book(_h["MQ0" + to_str(i)], "_MQ0" + to_str(i+12), refData("d" + to_str(i+12) + "-x01-y01")); book(_h["Qsum" + to_str(i)], "_Qsum" + to_str(i+16), refData("d" + to_str(i+16) + "-x01-y01")); book(_h["MQsum" + to_str(i)], "_MQsum" + to_str(i+20), refData("d" + to_str(i+20) + "-x01-y01")); book(_s["gapFracQ0" + to_str(i)], 8+i, 1 ,1, true); book(_s["gapFracMQ0" + to_str(i)], 12+i, 1, 1, true); book(_s["gapFracQsum" + to_str(i)], 16+i, 1, 1, true); book(_s["gapFracMQsum" + to_str(i)], 20+i, 1, 1, true); } } void analyze(const Event& event) { // Get the selected objects, using the projections. Jets all_jets = apply(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); const vector electrons = filter_discard(apply(event, "electrons").dressedLeptons(), [&](const DressedLepton &e) { return any(all_jets, deltaRLess(e, 0.4)); }); const vector muons = filter_discard(apply(event, "muons").dressedLeptons(), [&](const DressedLepton &m) { return any(all_jets, deltaRLess(m, 0.4)); }); if (electrons.size() != 1 || muons.size() != 1) vetoEvent; if (electrons[0].charge() == muons[0].charge()) vetoEvent; Jets bjets, extrajets; for (Jet j : all_jets) { size_t b_tagged = j.bTags(Cuts::pT > 5*GeV).size(); if (bjets.size() < 2 && b_tagged) bjets += j; else extrajets += j; } if (bjets.size() < 2) vetoEvent; double bjetpt = bjets[0].pt(); if (bjetpt > 250*GeV) bjetpt = 275*GeV; _h["bjet_pt"]->fill(bjetpt); double b2jetpt = bjets[1].pt(); if (b2jetpt > 150*GeV) b2jetpt = 175*GeV; _h["2bjet_pt"]->fill(b2jetpt); if (extrajets.size()) { double ljetpt = extrajets[0].pt(); if (ljetpt > 250*GeV) ljetpt = 275*GeV; _h["ljet_pt"]->fill(ljetpt); } double Memubb = (electrons[0].momentum() + muons[0].momentum() + bjets[0].momentum() + bjets[1].momentum()).mass(); vector leadpt = { 0., 0., 0., 0. }, ptsum = { 0., 0., 0., 0. }; vector njetcount = { 0, 0, 0, 0 }; for (size_t i = 0; i < extrajets.size(); ++i) { double absrap = extrajets[i].absrap(), pt = extrajets[i].pT(); if (pt > 25*GeV) ++njetcount[0]; if (pt > 40*GeV) ++njetcount[1]; if (pt > 60*GeV) ++njetcount[2]; if (pt > 80*GeV) ++njetcount[3]; if (absrap < 0.8 && pt > leadpt[0]) leadpt[0] = pt; else if (absrap > 0.8 && absrap < 1.5 && pt > leadpt[1]) leadpt[1] = pt; else if (absrap > 1.5 && absrap < 2.1 && pt > leadpt[2]) leadpt[2] = pt; if (absrap < 2.1 && pt > leadpt[3]) leadpt[3] = pt; if (absrap < 0.8) ptsum[0] += pt; else if (absrap > 0.8 && absrap < 1.5) ptsum[1] += pt; else if (absrap > 1.5 && absrap < 2.1) ptsum[2] += pt; if (absrap < 2.1) ptsum[3] += pt; } for (size_t i = 0; i < 4; ++i) { size_t cutoff = i? 3 : 4; if (njetcount[i] > cutoff) njetcount[i] = cutoff; _h["njet" + to_str(i)]->fill(njetcount[i]); if (leadpt[i] > 305*GeV) leadpt[i] = 305*GeV; _h["Q0" + to_str(i)]->fill(leadpt[i]); if (ptsum[i] > 505*GeV) ptsum[i] = 505*GeV; _h["Qsum" + to_str(i)]->fill(ptsum[i]); } for (size_t i = 0; i < 4; ++i) { if (i == 0 && !(Memubb < 300*GeV)) continue; if (i == 1 && !(Memubb > 300*GeV && Memubb < 425*GeV)) continue; if (i == 2 && !(Memubb > 425*GeV && Memubb < 600*GeV)) continue; if (i == 3 && !(Memubb > 600*GeV)) continue; _h["MQ0" + to_str(i)]->fill(leadpt[3]); _h["MQsum" + to_str(i)]->fill(ptsum[3]); } } void constructGapFraction(Scatter2DPtr out, Histo1DPtr in) { bool hasWeights = in->effNumEntries() != in->numEntries(); double denW = in->sumW(); double denW2 = in->sumW2(); size_t nEnd = out->numPoints(); for (size_t i = 0; i < nEnd; ++i) { double numW = in->sumW(), numW2 = in->sumW2(); for (size_t j = i; j < nEnd; ++j) { numW -= in->bin(j).sumW(); numW2 -= in->bin(j).sumW2(); } double yval = safediv(numW, denW); double yerr = sqrt(safediv(yval * (1 - yval), denW)); if (hasWeights) { // use F. James's approximation for weighted events yerr = sqrt( safediv((1 - 2 * yval) * numW2 + yval * yval * denW2, denW * denW) ); } out->point(i).setY(yval, yerr); } } void finalize() { // Build gap fraction plots for (size_t i = 0; i < 4; ++i) { constructGapFraction(_s["gapFracQ0" + to_str(i)], _h["Q0" + to_str(i)]); constructGapFraction(_s["gapFracMQ0" + to_str(i)], _h["MQ0" + to_str(i)]); constructGapFraction(_s["gapFracQsum" + to_str(i)], _h["Qsum" + to_str(i)]); constructGapFraction(_s["gapFracMQsum" + to_str(i)], _h["MQsum" + to_str(i)]); } // Normalize to cross-section for (map::iterator hit = _h.begin(); hit != _h.end(); ++hit) { if (hit->first.find("jet") != string::npos) normalize(hit->second); } } private: /// @name Histogram helper functions map _h; map _s; }; // Declare the class as a hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1495243); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1514251.cc b/analyses/pluginATLAS/ATLAS_2017_I1514251.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1514251.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1514251.cc @@ -1,199 +1,198 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { /// Z + jets in pp at 13 TeV /// @note This base class contains a "mode" variable for combined, e, and mu channel derived classes class ATLAS_2017_I1514251 : public Analysis { public: /// Constructor ATLAS_2017_I1514251(string name="ATLAS_2017_I1514251") : Analysis(name) { // This class uses the combined e+mu mode _mode = 0; } /// Book histograms and initialise projections before the run void init() { const FinalState fs; Cut cuts = (Cuts::pT > 25*GeV) && (Cuts::abseta < 2.5); ZFinder zeefinder(fs, cuts, PID::ELECTRON, 71*GeV, 111*GeV); ZFinder zmumufinder(fs, cuts, PID::MUON, 71*GeV, 111*GeV); declare(zeefinder, "zeefinder"); declare(zmumufinder, "zmumufinder"); // Define veto FS in order to prevent Z-decay products entering the jet algorithm VetoedFinalState had_fs; had_fs.addVetoOnThisFinalState(zeefinder); had_fs.addVetoOnThisFinalState(zmumufinder); - FastJets jets(had_fs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(had_fs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); // individual channels book(_h_Njets , 1, 1, _mode + 1); book(_h_Njets_Ratio, 1, 2, _mode + 1, true); book(_h_Njets_excl , 1, 3, _mode + 1); book(_h_HT , 1, 4, _mode + 1); book(_h_leading_jet_rap , 1, 5, _mode + 1); book(_h_leading_jet_pT_eq1jet, 1, 6, _mode + 1); book(_h_leading_jet_pT , 1, 7, _mode + 1); book(_h_leading_jet_pT_2jet , 1, 8, _mode + 1); book(_h_leading_jet_pT_3jet , 1, 9, _mode + 1); book(_h_leading_jet_pT_4jet , 1, 10, _mode + 1); book(_h_jet_dphi , 1, 11, _mode + 1); book(_h_jet_mass , 1, 12, _mode + 1); } /// Perform the per-event analysis void analyze(const Event& event) { const ZFinder& zeefinder = apply(event, "zeefinder"); const ZFinder& zmumufinder = apply(event, "zmumufinder"); const Particles& zees = zeefinder.bosons(); const Particles& zmumus = zmumufinder.bosons(); //Veto Z->mumu in electron mode, and vice versa: if (_mode==1 && (zees.size()!=1 || zmumus.size() ) ) vetoEvent; else if (_mode==2 && (zees.size() || zmumus.size()!=1 ) ) vetoEvent; else if (zees.size() + zmumus.size() != 1) { // Running in combined mode, we did not find exactly one Z. Not good. MSG_DEBUG("Did not find exactly one good Z candidate"); vetoEvent; } // Find the (dressed!) leptons const Particles& leptons = zees.size() ? zeefinder.constituents() : zmumufinder.constituents(); if (leptons.size() != 2) vetoEvent; Jets jets = apply(event, "jets").jetsByPt(Cuts::pT > 30*GeV && Cuts::absrap < 2.5); bool veto = false; for(const Jet& j : jets) { for(const Particle& l : leptons) { veto |= deltaR(j, l) < 0.4; } } if (veto) vetoEvent; double HT=0; for(const Particle& l : leptons) { HT += l.pT(); } const size_t Njets = jets.size(); _h_Njets_excl->fill(Njets); for(size_t i = 0; i <= Njets; ++i) { _h_Njets->fill(i); } if (Njets < 1) vetoEvent; for(size_t i = 0; i < Njets; ++i) { HT += jets[i].pT(); } const double pT = jets[0].pT(); const double rap = jets[0].rapidity(); _h_HT->fill(HT); _h_leading_jet_rap->fill(fabs(rap)); _h_leading_jet_pT->fill(pT); if (Njets == 1) _h_leading_jet_pT_eq1jet->fill(pT); if (Njets > 1) { _h_leading_jet_pT_2jet->fill(pT); _h_jet_dphi->fill( deltaPhi(jets[0], jets[1])); _h_jet_mass->fill( (jets[0].momentum()+jets[1].momentum()).mass() ); } if (Njets > 2) _h_leading_jet_pT_3jet->fill(pT); if (Njets > 3) _h_leading_jet_pT_4jet->fill(pT); } void finalize() { for (size_t i = 0; i < _h_Njets->numBins()-2; ++i) { double n = _h_Njets->bin(i + 1).sumW(); double dN = _h_Njets->bin(i + 1).sumW2(); double d = _h_Njets->bin(i).sumW(); double dD = _h_Njets->bin(i).sumW2(); double r = safediv(n, d); double e = sqrt( safediv(r * (1 - r), d) ); if ( _h_Njets->effNumEntries() != _h_Njets->numEntries() ) { // use F. James's approximation for weighted events: e = sqrt( safediv((1 - 2 * r) * dN + r * r * dD, d * d) ); } _h_Njets_Ratio->point(i).setY(r, e); } scale(_h_Njets, crossSectionPerEvent() ); scale(_h_Njets_excl, crossSectionPerEvent() ); scale(_h_HT, crossSectionPerEvent() ); scale(_h_leading_jet_rap, crossSectionPerEvent() ); scale(_h_leading_jet_pT, crossSectionPerEvent() ); scale(_h_leading_jet_pT_eq1jet, crossSectionPerEvent() ); scale(_h_leading_jet_pT_2jet, crossSectionPerEvent() ); scale(_h_leading_jet_pT_3jet, crossSectionPerEvent() ); scale(_h_leading_jet_pT_4jet, crossSectionPerEvent() ); scale(_h_jet_dphi, crossSectionPerEvent() ); scale(_h_jet_mass, crossSectionPerEvent() ); } //@} protected: size_t _mode; private: Scatter2DPtr _h_Njets_Ratio; Histo1DPtr _h_Njets; Scatter2DPtr _h_Njets_excl_Ratio; Histo1DPtr _h_Njets_excl; Histo1DPtr _h_HT; Histo1DPtr _h_leading_jet_rap; Histo1DPtr _h_leading_jet_pT; Histo1DPtr _h_leading_jet_pT_eq1jet; Histo1DPtr _h_leading_jet_pT_2jet; Histo1DPtr _h_leading_jet_pT_3jet; Histo1DPtr _h_leading_jet_pT_4jet; Histo1DPtr _h_jet_dphi; Histo1DPtr _h_jet_mass; }; class ATLAS_2017_I1514251_EL : public ATLAS_2017_I1514251 { public: ATLAS_2017_I1514251_EL() : ATLAS_2017_I1514251("ATLAS_2017_I1514251_EL") { _mode = 1; } }; class ATLAS_2017_I1514251_MU : public ATLAS_2017_I1514251 { public: ATLAS_2017_I1514251_MU() : ATLAS_2017_I1514251("ATLAS_2017_I1514251_MU") { _mode = 2; } }; DECLARE_RIVET_PLUGIN(ATLAS_2017_I1514251); DECLARE_RIVET_PLUGIN(ATLAS_2017_I1514251_EL); DECLARE_RIVET_PLUGIN(ATLAS_2017_I1514251_MU); } diff --git a/analyses/pluginATLAS/ATLAS_2017_I1614149.cc b/analyses/pluginATLAS/ATLAS_2017_I1614149.cc --- a/analyses/pluginATLAS/ATLAS_2017_I1614149.cc +++ b/analyses/pluginATLAS/ATLAS_2017_I1614149.cc @@ -1,351 +1,350 @@ // -*- C++ -* #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/IdentifiedFinalState.hh" #include "Rivet/Projections/PromptFinalState.hh" #include "Rivet/Projections/DressedLeptons.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/MissingMomentum.hh" #include "fastjet/tools/Filter.hh" // substructure includes included in fjcontrib-1.021 (http://fastjet.hepforge.org/contrib/) #include "Rivet/Tools/Nsubjettiness/Njettiness.hh" #include "Rivet/Tools/Nsubjettiness/Nsubjettiness.hh" #include "Rivet/Tools/Nsubjettiness/NjettinessPlugin.hh" namespace Rivet { class ATLAS_2017_I1614149 : public Analysis { public: /// Constructor ///@brief: Resolved and boosted ttbar l+jets cross sections at 13 TeV DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2017_I1614149); void init() { // Eta ranges Cut eta_full = (Cuts::abseta < 5.0); Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV); // All final state particles FinalState fs(eta_full); IdentifiedFinalState all_photons(fs); all_photons.acceptIdPair(PID::PHOTON); // Get photons to dress leptons IdentifiedFinalState ph_id(fs); ph_id.acceptIdPair(PID::PHOTON); // Projection to find the electrons IdentifiedFinalState el_id(fs); el_id.acceptIdPair(PID::ELECTRON); PromptFinalState photons(ph_id); photons.acceptTauDecays(true); declare(photons, "photons"); PromptFinalState electrons(el_id); electrons.acceptTauDecays(true); DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts); declare(dressedelectrons, "elecs"); DressedLeptons ewdressedelectrons(all_photons, electrons, 0.1, eta_full); // Projection to find the muons IdentifiedFinalState mu_id(fs); mu_id.acceptIdPair(PID::MUON); PromptFinalState muons(mu_id); muons.acceptTauDecays(true); DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts); declare(dressedmuons, "muons"); DressedLeptons ewdressedmuons(all_photons, muons, 0.1, eta_full); // Projection to find MET declare(MissingMomentum(fs), "MET"); // remove prompt neutrinos from jet clustering IdentifiedFinalState nu_id(fs); nu_id.acceptNeutrinos(); PromptFinalState neutrinos(nu_id); neutrinos.acceptTauDecays(true); // Jet clustering. VetoedFinalState vfs(fs); vfs.addVetoOnThisFinalState(ewdressedelectrons); vfs.addVetoOnThisFinalState(ewdressedmuons); vfs.addVetoOnThisFinalState(neutrinos); - FastJets jets(vfs, FastJets::ANTIKT, 0.4); - jets.useInvisibles(true); + FastJets jets(vfs, FastJets::ANTIKT, 0.4, JetAlg::Muons::ALL, JetAlg::Invisibles::DECAY); declare(jets, "jets"); // Addition of the large-R jets VetoedFinalState vfs1(fs); vfs1.addVetoOnThisFinalState(neutrinos); FastJets fjets(vfs1, FastJets::ANTIKT, 1.); fjets.useInvisibles(JetAlg::Invisibles::NONE); fjets.useMuons(JetAlg::Muons::NONE); declare(fjets, "fjets"); bookHists("top_pt_res", 15); bookHists("top_absrap_res", 17); bookHists("ttbar_pt_res", 19); bookHists("ttbar_absrap_res", 21); bookHists("ttbar_m_res", 23); bookHists("top_pt_boost", 25); bookHists("top_absrap_boost", 27); } void analyze(const Event& event) { // Get the selected objects, using the projections. vector electrons = apply(event, "elecs").dressedLeptons(); vector muons = apply(event, "muons").dressedLeptons(); const Jets& jets = apply(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); const PseudoJets& all_fjets = apply(event, "fjets").pseudoJetsByPt(); // get MET const Vector3 met = apply(event, "MET").vectorMPT(); Jets bjets, lightjets; for (const Jet& jet : jets) { bool b_tagged = jet.bTags(Cuts::pT > 5*GeV).size(); if ( b_tagged && bjets.size() < 2) bjets +=jet; else lightjets += jet; } // Implementing large-R jets definition // trim the jets PseudoJets trimmed_fatJets; float Rfilt = 0.2; float pt_fraction_min = 0.05; fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, Rfilt), fastjet::SelectorPtFractionMin(pt_fraction_min)); for (PseudoJet pjet : all_fjets) trimmed_fatJets += trimmer(pjet); trimmed_fatJets = fastjet::sorted_by_pt(trimmed_fatJets); PseudoJets trimmed_jets; for (unsigned int i = 0; i < trimmed_fatJets.size(); ++i) { FourMomentum tj_mom = momentum(trimmed_fatJets[i]); if (tj_mom.pt() <= 300*GeV) continue; if (tj_mom.abseta() >= 2.0) continue; trimmed_jets.push_back(trimmed_fatJets[i]); } bool single_electron = (electrons.size() == 1) && (muons.empty()); bool single_muon = (muons.size() == 1) && (electrons.empty()); DressedLepton *lepton = NULL; if (single_electron) lepton = &electrons[0]; else if (single_muon) lepton = &muons[0]; if (!single_electron && !single_muon) vetoEvent; bool pass_resolved = true; bool num_b_tagged_jets = (bjets.size() == 2); if (!num_b_tagged_jets) pass_resolved = false; if (jets.size() < 4) pass_resolved = false; bool pass_boosted = true; int fatJetIndex = -1; bool passTopTag = false; bool passDphi = false; bool passAddJet = false; bool goodLepJet = false; bool lepbtag = false; bool hadbtag=false; vector lepJetIndex; vector jet_farFromHadTopJetCandidate; if (met.mod() < 20*GeV) pass_boosted = false; if (pass_boosted) { double transmass = _mT(lepton->momentum(), met); if (transmass + met.mod() < 60*GeV) pass_boosted = false; } if (pass_boosted) { if (trimmed_jets.size() >= 1) { for (unsigned int j = 0; j 100*GeV && momentum(trimmed_jets.at(j)).pt() > 300*GeV && momentum(trimmed_jets.at(j)).pt() < 1500*GeV && fabs(momentum(trimmed_jets.at(j)).eta()) < 2.) { passTopTag = true; fatJetIndex = j; break; } } } } if(!passTopTag && fatJetIndex == -1) pass_boosted = false; if (pass_boosted) { double dPhi_fatjet = deltaPhi(lepton->phi(), momentum(trimmed_jets.at(fatJetIndex)).phi()); double dPhi_fatjet_lep_cut = 1.0; //2.3 if (dPhi_fatjet > dPhi_fatjet_lep_cut ) { passDphi = true; } } if (!passDphi) pass_boosted = false; if (bjets.empty()) pass_boosted = false; if (pass_boosted) { for (unsigned int sj = 0; sj < jets.size(); ++sj) { double dR = deltaR(jets.at(sj).momentum(), momentum(trimmed_jets.at(fatJetIndex))); if(dR > 1.5) { passAddJet = true; jet_farFromHadTopJetCandidate.push_back(sj); } } } if (!passAddJet) pass_boosted = false; if (pass_boosted) { for (int ltj : jet_farFromHadTopJetCandidate) { double dR_jet_lep = deltaR(jets.at(ltj).momentum(), lepton->momentum()); double dR_jet_lep_cut = 2.0;//1.5 if (dR_jet_lep < dR_jet_lep_cut) { lepJetIndex.push_back(ltj); goodLepJet = true; } } } if(!goodLepJet) pass_boosted = false; if (pass_boosted) { for (int lepj : lepJetIndex) { lepbtag = jets.at(lepj).bTags(Cuts::pT > 5*GeV).size(); if (lepbtag) break; } } double dR_fatBjet_cut = 1.0; if (pass_boosted) { for (const Jet& bjet : bjets) { hadbtag |= deltaR(momentum(trimmed_jets.at(fatJetIndex)), bjet) < dR_fatBjet_cut; } } if (!(lepbtag || hadbtag)) pass_boosted = false; FourMomentum pbjet1; //Momentum of bjet1 FourMomentum pbjet2; //Momentum of bjet int Wj1index = -1, Wj2index = -1; if (pass_resolved) { if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) { pbjet1 = bjets[0].momentum(); pbjet2 = bjets[1].momentum(); } else { pbjet1 = bjets[1].momentum(); pbjet2 = bjets[0].momentum(); } double bestWmass = 1000.0*TeV; double mWPDG = 80.399*GeV; for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) { for (unsigned int j = i + 1; j < lightjets.size(); ++j) { double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass(); if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) { bestWmass = wmass; Wj1index = i; Wj2index = j; } } } FourMomentum pjet1 = lightjets[Wj1index].momentum(); FourMomentum pjet2 = lightjets[Wj2index].momentum(); // compute hadronic W boson FourMomentum pWhadron = pjet1 + pjet2; double pz = computeneutrinoz(lepton->momentum(), met); FourMomentum ppseudoneutrino( sqrt(sqr(met.x()) + sqr(met.y()) + sqr(pz)), met.x(), met.y(), pz); //compute leptonic, hadronic, combined pseudo-top FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1; FourMomentum ppseudotophadron = pbjet2 + pWhadron; FourMomentum pttbar = ppseudotoplepton + ppseudotophadron; fillHists("top_pt_res", ppseudotophadron.pt()/GeV); fillHists("top_absrap_res", ppseudotophadron.absrap()); fillHists("ttbar_pt_res", pttbar.pt()/GeV); fillHists("ttbar_absrap_res", pttbar.absrap()); fillHists("ttbar_m_res", pttbar.mass()/GeV); } - + if (pass_boosted) {// Boosted selection double hadtop_pt= momentum(trimmed_jets.at(fatJetIndex)).pt() / GeV; double hadtop_absrap= momentum(trimmed_jets.at(fatJetIndex)).absrap(); fillHists("top_pt_boost", hadtop_pt); fillHists("top_absrap_boost", hadtop_absrap); } } void finalize() { // Normalize to cross-section const double sf = (crossSection() / sumOfWeights()); for (HistoMap::value_type& hist : _h) { scale(hist.second, sf); if (hist.first.find("_norm") != string::npos) normalize(hist.second); } } void bookHists(std::string name, unsigned int index) { book(_h[name], index, 1 ,1); book(_h[name + "_norm"], index + 1, 1 ,1); } void fillHists(std::string name, double value) { _h[name]->fill(value); _h[name + "_norm"]->fill(value); } double _mT(const FourMomentum &l, const Vector3 &met) const { return sqrt(2.0 * l.pT() * met.mod() * (1 - cos(deltaPhi(l, met))) ); } double tau32(const fastjet::PseudoJet &jet, double jet_rad) const { double alpha = 1.0; Nsubjettiness::NormalizedCutoffMeasure normalized_measure(alpha, jet_rad, 1000000); // WTA definition // Nsubjettiness::OnePass_WTA_KT_Axes wta_kt_axes; - // as in JetSubStructure recommendations + // as in JetSubStructure recommendations Nsubjettiness::KT_Axes kt_axes; /// NsubjettinessRatio uses the results from Nsubjettiness to calculate the ratio /// tau_N/tau_M, where N and M are specified by the user. The ratio of different tau values /// is often used in analyses, so this class is helpful to streamline code. Nsubjettiness::NsubjettinessRatio tau32_kt(3, 2, kt_axes, normalized_measure); double tau32 = tau32_kt.result(jet); return tau32; } double computeneutrinoz(const FourMomentum& lepton, const Vector3 &met) const { //computing z component of neutrino momentum given lepton and met double pzneutrino; double m_W = 80.399; // in GeV, given in the paper double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.x() + lepton.py() * met.y()); double a = sqr ( lepton.E() )- sqr ( lepton.pz() ); double b = -2*k*lepton.pz(); double c = sqr( lepton.E() ) * sqr( met.mod() ) - sqr( k ); double discriminant = sqr(b) - 4 * a * c; double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns if (discriminant < 0) pzneutrino = - b / (2 * a); //if the discriminant is negative else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value double absquad[2]; for (int n=0; n<2; ++n) absquad[n] = fabs(quad[n]); if (absquad[0] < absquad[1]) pzneutrino = quad[0]; else pzneutrino = quad[1]; } return pzneutrino; } private: /// @name Objects that are used by the event selection decisions typedef map HistoMap; HistoMap _h; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2017_I1614149); }