diff --git a/analyses/pluginCDF/CDF_1996_S3108457.cc b/analyses/pluginCDF/CDF_1996_S3108457.cc --- a/analyses/pluginCDF/CDF_1996_S3108457.cc +++ b/analyses/pluginCDF/CDF_1996_S3108457.cc @@ -1,121 +1,110 @@ // -*- 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_S3108457 : public Analysis { public: - /// @name Constructors etc. - //@{ - /// Constructor - CDF_1996_S3108457() - : Analysis("CDF_1996_S3108457") - { } - - //@} + DEFAULT_RIVET_ANALYSIS_CTOR(CDF_1996_S3108457); /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { /// Initialise and register projections here - const FinalState fs((Cuts::etaIn(-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 + // 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_E"); /// Book histograms here, e.g.: for (size_t i=0; i<5; ++i) { - book(_h_m[i] ,1+i, 1, 1); - book(_h_costheta[i] ,10+i, 1, 1); - book(_h_pT[i] ,15+i, 1, 1); + book(_h_m[i], 1+i, 1, 1); + book(_h_costheta[i], 10+i, 1, 1); + book(_h_pT[i], 15+i, 1, 1); } /// @todo Ratios of mass histograms left out: Binning doesn't work out } /// Perform the per-event analysis void analyze(const Event& event) { // Get the smeared jets Jets SJets = apply(event, "SmearedJets_E").jets(Cuts::Et > 20.0*GeV, cmpMomByEt); if (SJets.size() < 2 || SJets.size() > 6) vetoEvent; // Calculate Et, total jet 4 Momentum double sumEt(0), sumE(0); FourMomentum JS(0,0,0,0); - for(const Jet& jet : SJets) { + for (const Jet& jet : SJets) { sumEt += jet.Et()*GeV; sumE += jet.E()*GeV; JS+=jet.momentum(); } if (sumEt < 420*GeV || sumE > 2000*GeV) vetoEvent; double mass = JS.mass(); LorentzTransform cms_boost = LorentzTransform::mkFrameTransformFromBeta(JS.betaVec()); FourMomentum jet0boosted(cms_boost.transform(SJets[0].momentum())); double costheta0 = fabs(cos(jet0boosted.theta())); - if (costheta0 < 2.0/3.0) { - _h_m[SJets.size()-2]->fill(mass); - } + if (costheta0 < 2.0/3.0) _h_m[SJets.size()-2]->fill(mass); if (mass > 600.0*GeV) _h_costheta[JS.size()-2]->fill(costheta0); if (costheta0 < 2.0/3.0 && mass > 600.0*GeV) { - for (const Jet& jet : SJets) { - _h_pT[SJets.size()-2]->fill(jet.pT()); - } + for (const Jet& jet : SJets) _h_pT[SJets.size()-2]->fill(jet.pT()); } } /// Normalise histograms etc., after the run void finalize() { /// Normalise, scale and otherwise manipulate histograms here for (size_t i=0; i<5; ++i) { normalize(_h_m[i], 40.0); normalize(_h_costheta[i], 2.0); normalize(_h_pT[i], 20.0); } } //@} private: /// @name Histograms //@{ Histo1DPtr _h_m[5]; Histo1DPtr _h_costheta[5]; Histo1DPtr _h_pT[5]; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(CDF_1996_S3108457); } 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,455 +1,447 @@ // -*- 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: - /// @name Constructors etc. - //@{ - /// Constructor - CDF_1996_S3349578() - : Analysis("CDF_1996_S3349578") - { - } - - //@} + 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((Cuts::etaIn(-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 + + // 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.: book(_h_3_mNJ ,1, 1, 1); book(_h_3_X3 ,2, 1, 1); book(_h_3_X4 ,3, 1, 1); book(_h_3_costheta3 ,8, 1, 1); book(_h_3_psi3 ,9, 1, 1); book(_h_3_f3 ,14, 1, 1); book(_h_3_f4 ,14, 1, 2); book(_h_3_f5 ,14, 1, 3); book(_h_4_mNJ ,1, 1, 2); book(_h_4_X3 ,4, 1, 1); book(_h_4_X4 ,5, 1, 1); book(_h_4_costheta3 ,10, 1, 1); book(_h_4_psi3 ,11, 1, 1); book(_h_4_f3 ,15, 1, 1); book(_h_4_f4 ,15, 1, 2); book(_h_4_f5 ,15, 1, 3); book(_h_4_XA ,17, 1, 1); book(_h_4_psiAB ,19, 1, 1); book(_h_4_fA ,21, 1, 1); book(_h_4_fB ,21, 1, 2); book(_h_5_mNJ ,1, 1, 3); book(_h_5_X3 ,6, 1, 1); book(_h_5_X4 ,7, 1, 1); book(_h_5_costheta3 ,12, 1, 1); book(_h_5_psi3 ,13, 1, 1); book(_h_5_f3 ,16, 1, 1); book(_h_5_f4 ,16, 1, 2); book(_h_5_f5 ,16, 1, 3); book(_h_5_XA ,18, 1, 1); book(_h_5_XC ,18, 1, 2); book(_h_5_psiAB ,20, 1, 1); book(_h_5_psiCD ,20, 1, 2); book(_h_5_fA ,22, 1, 1); book(_h_5_fB ,23, 1, 1); book(_h_5_fC ,24, 1, 1); book(_h_5_fD ,25, 1, 1); } void analyze(const Event& event) { Jets jets; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); for (const Jet& jet : apply(event, "SmearedJets").jets(Cuts::Et > 20.0*GeV, cmpMomByEt)) { bool separated = true; for (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; } if (jets.size() > 4) { _fiveJetAnalysis(jets); jets.resize(4); } if (jets.size() > 3) { _fourJetAnalysis(jets); jets.resize(3); } if (jets.size() > 2) { _threeJetAnalysis(jets); } } void _threeJetAnalysis(const Jets& jets) { MSG_DEBUG("3 jet analysis"); double sumEt = 0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); for (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; for (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); _h_3_X3->fill(X3); _h_3_X4->fill(X4); _h_3_costheta3->fill(costheta3); _h_3_psi3->fill(psi3); _h_3_f3->fill(f3); _h_3_f4->fill(f4); _h_3_f5->fill(f5); } void _fourJetAnalysis(const Jets& jets) { MSG_DEBUG("4 jet analysis"); double sumEt=0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); for (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; for (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); _h_4_X3->fill(X3); _h_4_X4->fill(X4); _h_4_costheta3->fill(costheta3); _h_4_psi3->fill(psi3); _h_4_f3->fill(f3); _h_4_f4->fill(f4); _h_4_f5->fill(f5); _h_4_XA->fill(XA); _h_4_psiAB->fill(psiAB); _h_4_fA->fill(fA); _h_4_fB->fill(fB); } void _fiveJetAnalysis(const Jets& jets) { MSG_DEBUG("5 jet analysis"); double sumEt=0.0; FourMomentum jetsystem(0.0, 0.0, 0.0, 0.0); for (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; for (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); _h_5_X3->fill(X3); _h_5_X4->fill(X4); _h_5_costheta3->fill(costheta3); _h_5_psi3->fill(psi3); _h_5_f3->fill(f3); _h_5_f4->fill(f4); _h_5_f5->fill(f5); _h_5_XA->fill(XA); _h_5_psiAB->fill(psiAB); _h_5_fA->fill(fA); _h_5_fB->fill(fB); _h_5_XC->fill(XC); _h_5_psiCD->fill(psiCD); _h_5_fC->fill(fC); _h_5_fD->fill(fD); } /// 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,249 +1,248 @@ // -*- 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((Cuts::etaIn(-4.2, 4.2))); + const FinalState fs(Cuts::abseta < 4.2); FastJets fj(fs, FastJets::CDFJETCLU, 0.7); - // declare(fj, "Jets"); // 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"); book(_h_m6J ,1, 1, 1); book(_h_X3ppp ,2, 1, 1); book(_h_X4ppp ,3, 1, 1); book(_h_costheta3ppp ,4, 1, 1); book(_h_psi3ppp ,5, 1, 1); book(_h_f3ppp ,6, 1, 1); book(_h_f4ppp ,6, 1, 2); book(_h_f5ppp ,6, 1, 3); book(_h_XApp ,7, 1, 1); book(_h_XCp ,8, 1, 1); book(_h_XE ,9, 1, 1); book(_h_psiAppBpp ,10, 1, 1); book(_h_psiCpDp ,11, 1, 1); book(_h_psiEF ,12, 1, 1); book(_h_fApp ,13, 1, 1); book(_h_fBpp ,14, 1, 1); book(_h_fCp ,15, 1, 1); book(_h_fDp ,16, 1, 1); book(_h_fE ,17, 1, 1); book(_h_fF ,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; // 3-jet-system variables _h_m6J->fill(m6J); _h_X3ppp->fill(X3ppp); _h_X4ppp->fill(2.0*p4ppp.E()/m6J); _h_costheta3ppp->fill(costheta3ppp); double psi3ppp = _psi(p3ppp, pAV, p4ppp, p5ppp); _h_psi3ppp->fill(psi3ppp); _h_f3ppp->fill(_safeMass(p3ppp)/m6J); _h_f4ppp->fill(_safeMass(p4ppp)/m6J); _h_f5ppp->fill(_safeMass(p5ppp)/m6J); // 4 -> 3 jet variables _h_fApp->fill(_safeMass(pApp)/m6J); _h_fBpp->fill(_safeMass(pApp)/m6J); _h_XApp->fill(pApp.E()/(pApp.E()+pBpp.E())); double psiAppBpp = _psi(pApp, pBpp, pApp+pBpp, pAV); _h_psiAppBpp->fill(psiAppBpp); // 5 -> 4 jet variables _h_fCp->fill(_safeMass(pCp)/m6J); _h_fDp->fill(_safeMass(pDp)/m6J); _h_XCp->fill(pCp.E()/(pCp.E()+pDp.E())); double psiCpDp = _psi(pCp, pDp, pCp+pDp, pAV); _h_psiCpDp->fill(psiCpDp); // 6 -> 5 jet variables _h_fE->fill(_safeMass(pE)/m6J); _h_fF->fill(_safeMass(pF)/m6J); _h_XE->fill(pE.E()/(pE.E()+pF.E())); double psiEF = _psi(pE, pF, pE+pF, pAV); _h_psiEF->fill(psiEF); } 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,264 +1,264 @@ // -*- 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") { } /// @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((Cuts::etaIn(-1.0, 1.0) && Cuts::pT >= 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"); book(_numvsDeltaPhi2 ,1, 1, 1); book(_numvsDeltaPhi5 ,1, 1, 2); book(_numvsDeltaPhi30 ,1, 1, 3); book(_pTvsDeltaPhi2 ,2, 1, 1); book(_pTvsDeltaPhi5 ,2, 1, 2); book(_pTvsDeltaPhi30 ,2, 1, 3); book(_numTowardMB ,3, 1, 1); book(_numTransMB ,3, 1, 2); book(_numAwayMB ,3, 1, 3); book(_numTowardJ20 ,4, 1, 1); book(_numTransJ20 ,4, 1, 2); book(_numAwayJ20 ,4, 1, 3); book(_ptsumTowardMB ,5, 1, 1); book(_ptsumTransMB ,5, 1, 2); book(_ptsumAwayMB ,5, 1, 3); book(_ptsumTowardJ20 ,6, 1, 1); book(_ptsumTransJ20 ,6, 1, 2); book(_ptsumAwayJ20 ,6, 1, 3); book(_ptTrans2 ,7, 1, 1); book(_ptTrans5 ,7, 1, 2); book(_ptTrans30 ,7, 1, 3); book(_totalNumTrans2, "totalNumTrans2"); book(_totalNumTrans5, "totalNumTrans5"); book(_totalNumTrans30, "totalNumTrans30"); book(_sumWeightsPtLead2, "sumWeightsPtLead2"); book(_sumWeightsPtLead5, "sumWeightsPtLead5"); book(_sumWeightsPtLead30, "sumWeightsPtLead30"); } /// 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 if (ptLead/GeV > 2.0) _sumWeightsPtLead2->fill(); if (ptLead/GeV > 5.0) _sumWeightsPtLead5->fill(); if (ptLead/GeV > 30.0) _sumWeightsPtLead30->fill(); // 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); _totalNumTrans2->fill(); } if (ptLead/GeV > 5.0) { _ptTrans5->fill(pT/GeV); _totalNumTrans5->fill(); } if (ptLead/GeV > 30.0) { _ptTrans30->fill(pT/GeV); _totalNumTrans30->fill(); } } 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); _pTvsDeltaPhi2->fill(x2, pt2); } 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); _pTvsDeltaPhi5->fill(x5, pt5); } 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); _pTvsDeltaPhi30->fill(x30, pt30); } } // 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); _ptsumTowardJ20->fill(ptLead/GeV, ptSumToward/GeV); _ptsumTransMB->fill(ptLead/GeV, ptSumTrans/GeV); _ptsumTransJ20->fill(ptLead/GeV, ptSumTrans/GeV); _ptsumAwayMB->fill(ptLead/GeV, ptSumAway/GeV); _ptsumAwayJ20->fill(ptLead/GeV, ptSumAway/GeV); // 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); _numTowardJ20->fill(ptLead/GeV, numToward); _numTransMB->fill(ptLead/GeV, numTrans); _numTransJ20->fill(ptLead/GeV, numTrans); _numAwayMB->fill(ptLead/GeV, numAway); _numAwayJ20->fill(ptLead/GeV, numAway); } /// 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. CounterPtr _totalNumTrans2, _totalNumTrans5, _totalNumTrans30; /// Sum the total number of events in 3 \f$ p_\perp^\text{lead} \f$ bins. CounterPtr _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); } diff --git a/analyses/pluginRHIC/STAR_2008_S7869363.cc b/analyses/pluginRHIC/STAR_2008_S7869363.cc --- a/analyses/pluginRHIC/STAR_2008_S7869363.cc +++ b/analyses/pluginRHIC/STAR_2008_S7869363.cc @@ -1,173 +1,170 @@ // -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/LossyFinalState.hh" namespace Rivet { /// @todo Replace with SmearedParticles class STARRandomFilter { public: - STARRandomFilter() { } - // Return true to throw away a particle bool operator()(const Particle& p) { - /// @todo Use a better RNG? size_t idx = int(floor(p.pT()/MeV/50)); if (idx > 11) idx = 11; - return (rand()/static_cast(RAND_MAX) > _trkeff[idx]); + return (rand01() > _trkeff[idx]); } CmpState compare(const STARRandomFilter& other) const { - return CmpState::GT; // @todo really? + return CmpState::EQ; } private: const static double _trkeff[12]; }; // Here we have the track reconstruction efficiencies for tracks with pT from 0 to 600 MeV // in steps of 50 MeV. The efficiency is assumed to be 0.88 for pT >= 600 MeV const double STARRandomFilter::_trkeff[12] = {0,0,0.38,0.72,0.78,0.81,0.82,0.84,0.85,0.86,0.87,0.88}; class STAR_2008_S7869363 : public Analysis { public: /// @name Constructors etc. //@{ /// Constructor STAR_2008_S7869363() : Analysis("STAR_2008_S7869363") { } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { - const ChargedFinalState cfs((Cuts::etaIn(-0.5, 0.5) && Cuts::pT >= 0.2*GeV)); + const ChargedFinalState cfs(Cuts::abseta < 0.5 && Cuts::pT > 0.2*GeV); const LossyFinalState lfs(cfs, STARRandomFilter()); declare(lfs, "FS"); book(_h_dNch ,1, 1, 1); book(_h_dpT_Pi ,2, 1, 1); book(_h_dpT_Piplus ,2, 1, 2); book(_h_dpT_Kaon ,2, 1, 3); book(_h_dpT_Kaonplus ,2, 1, 4); book(_h_dpT_AntiProton ,2, 1, 5); book(_h_dpT_Proton ,2, 1, 6); // book(nCutsPassed, "nCutsPassed"); // book(nPi, "nPi"); // book(nPiPlus, "nPiPlus"); // book(nKaon, "nKaon"); // book(nKaonPlus, "nKaonPlus"); // book(nProton, "nProton"); // book(nAntiProton, "nAntiProton"); } /// Perform the per-event analysis void analyze(const Event& event) { const FinalState& charged = apply(event, "FS"); // Vertex reconstruction efficiencies as a function of charged multiplicity. // For events with more than 23 reconstructed tracks the efficiency is 100%. double vtxeffs[24] = { 0.000000,0.512667,0.739365,0.847131,0.906946,0.940922,0.959328,0.96997, 0.975838,0.984432,0.988311,0.990327,0.990758,0.995767,0.99412,0.992271, 0.996631,0.994802,0.99635,0.997384,0.998986,0.996441,0.994513,1.000000 }; double vtxeff = 1.0; if (charged.particles().size() < 24) { vtxeff = vtxeffs[charged.particles().size()]; } const double weight = vtxeff; for (const Particle& p : charged.particles()) { double pT = p.pT()/GeV; double y = p.rapidity(); if (fabs(y) < 0.1) { // nCutsPassed->fill(weight); const PdgId id = p.pid(); switch (id) { case -211: _h_dpT_Pi->fill(pT, weight/(TWOPI*pT*0.2)); // nPi->fill(weight); break; case 211: _h_dpT_Piplus->fill(pT, weight/(TWOPI*pT*0.2)); // nPiPlus->fill(weight); break; case -321: _h_dpT_Kaon->fill(pT, weight/(TWOPI*pT*0.2)); // nKaon->fill(weight); break; case 321: _h_dpT_Kaonplus->fill(pT, weight/(TWOPI*pT*0.2)); // nKaonPlus->fill(weight); break; case -2212: _h_dpT_AntiProton->fill(pT, weight/(TWOPI*pT*0.2)); // nAntiProton->fill(weight); break; case 2212: _h_dpT_Proton->fill(pT, weight/(TWOPI*pT*0.2)); // nProton->fill(weight); break; } } else { continue; } } _h_dNch->fill(charged.particles().size(), weight); } /// Normalise histograms etc., after the run void finalize() { //double nTot = nPi + nPiPlus + nKaon + nKaonPlus + nProton + nAntiProton; normalize(_h_dNch); /// @todo Norm to data! normalize(_h_dpT_Pi , 0.389825 ); normalize(_h_dpT_Piplus , 0.396025 ); normalize(_h_dpT_Kaon , 0.03897 ); normalize(_h_dpT_Kaonplus , 0.04046 ); normalize(_h_dpT_AntiProton, 0.0187255); normalize(_h_dpT_Proton , 0.016511 ); } private: Histo1DPtr _h_dNch; Histo1DPtr _h_dpT_Pi, _h_dpT_Piplus; Histo1DPtr _h_dpT_Kaon, _h_dpT_Kaonplus; Histo1DPtr _h_dpT_AntiProton, _h_dpT_Proton; Profile1DPtr _h_pT_vs_Nch; //CounterPtr nCutsPassed, nPi, nPiPlus, nKaon, nKaonPlus, nProton, nAntiProton; }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(STAR_2008_S7869363); }