diff --git a/analyses/pluginCDF/CDF_1996_S3349578.cc b/analyses/pluginCDF/CDF_1996_S3349578.cc --- a/analyses/pluginCDF/CDF_1996_S3349578.cc +++ b/analyses/pluginCDF/CDF_1996_S3349578.cc @@ -1,448 +1,448 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/SmearedJets.hh" namespace Rivet { /// @brief CDF properties of high-mass multi-jet events class CDF_1996_S3349578 : public Analysis { public: /// Constructor DEFAULT_RIVET_ANALYSIS_CTOR(CDF_1996_S3349578); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { /// Initialise and register projections here - const FinalState fs(-4.2, 4.2); + const FinalState fs(Cuts::abseta < 4.2); FastJets fj(fs, FastJets::CDFJETCLU, 0.7); declare(fj, "Jets"); // Smear energy and mass with the 10% uncertainty quoted in the paper SmearedJets sj_E(fj, [](const Jet& jet){ return P4_SMEAR_MASS_GAUSS(P4_SMEAR_E_GAUSS(jet, 0.1*jet.E()), 0.1*jet.mass()); }); declare(sj_E, "SmearedJets"); /// Book histograms here, e.g.: _h_3_mNJ = bookHisto1D(1, 1, 1); _h_3_X3 = bookHisto1D(2, 1, 1); _h_3_X4 = bookHisto1D(3, 1, 1); _h_3_costheta3 = bookHisto1D(8, 1, 1); _h_3_psi3 = bookHisto1D(9, 1, 1); _h_3_f3 = bookHisto1D(14, 1, 1); _h_3_f4 = bookHisto1D(14, 1, 2); _h_3_f5 = bookHisto1D(14, 1, 3); _h_4_mNJ = bookHisto1D(1, 1, 2); _h_4_X3 = bookHisto1D(4, 1, 1); _h_4_X4 = bookHisto1D(5, 1, 1); _h_4_costheta3 = bookHisto1D(10, 1, 1); _h_4_psi3 = bookHisto1D(11, 1, 1); _h_4_f3 = bookHisto1D(15, 1, 1); _h_4_f4 = bookHisto1D(15, 1, 2); _h_4_f5 = bookHisto1D(15, 1, 3); _h_4_XA = bookHisto1D(17, 1, 1); _h_4_psiAB = bookHisto1D(19, 1, 1); _h_4_fA = bookHisto1D(21, 1, 1); _h_4_fB = bookHisto1D(21, 1, 2); _h_5_mNJ = bookHisto1D(1, 1, 3); _h_5_X3 = bookHisto1D(6, 1, 1); _h_5_X4 = bookHisto1D(7, 1, 1); _h_5_costheta3 = bookHisto1D(12, 1, 1); _h_5_psi3 = bookHisto1D(13, 1, 1); _h_5_f3 = bookHisto1D(16, 1, 1); _h_5_f4 = bookHisto1D(16, 1, 2); _h_5_f5 = bookHisto1D(16, 1, 3); _h_5_XA = bookHisto1D(18, 1, 1); _h_5_XC = bookHisto1D(18, 1, 2); _h_5_psiAB = bookHisto1D(20, 1, 1); _h_5_psiCD = bookHisto1D(20, 1, 2); _h_5_fA = bookHisto1D(22, 1, 1); _h_5_fB = bookHisto1D(23, 1, 1); _h_5_fC = bookHisto1D(24, 1, 1); _h_5_fD = bookHisto1D(25, 1, 1); } void analyze(const Event& event) { Jets jets; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); foreach (const Jet& jet, apply(event, "SmearedJets").jets(Cuts::Et > 20.0*GeV, cmpMomByEt)) { bool separated = true; foreach (const Jet& ref, jets) { if (deltaR(jet, ref) < 0.9) { separated = false; break; } } if (!separated) continue; jets.push_back(jet); jetsystem += jet.momentum(); if (jets.size() >= 5) break; } const double weight = event.weight(); if (jets.size() > 4) { _fiveJetAnalysis(jets, weight); jets.resize(4); } if (jets.size() > 3) { _fourJetAnalysis(jets, weight); jets.resize(3); } if (jets.size() > 2) { _threeJetAnalysis(jets, weight); } } void _threeJetAnalysis(const Jets& jets, const double& weight) { MSG_DEBUG("3 jet analysis"); double sumEt = 0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); foreach (const Jet& jet, jets) { sumEt += jet.Et(); jetsystem += jet.momentum(); } if (sumEt < 420.0*GeV) return; const double m3J = _safeMass(jetsystem); if (m3J < 600*GeV) return; const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec()); vector jets3; foreach (Jet jet, jets) { jets3.push_back(cms_boost.transform(jet.momentum())); } std::sort(jets3.begin(), jets3.end(), FourMomentum::byEDescending()); FourMomentum p3(jets3[0]), p4(jets3[1]), p5(jets3[2]); FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m3J, jetsystem.rapidity())); double costheta3 = pAV.p3().unit().dot(p3.p3().unit()); if (fabs(costheta3) > 0.6) return; double X3 = 2.0*p3.E()/m3J; if (X3 > 0.9) return; const double X4 = 2.0*p4.E()/m3J; const double psi3 = _psi(p3, pAV, p4, p5); const double f3 = _safeMass(p3)/m3J; const double f4 = _safeMass(p4)/m3J; const double f5 = _safeMass(p5)/m3J; _h_3_mNJ->fill(m3J, weight); _h_3_X3->fill(X3, weight); _h_3_X4->fill(X4, weight); _h_3_costheta3->fill(costheta3, weight); _h_3_psi3->fill(psi3, weight); _h_3_f3->fill(f3, weight); _h_3_f4->fill(f4, weight); _h_3_f5->fill(f5, weight); } void _fourJetAnalysis(const Jets& jets, const double& weight) { MSG_DEBUG("4 jet analysis"); double sumEt=0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); foreach (const Jet& jet, jets) { sumEt+=jet.Et(); jetsystem+=jet.momentum(); } if (sumEt < 420.0*GeV) return; const double m4J = _safeMass(jetsystem); if (m4J < 650*GeV) return; const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec()); vector jets4; foreach (Jet jet, jets) { jets4.push_back(cms_boost.transform(jet.momentum())); } std::sort(jets4.begin(), jets4.end(), FourMomentum::byEDescending()); FourMomentum pA, pB; vector jets3(_reduce(jets4, pA, pB)); std::sort(jets3.begin(), jets3.end(), FourMomentum::byEDescending()); FourMomentum p3(jets3[0]); FourMomentum p4(jets3[1]); FourMomentum p5(jets3[2]); FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m4J, jetsystem.rapidity())); double costheta3=pAV.p3().unit().dot(p3.p3().unit()); if (fabs(costheta3)>0.8) { return; } const double X3 = 2.0*p3.E()/m4J; if (X3>0.9) { return; } // fill histograms const double X4 = 2.0*p4.E()/m4J; const double psi3 = _psi(p3, pAV, p4, p5); const double f3 = _safeMass(p3)/m4J; const double f4 = _safeMass(p4)/m4J; const double f5 = _safeMass(p5)/m4J; const double fA = _safeMass(pA)/m4J; const double fB = _safeMass(pB)/m4J; const double XA = pA.E()/(pA.E()+pB.E()); const double psiAB = _psi(pA, pB, pA+pB, pAV); _h_4_mNJ->fill(m4J, weight); _h_4_X3->fill(X3, weight); _h_4_X4->fill(X4, weight); _h_4_costheta3->fill(costheta3, weight); _h_4_psi3->fill(psi3, weight); _h_4_f3->fill(f3, weight); _h_4_f4->fill(f4, weight); _h_4_f5->fill(f5, weight); _h_4_XA->fill(XA, weight); _h_4_psiAB->fill(psiAB, weight); _h_4_fA->fill(fA, weight); _h_4_fB->fill(fB, weight); } void _fiveJetAnalysis(const Jets& jets, const double& weight) { MSG_DEBUG("5 jet analysis"); double sumEt=0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); foreach (const Jet& jet, jets) { sumEt+=jet.Et(); jetsystem+=jet.momentum(); } if (sumEt < 420.0*GeV) return; const double m5J = _safeMass(jetsystem); if (m5J < 750*GeV) return; const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec()); vector jets5; foreach (Jet jet, jets) { jets5.push_back(cms_boost.transform(jet.momentum())); } std::sort(jets5.begin(), jets5.end(), FourMomentum::byEDescending()); FourMomentum pC, pD; vector jets4(_reduce(jets5, pC, pD)); std::sort(jets4.begin(), jets4.end(), FourMomentum::byEDescending()); FourMomentum pA, pB; vector jets3(_reduce(jets4, pA, pB)); std::sort(jets3.begin(), jets3.end(), FourMomentum::byEDescending()); FourMomentum p3(jets3[0]); FourMomentum p4(jets3[1]); FourMomentum p5(jets3[2]); // fill histograms FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m5J, jetsystem.rapidity())); const double costheta3 = pAV.p3().unit().dot(p3.p3().unit()); const double X3 = 2.0*p3.E()/m5J; const double X4 = 2.0*p4.E()/m5J; const double psi3 = _psi(p3, pAV, p4, p5); const double f3 = _safeMass(p3)/m5J; const double f4 = _safeMass(p4)/m5J; const double f5 = _safeMass(p5)/m5J; const double fA = _safeMass(pA)/m5J; const double fB = _safeMass(pB)/m5J; const double XA = pA.E()/(pA.E()+pB.E()); const double psiAB = _psi(pA, pB, pA+pB, pAV); const double fC = _safeMass(pC)/m5J; const double fD = _safeMass(pD)/m5J; const double XC = pC.E()/(pC.E()+pD.E()); const double psiCD = _psi(pC, pD, pC+pD, pAV); _h_5_mNJ->fill(m5J, weight); _h_5_X3->fill(X3, weight); _h_5_X4->fill(X4, weight); _h_5_costheta3->fill(costheta3, weight); _h_5_psi3->fill(psi3, weight); _h_5_f3->fill(f3, weight); _h_5_f4->fill(f4, weight); _h_5_f5->fill(f5, weight); _h_5_XA->fill(XA, weight); _h_5_psiAB->fill(psiAB, weight); _h_5_fA->fill(fA, weight); _h_5_fB->fill(fB, weight); _h_5_XC->fill(XC, weight); _h_5_psiCD->fill(psiCD, weight); _h_5_fC->fill(fC, weight); _h_5_fD->fill(fD, weight); } /// Normalise histograms etc., after the run void finalize() { /// Normalise, scale and otherwise manipulate histograms here normalize(_h_3_mNJ, 1.0); normalize(_h_3_X3, 1.0); normalize(_h_3_X4, 1.0); normalize(_h_3_costheta3, 1.0); normalize(_h_3_psi3, 1.0); normalize(_h_3_f3, 1.0); normalize(_h_3_f4, 1.0); normalize(_h_3_f5, 1.0); normalize(_h_4_mNJ, 1.0); normalize(_h_4_X3, 1.0); normalize(_h_4_X4, 1.0); normalize(_h_4_costheta3, 1.0); normalize(_h_4_psi3, 1.0); normalize(_h_4_f3, 1.0); normalize(_h_4_f4, 1.0); normalize(_h_4_f5, 1.0); normalize(_h_4_XA, 1.0); normalize(_h_4_psiAB, 1.0); normalize(_h_4_fA, 1.0); normalize(_h_4_fB, 1.0); normalize(_h_5_mNJ, 1.0); normalize(_h_5_X3, 1.0); normalize(_h_5_X4, 1.0); normalize(_h_5_costheta3, 1.0); normalize(_h_5_psi3, 1.0); normalize(_h_5_f3, 1.0); normalize(_h_5_f4, 1.0); normalize(_h_5_f5, 1.0); normalize(_h_5_XA, 1.0); normalize(_h_5_XC, 1.0); normalize(_h_5_psiAB, 1.0); normalize(_h_5_psiCD, 1.0); normalize(_h_5_fA, 1.0); normalize(_h_5_fB, 1.0); normalize(_h_5_fC, 1.0); normalize(_h_5_fD, 1.0); } //@} private: vector _reduce(const vector& jets, FourMomentum& combined1, FourMomentum& combined2) { double minMass2 = 1e9; size_t idx1(jets.size()), idx2(jets.size()); for (size_t i=0; i newjets; for (size_t i=0; i1e-3) { FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y)); const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(boostvec.betaVec()).inverse(); beam1 = cms_boost.transform(beam1); beam2 = cms_boost.transform(beam2); } return (beam1.E() > beam2.E()) ? beam1-beam2 : beam2-beam1; } double _psi(const FourMomentum& p1, const FourMomentum& p2, const FourMomentum& p3, const FourMomentum& p4) { Vector3 p1xp2 = p1.p3().cross(p2.p3()); Vector3 p3xp4 = p3.p3().cross(p4.p3()); return mapAngle0ToPi(acos(p1xp2.unit().dot(p3xp4.unit()))); } double _safeMass(const FourMomentum& p) { double mass2=p.mass2(); if (mass2>0.0) return sqrt(mass2); else if (mass2<-1.0e-5) { MSG_WARNING("m2 = " << m2 << ". Assuming m2=0."); return 0.0; } else return 0.0; } private: /// @name Histograms //@{ Histo1DPtr _h_3_mNJ; Histo1DPtr _h_3_X3; Histo1DPtr _h_3_X4; Histo1DPtr _h_3_costheta3; Histo1DPtr _h_3_psi3; Histo1DPtr _h_3_f3; Histo1DPtr _h_3_f4; Histo1DPtr _h_3_f5; Histo1DPtr _h_4_mNJ; Histo1DPtr _h_4_X3; Histo1DPtr _h_4_X4; Histo1DPtr _h_4_costheta3; Histo1DPtr _h_4_psi3; Histo1DPtr _h_4_f3; Histo1DPtr _h_4_f4; Histo1DPtr _h_4_f5; Histo1DPtr _h_4_XA; Histo1DPtr _h_4_psiAB; Histo1DPtr _h_4_fA; Histo1DPtr _h_4_fB; Histo1DPtr _h_5_mNJ; Histo1DPtr _h_5_X3; Histo1DPtr _h_5_X4; Histo1DPtr _h_5_costheta3; Histo1DPtr _h_5_psi3; Histo1DPtr _h_5_f3; Histo1DPtr _h_5_f4; Histo1DPtr _h_5_f5; Histo1DPtr _h_5_XA; Histo1DPtr _h_5_XC; Histo1DPtr _h_5_psiAB; Histo1DPtr _h_5_psiCD; Histo1DPtr _h_5_fA; Histo1DPtr _h_5_fB; Histo1DPtr _h_5_fC; Histo1DPtr _h_5_fD; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CDF_1996_S3349578); } diff --git a/analyses/pluginCDF/CDF_1997_S3541940.cc b/analyses/pluginCDF/CDF_1997_S3541940.cc --- a/analyses/pluginCDF/CDF_1997_S3541940.cc +++ b/analyses/pluginCDF/CDF_1997_S3541940.cc @@ -1,250 +1,250 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/SmearedJets.hh" namespace Rivet { /// @brief CDF properties of 6-jet events with large 6-jet mass class CDF_1997_S3541940 : public Analysis { public: DEFAULT_RIVET_ANALYSIS_CTOR(CDF_1997_S3541940); void init() { // Find true jets - const FinalState fs(-4.2, 4.2); + const FinalState fs(Cuts::abseta < 4.2); FastJets fj(fs, FastJets::CDFJETCLU, 0.7); // Smear jet energy and mass with the 10% uncertainty quoted in the paper SmearedJets sj_E(fj, [](const Jet& jet){ return P4_SMEAR_MASS_GAUSS(P4_SMEAR_E_GAUSS(jet, 0.1*jet.E()), 0.1*jet.mass()); }); declare(sj_E, "Jets"); _h_m6J = bookHisto1D(1, 1, 1); _h_X3ppp = bookHisto1D(2, 1, 1); _h_X4ppp = bookHisto1D(3, 1, 1); _h_costheta3ppp = bookHisto1D(4, 1, 1); _h_psi3ppp = bookHisto1D(5, 1, 1); _h_f3ppp = bookHisto1D(6, 1, 1); _h_f4ppp = bookHisto1D(6, 1, 2); _h_f5ppp = bookHisto1D(6, 1, 3); _h_XApp = bookHisto1D(7, 1, 1); _h_XCp = bookHisto1D(8, 1, 1); _h_XE = bookHisto1D(9, 1, 1); _h_psiAppBpp = bookHisto1D(10, 1, 1); _h_psiCpDp = bookHisto1D(11, 1, 1); _h_psiEF = bookHisto1D(12, 1, 1); _h_fApp = bookHisto1D(13, 1, 1); _h_fBpp = bookHisto1D(14, 1, 1); _h_fCp = bookHisto1D(15, 1, 1); _h_fDp = bookHisto1D(16, 1, 1); _h_fE = bookHisto1D(17, 1, 1); _h_fF = bookHisto1D(18, 1, 1); } void analyze(const Event& event) { const Jets alljets = apply(event, "Jets").jets(Cuts::Et > 20*GeV && Cuts::abseta < 3, cmpMomByEt); Jets jets; double sumEt = 0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); for (const Jet& jet : alljets) { double Et = jet.Et(); bool separated = true; for (const Jet& ref : jets) { if (deltaR(jet, ref) < 0.9) { separated = false; break; } } if (!separated) continue; jets.push_back(jet); sumEt += Et; jetsystem += jet.momentum(); if (jets.size() >= 6) break; } if (jets.size() < 6) vetoEvent; if (sumEt < 320.0*GeV) vetoEvent; double m6J = _safeMass(jetsystem); if (m6J < 520.0*GeV) vetoEvent; if (getLog().isActive(Log::DEBUG)) { stringstream ss; ss << "Jets:\n"; for (const Jet& j : jets) ss << j << "\n"; MSG_DEBUG(ss.str()); } const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(jetsystem.betaVec()); vector jets6; for (Jet jet : jets) { jets6.push_back(cms_boost.transform(jet.momentum())); } std::sort(jets6.begin(), jets6.end(), FourMomentum::byEDescending()); FourMomentum pE, pF; vector jets5(_reduce(jets6, pE, pF)); std::sort(jets5.begin(), jets5.end(), FourMomentum::byEDescending()); FourMomentum pCp, pDp; vector jets4(_reduce(jets5, pCp, pDp)); std::sort(jets4.begin(), jets4.end(), FourMomentum::byEDescending()); FourMomentum pApp, pBpp; vector jets3(_reduce(jets4, pApp, pBpp)); std::sort(jets3.begin(), jets3.end(), FourMomentum::byEDescending()); FourMomentum p3ppp(jets3[0]); FourMomentum p4ppp(jets3[1]); FourMomentum p5ppp(jets3[2]); double X3ppp = 2.0*p3ppp.E()/m6J; if (X3ppp > 0.9) vetoEvent; FourMomentum pAV = cms_boost.transform(_avg_beam_in_lab(m6J, jetsystem.rapidity())); double costheta3ppp = pAV.p3().unit().dot(p3ppp.p3().unit()); if (fabs(costheta3ppp) > 0.9) vetoEvent; const double weight = event.weight(); // 3-jet-system variables _h_m6J->fill(m6J, weight); _h_X3ppp->fill(X3ppp, weight); _h_X4ppp->fill(2.0*p4ppp.E()/m6J, weight); _h_costheta3ppp->fill(costheta3ppp, weight); double psi3ppp = _psi(p3ppp, pAV, p4ppp, p5ppp); _h_psi3ppp->fill(psi3ppp, weight); _h_f3ppp->fill(_safeMass(p3ppp)/m6J, weight); _h_f4ppp->fill(_safeMass(p4ppp)/m6J, weight); _h_f5ppp->fill(_safeMass(p5ppp)/m6J, weight); // 4 -> 3 jet variables _h_fApp->fill(_safeMass(pApp)/m6J, weight); _h_fBpp->fill(_safeMass(pApp)/m6J, weight); _h_XApp->fill(pApp.E()/(pApp.E()+pBpp.E()), weight); double psiAppBpp = _psi(pApp, pBpp, pApp+pBpp, pAV); _h_psiAppBpp->fill(psiAppBpp, weight); // 5 -> 4 jet variables _h_fCp->fill(_safeMass(pCp)/m6J, weight); _h_fDp->fill(_safeMass(pDp)/m6J, weight); _h_XCp->fill(pCp.E()/(pCp.E()+pDp.E()), weight); double psiCpDp = _psi(pCp, pDp, pCp+pDp, pAV); _h_psiCpDp->fill(psiCpDp, weight); // 6 -> 5 jet variables _h_fE->fill(_safeMass(pE)/m6J, weight); _h_fF->fill(_safeMass(pF)/m6J, weight); _h_XE->fill(pE.E()/(pE.E()+pF.E()), weight); double psiEF = _psi(pE, pF, pE+pF, pAV); _h_psiEF->fill(psiEF, weight); } void finalize() { normalize(_h_m6J); normalize(_h_X3ppp); normalize(_h_X4ppp); normalize(_h_costheta3ppp); normalize(_h_psi3ppp); normalize(_h_f3ppp); normalize(_h_f4ppp); normalize(_h_f5ppp); normalize(_h_XApp); normalize(_h_XCp); normalize(_h_XE); normalize(_h_psiAppBpp); normalize(_h_psiCpDp); normalize(_h_psiEF); normalize(_h_fApp); normalize(_h_fBpp); normalize(_h_fCp); normalize(_h_fDp); normalize(_h_fE); normalize(_h_fF); } private: vector _reduce(const vector& jets, FourMomentum& combined1, FourMomentum& combined2) { double minMass2 = 1e9; size_t idx1(jets.size()), idx2(jets.size()); for (size_t i = 0; i < jets.size(); ++i) { for (size_t j = i+1; j < jets.size(); ++j) { double mass2 = FourMomentum(jets[i] + jets[j]).mass2(); if (mass2 < minMass2) { idx1 = i; idx2 = j; } } } vector newjets; for (size_t i = 0; i < jets.size(); ++i) { if (i != idx1 && i != idx2) newjets.push_back(jets[i]); } newjets.push_back(jets[idx1] + jets[idx2]); combined1 = jets[idx1]; combined2 = jets[idx2]; return newjets; } FourMomentum _avg_beam_in_lab(const double& m, const double& y) { const double mt = m/2.0; FourMomentum beam1(mt, 0, 0, mt); FourMomentum beam2(mt, 0, 0, -mt); if (fabs(y) > 1e-3) { FourMomentum boostvec(cosh(y), 0.0, 0.0, sinh(y)); const LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(boostvec.betaVec()).inverse(); beam1 = cms_boost.transform(beam1); beam2 = cms_boost.transform(beam2); } return (beam1.E() > beam2.E()) ? beam1 - beam2 : beam2 - beam1; } double _psi(const FourMomentum& p1, const FourMomentum& p2, const FourMomentum& p3, const FourMomentum& p4) { Vector3 p1xp2 = p1.p3().cross(p2.p3()); Vector3 p3xp4 = p3.p3().cross(p4.p3()); return mapAngle0ToPi(acos(p1xp2.unit().dot(p3xp4.unit()))); } double _safeMass(const FourMomentum& p) { double mass2 = p.mass2(); if (mass2 > 0.0) return sqrt(mass2); if (mass2 < -1e-5) MSG_WARNING("m2 = " << m2 << ". Assuming m2=0."); return 0.0; } private: Histo1DPtr _h_m6J; Histo1DPtr _h_X3ppp, _h_X4ppp; Histo1DPtr _h_costheta3ppp; Histo1DPtr _h_psi3ppp; Histo1DPtr _h_f3ppp, _h_f4ppp, _h_f5ppp; Histo1DPtr _h_XApp, _h_XCp, _h_XE; Histo1DPtr _h_psiAppBpp, _h_psiCpDp, _h_psiEF; Histo1DPtr _h_fApp, _h_fBpp, _h_fCp, _h_fDp, _h_fE, _h_fF; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CDF_1997_S3541940); } diff --git a/analyses/pluginCDF/CDF_2001_S4751469.cc b/analyses/pluginCDF/CDF_2001_S4751469.cc --- a/analyses/pluginCDF/CDF_2001_S4751469.cc +++ b/analyses/pluginCDF/CDF_2001_S4751469.cc @@ -1,260 +1,260 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/TriggerCDFRun0Run1.hh" #include "Rivet/Projections/ConstLossyFinalState.hh" //#include "Rivet/Projections/SmearedParticles.hh" namespace Rivet { /// @brief Field-Stuart CDF Run I track-jet underlying event analysis /// /// @author Andy Buckley /// /// The "original" underlying event analysis, using a non-standard track-jet algorithm. /// /// @par Run conditions /// /// @arg \f$ \sqrt{s} = \f$ 1800 GeV /// @arg Run with generic QCD events. /// @arg Several \f$ p_\perp^\text{min} \f$ cutoffs are probably required to fill the profile histograms: /// @arg \f$ p_\perp^\text{min} = \f$ 0 (min bias), 10, 20 GeV class CDF_2001_S4751469 : public Analysis { public: /// Constructor: cuts on final state are \f$ -1 < \eta < 1 \f$ /// and \f$ p_T > 0.5 \f$ GeV. CDF_2001_S4751469() : Analysis("CDF_2001_S4751469"), _totalNumTrans2(0), _totalNumTrans5(0), _totalNumTrans30(0), _sumWeightsPtLead2(0),_sumWeightsPtLead5(0), _sumWeightsPtLead30(0) { } /// @name Analysis methods //@{ // Book histograms void init() { declare(TriggerCDFRun0Run1(), "Trigger"); // Randomly discard 8% of charged particles as a kind of hacky detector correction. - const ChargedFinalState cfs(-1.0, 1.0, 0.5*GeV); + const ChargedFinalState cfs(Cuts::abseta < 1.0 && Cuts::pT > 0.5*GeV); /// @todo Replace ConstLossyFinalState with SmearedParticles const ConstLossyFinalState lossyfs(cfs, 0.08); //const SmearedParticles lossyfs(cfs, [](const Particle&){ return 0.92; }); declare(lossyfs, "FS"); declare(FastJets(lossyfs, FastJets::TRACKJET, 0.7), "TrackJet"); _numvsDeltaPhi2 = bookProfile1D(1, 1, 1); _numvsDeltaPhi5 = bookProfile1D(1, 1, 2); _numvsDeltaPhi30 = bookProfile1D(1, 1, 3); _pTvsDeltaPhi2 = bookProfile1D(2, 1, 1); _pTvsDeltaPhi5 = bookProfile1D(2, 1, 2); _pTvsDeltaPhi30 = bookProfile1D(2, 1, 3); _numTowardMB = bookProfile1D(3, 1, 1); _numTransMB = bookProfile1D(3, 1, 2); _numAwayMB = bookProfile1D(3, 1, 3); _numTowardJ20 = bookProfile1D(4, 1, 1); _numTransJ20 = bookProfile1D(4, 1, 2); _numAwayJ20 = bookProfile1D(4, 1, 3); _ptsumTowardMB = bookProfile1D(5, 1, 1); _ptsumTransMB = bookProfile1D(5, 1, 2); _ptsumAwayMB = bookProfile1D(5, 1, 3); _ptsumTowardJ20 = bookProfile1D(6, 1, 1); _ptsumTransJ20 = bookProfile1D(6, 1, 2); _ptsumAwayJ20 = bookProfile1D(6, 1, 3); _ptTrans2 = bookHisto1D(7, 1, 1); _ptTrans5 = bookHisto1D(7, 1, 2); _ptTrans30 = bookHisto1D(7, 1, 3); } /// Do the analysis void analyze(const Event& event) { // Trigger const bool trigger = apply(event, "Trigger").minBiasDecision(); if (!trigger) vetoEvent; // Get jets, sorted by pT const Jets jets = apply(event, "TrackJet").jetsByPt(); if (jets.empty()) vetoEvent; const Jet jet1 = jets.front(); const double ptLead = jet1.pT(); // Cut on highest pT jet: combined 0.5 GeV < pT(lead) < 50 GeV if (ptLead/GeV < 0.5) vetoEvent; if (ptLead/GeV > 50.0) vetoEvent; // Count sum of all event weights in three pT_lead regions const double weight = event.weight(); if (ptLead/GeV > 2.0) _sumWeightsPtLead2 += weight; if (ptLead/GeV > 5.0) _sumWeightsPtLead5 += weight; if (ptLead/GeV > 30.0) _sumWeightsPtLead30 += weight; // Run over tracks double ptSumToward(0.0), ptSumAway(0.0), ptSumTrans(0.0); size_t numToward(0), numTrans(0), numAway(0); // Temporary histos that bin N and pT in dphi Profile1D htmp_num_dphi_2(refData(1, 1, 1)), htmp_num_dphi_5(refData(1, 1, 2)), htmp_num_dphi_30(refData(1, 1, 3)); Profile1D htmp_pt_dphi_2(refData(2, 1, 1)), htmp_pt_dphi_5(refData(2, 1, 2)), htmp_pt_dphi_30(refData(2, 1, 3)); // Final state charged particles /// @todo Non-trackjet track efficiencies are corrected? const Particles& tracks = apply(event, "FS").particles(); for (const Particle& p : tracks) { const double dPhi = deltaPhi(p, jet1); const double pT = p.pT(); if (dPhi < PI/3.0) { ptSumToward += pT; ++numToward; } else if (dPhi < 2*PI/3.0) { ptSumTrans += pT; ++numTrans; // Fill transverse pT distributions if (ptLead/GeV > 2.0) { _ptTrans2->fill(pT/GeV, weight); _totalNumTrans2 += weight; } if (ptLead/GeV > 5.0) { _ptTrans5->fill(pT/GeV, weight); _totalNumTrans5 += weight; } if (ptLead/GeV > 30.0) { _ptTrans30->fill(pT/GeV, weight); _totalNumTrans30 += weight; } } else { ptSumAway += pT; ++numAway; } // Fill tmp histos to bin event's track Nch & pT in dphi const double dPhideg = 180*dPhi/M_PI; if (ptLead/GeV > 2.0) { htmp_num_dphi_2.fill(dPhideg, 1); htmp_pt_dphi_2.fill (dPhideg, pT/GeV); } if (ptLead/GeV > 5.0) { htmp_num_dphi_5.fill(dPhideg, 1); htmp_pt_dphi_5.fill (dPhideg, pT/GeV); } if (ptLead/GeV > 30.0) { htmp_num_dphi_30.fill(dPhideg, 1); htmp_pt_dphi_30.fill (dPhideg, pT/GeV); } } // Update the "proper" dphi profile histograms for (int i = 0; i < 50; i++) { ///< @todo Should really explicitly iterate over nbins for each temp histo if (ptLead/GeV > 2.0) { const double x2 = htmp_pt_dphi_2.bin(i).xMid(); const double num2 = (htmp_num_dphi_2.bin(i).numEntries() > 0) ? htmp_num_dphi_2.bin(i).mean() : 0.0; const double pt2 = (htmp_num_dphi_2.bin(i).numEntries() > 0) ? htmp_pt_dphi_2.bin(i).mean() : 0.0; _numvsDeltaPhi2->fill(x2, num2, weight); _pTvsDeltaPhi2->fill(x2, pt2, weight); } if (ptLead/GeV > 5.0) { const double x5 = htmp_pt_dphi_5.bin(i).xMid(); const double num5 = (htmp_num_dphi_5.bin(i).numEntries() > 0) ? htmp_num_dphi_5.bin(i).mean() : 0.0; const double pt5 = (htmp_num_dphi_5.bin(i).numEntries() > 0) ? htmp_pt_dphi_5.bin(i).mean() : 0.0; _numvsDeltaPhi5->fill(x5, num5, weight); _pTvsDeltaPhi5->fill(x5, pt5, weight); } if (ptLead/GeV > 30.0) { const double x30 = htmp_pt_dphi_30.bin(i).xMid(); const double num30 = (htmp_num_dphi_30.bin(i).numEntries() > 0) ? htmp_num_dphi_30.bin(i).mean() : 0.0; const double pt30 = (htmp_num_dphi_30.bin(i).numEntries() > 0) ? htmp_pt_dphi_30.bin(i).mean() : 0.0; _numvsDeltaPhi30->fill(x30, num30, weight); _pTvsDeltaPhi30->fill(x30, pt30, weight); } } // Log some event details about pT MSG_DEBUG("pT [lead; twd, away, trans] = [" << ptLead << "; " << ptSumToward << ", " << ptSumAway << ", " << ptSumTrans << "]"); // Update the pT profile histograms _ptsumTowardMB->fill(ptLead/GeV, ptSumToward/GeV, weight); _ptsumTowardJ20->fill(ptLead/GeV, ptSumToward/GeV, weight); _ptsumTransMB->fill(ptLead/GeV, ptSumTrans/GeV, weight); _ptsumTransJ20->fill(ptLead/GeV, ptSumTrans/GeV, weight); _ptsumAwayMB->fill(ptLead/GeV, ptSumAway/GeV, weight); _ptsumAwayJ20->fill(ptLead/GeV, ptSumAway/GeV, weight); // Log some event details about Nch MSG_DEBUG("N [twd, away, trans] = [" << ptLead << "; " << numToward << ", " << numTrans << ", " << numAway << "]"); // Update the N_track profile histograms _numTowardMB->fill(ptLead/GeV, numToward, weight); _numTowardJ20->fill(ptLead/GeV, numToward, weight); _numTransMB->fill(ptLead/GeV, numTrans, weight); _numTransJ20->fill(ptLead/GeV, numTrans, weight); _numAwayMB->fill(ptLead/GeV, numAway, weight); _numAwayJ20->fill(ptLead/GeV, numAway, weight); } /// Normalize histos void finalize() { normalize(_ptTrans2, _totalNumTrans2 / _sumWeightsPtLead2); normalize(_ptTrans5, _totalNumTrans5 / _sumWeightsPtLead5); normalize(_ptTrans30, _totalNumTrans30 / _sumWeightsPtLead30); } //@} private: /// Sum total number of charged particles in the trans region, in 3 \f$ p_\perp^\text{lead} \f$ bins. double _totalNumTrans2, _totalNumTrans5, _totalNumTrans30; /// Sum the total number of events in 3 \f$ p_\perp^\text{lead} \f$ bins. double _sumWeightsPtLead2,_sumWeightsPtLead5, _sumWeightsPtLead30; /// @name Histogram collections //@{ // The sumpt vs. dphi and Nch vs. dphi histos Profile1DPtr _numvsDeltaPhi2, _numvsDeltaPhi5, _numvsDeltaPhi30; Profile1DPtr _pTvsDeltaPhi2, _pTvsDeltaPhi5, _pTvsDeltaPhi30; /// Profile histograms, binned in the \f$ p_T \f$ of the leading jet, for /// the \f$ p_T \f$ sum in the toward, transverse and away regions. Profile1DPtr _ptsumTowardMB, _ptsumTransMB, _ptsumAwayMB; Profile1DPtr _ptsumTowardJ20, _ptsumTransJ20, _ptsumAwayJ20; /// Profile histograms, binned in the \f$ p_T \f$ of the leading jet, for /// the number of charged particles per jet in the toward, transverse and /// away regions. Profile1DPtr _numTowardMB, _numTransMB, _numAwayMB; Profile1DPtr _numTowardJ20, _numTransJ20, _numAwayJ20; /// Histogram of \f$ p_T \f$ distribution for 3 different \f$ p_{T1} \f$ IR cutoffs. Histo1DPtr _ptTrans2, _ptTrans5, _ptTrans30; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CDF_2001_S4751469); }