|
[Rivet-svn] r1785 - in trunk: . bin include/Rivet include/Rivet/Analyses src/Analysesblackhole at projects.hepforge.org blackhole at projects.hepforge.orgMon Aug 31 09:50:20 BST 2009
Author: buckley Date: Mon Aug 31 09:50:19 2009 New Revision: 1785 Log: Removing/merging headers for all D0 analyses Deleted: trunk/include/Rivet/Analyses/D0_1996_S3214044.hh trunk/include/Rivet/Analyses/D0_1996_S3324664.hh trunk/include/Rivet/Analyses/D0_2001_S4674421.hh trunk/include/Rivet/Analyses/D0_2004_S5992206.hh trunk/include/Rivet/Analyses/D0_2006_S6438750.hh trunk/include/Rivet/Analyses/D0_2007_S7075677.hh trunk/include/Rivet/Analyses/D0_2008_S6879055.hh trunk/include/Rivet/Analyses/D0_2008_S7554427.hh trunk/include/Rivet/Analyses/D0_2008_S7662670.hh trunk/include/Rivet/Analyses/D0_2008_S7719523.hh trunk/include/Rivet/Analyses/D0_2008_S7837160.hh trunk/include/Rivet/Analyses/D0_2008_S7863608.hh trunk/include/Rivet/Analyses/D0_2009_S8202443.hh trunk/include/Rivet/Analyses/D0_2009_S8320160.hh trunk/include/Rivet/Analyses/D0_2009_S8349509.hh Modified: trunk/ChangeLog trunk/bin/rivet-mkanalysis trunk/include/Rivet/Makefile.am trunk/src/Analyses/D0_1996_S3214044.cc trunk/src/Analyses/D0_1996_S3324664.cc trunk/src/Analyses/D0_2001_S4674421.cc trunk/src/Analyses/D0_2004_S5992206.cc trunk/src/Analyses/D0_2006_S6438750.cc trunk/src/Analyses/D0_2007_S7075677.cc trunk/src/Analyses/D0_2008_S6879055.cc trunk/src/Analyses/D0_2008_S7554427.cc trunk/src/Analyses/D0_2008_S7662670.cc trunk/src/Analyses/D0_2008_S7719523.cc trunk/src/Analyses/D0_2008_S7837160.cc trunk/src/Analyses/D0_2008_S7863608.cc trunk/src/Analyses/D0_2009_S8202443.cc trunk/src/Analyses/D0_2009_S8320160.cc trunk/src/Analyses/D0_2009_S8349509.cc Modified: trunk/ChangeLog ============================================================================== --- trunk/ChangeLog Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/ChangeLog Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,7 @@ 2009-08-31 Andy Buckley <andy at insectnation.org> + * Removing headers for D0 analyses. + * Exit with an error message if addProjection is used twice from the same parent with distinct projections. Modified: trunk/bin/rivet-mkanalysis ============================================================================== --- trunk/bin/rivet-mkanalysis Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/bin/rivet-mkanalysis Mon Aug 31 09:50:19 2009 (r1785) @@ -122,6 +122,7 @@ /// Perform the per-event analysis void analyze(const Event& event) { + const double weight = event.weight(); /// @todo Do the event by event analysis here Modified: trunk/include/Rivet/Makefile.am ============================================================================== --- trunk/include/Rivet/Makefile.am Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/include/Rivet/Makefile.am Mon Aug 31 09:50:19 2009 (r1785) @@ -29,26 +29,9 @@ ## Standard analyses -# nobase_pkginclude_HEADERS += \ -# Analyses/ExampleAnalysis.hh \ -# Analyses/ExampleTree.hh nobase_dist_noinst_HEADERS += \ Analyses/ExampleAnalysis.hh \ Analyses/ExampleTree.hh \ - Analyses/D0_1996_S3214044.hh \ - Analyses/D0_1996_S3324664.hh \ - Analyses/D0_2001_S4674421.hh \ - Analyses/D0_2004_S5992206.hh \ - Analyses/D0_2006_S6438750.hh \ - Analyses/D0_2007_S7075677.hh \ - Analyses/D0_2008_S6879055.hh \ - Analyses/D0_2008_S7554427.hh \ - Analyses/D0_2008_S7662670.hh \ - Analyses/D0_2008_S7719523.hh \ - Analyses/D0_2008_S7837160.hh \ - Analyses/D0_2008_S7863608.hh \ - Analyses/D0_2009_S8202443.hh \ - Analyses/D0_2009_S8320160.hh \ Analyses/DELPHI_1995_S3137023.hh \ Analyses/DELPHI_1996_S3430090.hh \ Analyses/DELPHI_2002_069_CONF_603.hh \ Modified: trunk/src/Analyses/D0_1996_S3214044.cc ============================================================================== --- trunk/src/Analyses/D0_1996_S3214044.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_1996_S3214044.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,8 +1,7 @@ // -*- C++ -*- -#include "Rivet/Rivet.hh" +#include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" -#include "Rivet/Analyses/D0_1996_S3214044.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Math/LorentzTrans.hh" @@ -12,209 +11,262 @@ namespace Rivet { - D0_1996_S3214044::D0_1996_S3214044() : Analysis("D0_1996_S3214044") { - setBeams(PROTON, ANTIPROTON); - setNeedsCrossSection(false); - - const FinalState fs(-4.5, 4.5); - addProjection(fs, "FS"); - /// @todo Use correct jet algorithm - addProjection(FastJets(fs, FastJets::D0ILCONE, 0.7, 20.0*GeV), "ConeJets"); - } - - - void D0_1996_S3214044::init() { - - _h_3j_x3 = bookHistogram1D(1, 1, 1); - _h_3j_x5 = bookHistogram1D(2, 1, 1); - _h_3j_costheta3 = bookHistogram1D(3, 1, 1); - _h_3j_psi = bookHistogram1D(4, 1, 1); - _h_3j_mu34 = bookHistogram1D(5, 1, 1); - _h_3j_mu35 = bookHistogram1D(6, 1, 1); - _h_3j_mu45 = bookHistogram1D(7, 1, 1); - - _h_4j_x3 = bookHistogram1D(8, 1, 1); - _h_4j_x4 = bookHistogram1D(9, 1, 1); - _h_4j_x5 = bookHistogram1D(10, 1, 1); - _h_4j_x6 = bookHistogram1D(11, 1, 1); - _h_4j_costheta3 = bookHistogram1D(12, 1, 1); - _h_4j_costheta4 = bookHistogram1D(13, 1, 1); - _h_4j_costheta5 = bookHistogram1D(14, 1, 1); - _h_4j_costheta6 = bookHistogram1D(15, 1, 1); - _h_4j_cosomega34 = bookHistogram1D(16, 1, 1); - _h_4j_cosomega35 = bookHistogram1D(17, 1, 1); - _h_4j_cosomega36 = bookHistogram1D(18, 1, 1); - _h_4j_cosomega45 = bookHistogram1D(19, 1, 1); - _h_4j_cosomega46 = bookHistogram1D(20, 1, 1); - _h_4j_cosomega56 = bookHistogram1D(21, 1, 1); - _h_4j_mu34 = bookHistogram1D(22, 1, 1); - _h_4j_mu35 = bookHistogram1D(23, 1, 1); - _h_4j_mu36 = bookHistogram1D(24, 1, 1); - _h_4j_mu45 = bookHistogram1D(25, 1, 1); - _h_4j_mu46 = bookHistogram1D(26, 1, 1); - _h_4j_mu56 = bookHistogram1D(27, 1, 1); - _h_4j_theta_BZ = bookHistogram1D(28, 1, 1); - _h_4j_costheta_NR = bookHistogram1D(29, 1, 1); - - } - - - void D0_1996_S3214044::analyze(const Event& event) { - - const double weight = event.weight(); - - Jets jets_in; - foreach (const Jet& jet, applyProjection<FastJets>(event, "ConeJets").jetsByEt()) { - if (fabs(jet.momentum().eta())<3.0) { - jets_in.push_back(jet); - } - } - - Jets jets_isolated; - for (size_t i=0; i<jets_in.size(); ++i) { - bool isolated=true; - for (size_t j=0; j<jets_in.size(); ++j) { - if (i!=j && deltaR(jets_in[i].momentum(), jets_in[j].momentum())<1.4) { - isolated=false; - break; + class D0_1996_S3214044 : public Analysis { + public: + + /// @name Constructors etc. + //@{ + + /// Constructor + D0_1996_S3214044() + : Analysis("D0_1996_S3214044") + { + setBeams(PROTON, ANTIPROTON); + setNeedsCrossSection(false); + + const FinalState fs(-4.5, 4.5); + addProjection(fs, "FS"); + /// @todo Use correct jet algorithm + addProjection(FastJets(fs, FastJets::D0ILCONE, 0.7, 20.0*GeV), "ConeJets"); + } + + + /// @name Analysis methods + //@{ + + /// Book histograms + void init() { + + _h_3j_x3 = bookHistogram1D(1, 1, 1); + _h_3j_x5 = bookHistogram1D(2, 1, 1); + _h_3j_costheta3 = bookHistogram1D(3, 1, 1); + _h_3j_psi = bookHistogram1D(4, 1, 1); + _h_3j_mu34 = bookHistogram1D(5, 1, 1); + _h_3j_mu35 = bookHistogram1D(6, 1, 1); + _h_3j_mu45 = bookHistogram1D(7, 1, 1); + + _h_4j_x3 = bookHistogram1D(8, 1, 1); + _h_4j_x4 = bookHistogram1D(9, 1, 1); + _h_4j_x5 = bookHistogram1D(10, 1, 1); + _h_4j_x6 = bookHistogram1D(11, 1, 1); + _h_4j_costheta3 = bookHistogram1D(12, 1, 1); + _h_4j_costheta4 = bookHistogram1D(13, 1, 1); + _h_4j_costheta5 = bookHistogram1D(14, 1, 1); + _h_4j_costheta6 = bookHistogram1D(15, 1, 1); + _h_4j_cosomega34 = bookHistogram1D(16, 1, 1); + _h_4j_cosomega35 = bookHistogram1D(17, 1, 1); + _h_4j_cosomega36 = bookHistogram1D(18, 1, 1); + _h_4j_cosomega45 = bookHistogram1D(19, 1, 1); + _h_4j_cosomega46 = bookHistogram1D(20, 1, 1); + _h_4j_cosomega56 = bookHistogram1D(21, 1, 1); + _h_4j_mu34 = bookHistogram1D(22, 1, 1); + _h_4j_mu35 = bookHistogram1D(23, 1, 1); + _h_4j_mu36 = bookHistogram1D(24, 1, 1); + _h_4j_mu45 = bookHistogram1D(25, 1, 1); + _h_4j_mu46 = bookHistogram1D(26, 1, 1); + _h_4j_mu56 = bookHistogram1D(27, 1, 1); + _h_4j_theta_BZ = bookHistogram1D(28, 1, 1); + _h_4j_costheta_NR = bookHistogram1D(29, 1, 1); + + } + + + void analyze(const Event& event) { + const double weight = event.weight(); + + Jets jets_in; + foreach (const Jet& jet, applyProjection<FastJets>(event, "ConeJets").jetsByEt()) { + if (fabs(jet.momentum().eta()) < 3.0) { + jets_in.push_back(jet); + } + } + + Jets jets_isolated; + for (size_t i = 0; i < jets_in.size(); ++i) { + bool isolated=true; + for (size_t j = 0; j < jets_in.size(); ++j) { + if (i != j && deltaR(jets_in[i].momentum(), jets_in[j].momentum()) < 1.4) { + isolated = false; + break; + } + } + if (isolated) { + jets_isolated.push_back(jets_in[i]); } } - if (isolated) { - jets_isolated.push_back(jets_in[i]); + + if (jets_isolated.size() == 0 || jets_isolated[0].momentum().Et() < 60.0*GeV) { + vetoEvent; } + + if (jets_isolated.size() > 2) _threeJetAnalysis(jets_isolated, weight); + if (jets_isolated.size() > 3) _fourJetAnalysis(jets_isolated, weight); + } + + + void finalize() { + normalize(_h_3j_x3, 1.0); + normalize(_h_3j_x5, 1.0); + normalize(_h_3j_costheta3, 1.0); + normalize(_h_3j_psi, 1.0); + normalize(_h_3j_mu34, 1.0); + normalize(_h_3j_mu35, 1.0); + normalize(_h_3j_mu45, 1.0); + normalize(_h_4j_x3, 1.0); + normalize(_h_4j_x4, 1.0); + normalize(_h_4j_x5, 1.0); + normalize(_h_4j_x6, 1.0); + normalize(_h_4j_costheta3, 1.0); + normalize(_h_4j_costheta4, 1.0); + normalize(_h_4j_costheta5, 1.0); + normalize(_h_4j_costheta6, 1.0); + normalize(_h_4j_cosomega34, 1.0); + normalize(_h_4j_cosomega35, 1.0); + normalize(_h_4j_cosomega36, 1.0); + normalize(_h_4j_cosomega45, 1.0); + normalize(_h_4j_cosomega46, 1.0); + normalize(_h_4j_cosomega56, 1.0); + normalize(_h_4j_mu34, 1.0); + normalize(_h_4j_mu35, 1.0); + normalize(_h_4j_mu36, 1.0); + normalize(_h_4j_mu45, 1.0); + normalize(_h_4j_mu46, 1.0); + normalize(_h_4j_mu56, 1.0); + normalize(_h_4j_theta_BZ, 1.0); + normalize(_h_4j_costheta_NR, 1.0); } - if (jets_isolated.size()==0 || jets_isolated[0].momentum().Et()<60.0*GeV) { - vetoEvent; + //@} + + + private: + + /// @name Helper functions + //@{ + + void _threeJetAnalysis(const Jets& jets, const double& weight) { + // >=3 jet events + FourMomentum jjj(jets[0].momentum()+jets[1].momentum()+jets[2].momentum()); + const double sqrts = jjj.mass(); + if (sqrts<200*GeV) { + return; + } + + LorentzTransform cms_boost(-jjj.boostVector()); + vector<FourMomentum> jets_boosted; + foreach (Jet jet, jets) { + jets_boosted.push_back(cms_boost.transform(jet.momentum())); + } + std::sort(jets_boosted.begin(), jets_boosted.end(), FourMomentum::byEDescending()); + FourMomentum p3(jets_boosted[0]); + FourMomentum p4(jets_boosted[1]); + FourMomentum p5(jets_boosted[2]); + + Vector3 beam1(0.0, 0.0, 1.0); + Vector3 p1xp3 = beam1.cross(p3.vector3()); + Vector3 p4xp5 = p4.vector3().cross(p5.vector3()); + const double cospsi = p1xp3.dot(p4xp5)/p1xp3.mod()/p4xp5.mod(); + + _h_3j_x3->fill(2.0*p3.E()/sqrts, weight); + _h_3j_x5->fill(2.0*p5.E()/sqrts, weight); + _h_3j_costheta3->fill(fabs(cos(p3.theta())), weight); + _h_3j_psi->fill(acos(cospsi)/degree, weight); + _h_3j_mu34->fill(FourMomentum(p3+p4).mass()/sqrts, weight); + _h_3j_mu35->fill(FourMomentum(p3+p5).mass()/sqrts, weight); + _h_3j_mu45->fill(FourMomentum(p4+p5).mass()/sqrts, weight); } - if (jets_isolated.size()>2) threeJetAnalysis(jets_isolated, weight); - if (jets_isolated.size()>3) fourJetAnalysis(jets_isolated, weight); - } - - - void D0_1996_S3214044::threeJetAnalysis(const Jets& jets, const double& weight) { - // >=3 jet events - FourMomentum jjj(jets[0].momentum()+jets[1].momentum()+jets[2].momentum()); - const double sqrts = jjj.mass(); - if (sqrts<200*GeV) { - return; - } - - LorentzTransform cms_boost(-jjj.boostVector()); - std::vector<FourMomentum> jets_boosted; - foreach (Jet jet, jets) { - jets_boosted.push_back(cms_boost.transform(jet.momentum())); - } - std::sort(jets_boosted.begin(), jets_boosted.end(), FourMomentum::byEDescending()); - FourMomentum p3(jets_boosted[0]); - FourMomentum p4(jets_boosted[1]); - FourMomentum p5(jets_boosted[2]); - - Vector3 beam1(0.0, 0.0, 1.0); - Vector3 p1xp3 = beam1.cross(p3.vector3()); - Vector3 p4xp5 = p4.vector3().cross(p5.vector3()); - const double cospsi = p1xp3.dot(p4xp5)/p1xp3.mod()/p4xp5.mod(); - - _h_3j_x3->fill(2.0*p3.E()/sqrts, weight); - _h_3j_x5->fill(2.0*p5.E()/sqrts, weight); - _h_3j_costheta3->fill(fabs(cos(p3.theta())), weight); - _h_3j_psi->fill(acos(cospsi)/degree, weight); - _h_3j_mu34->fill(FourMomentum(p3+p4).mass()/sqrts, weight); - _h_3j_mu35->fill(FourMomentum(p3+p5).mass()/sqrts, weight); - _h_3j_mu45->fill(FourMomentum(p4+p5).mass()/sqrts, weight); - } - - - void D0_1996_S3214044::fourJetAnalysis(const Jets& jets, const double& weight) { - // >=4 jet events - FourMomentum jjjj(jets[0].momentum()+jets[1].momentum()+jets[2].momentum()+ - jets[3].momentum()); - const double sqrts = jjjj.mass(); - if (sqrts<200*GeV) { - return; - } - - LorentzTransform cms_boost(-jjjj.boostVector()); - std::vector<FourMomentum> jets_boosted; - foreach (Jet jet, jets) { - jets_boosted.push_back(cms_boost.transform(jet.momentum())); - } - sort(jets_boosted.begin(), jets_boosted.end(), FourMomentum::byEDescending()); - FourMomentum p3(jets_boosted[0]); - FourMomentum p4(jets_boosted[1]); - FourMomentum p5(jets_boosted[2]); - FourMomentum p6(jets_boosted[3]); - - Vector3 p3xp4 = p3.vector3().cross(p4.vector3()); - Vector3 p5xp6 = p5.vector3().cross(p6.vector3()); - const double costheta_BZ = p3xp4.dot(p5xp6)/p3xp4.mod()/p5xp6.mod(); - const double costheta_NR = (p3.vector3()-p4.vector3()).dot(p5.vector3()-p6.vector3())/ - (p3.vector3()-p4.vector3()).mod()/(p5.vector3()-p6.vector3()).mod(); - - _h_4j_x3->fill(2.0*p3.E()/sqrts, weight); - _h_4j_x4->fill(2.0*p4.E()/sqrts, weight); - _h_4j_x5->fill(2.0*p5.E()/sqrts, weight); - _h_4j_x6->fill(2.0*p6.E()/sqrts, weight); - _h_4j_costheta3->fill(fabs(cos(p3.theta())), weight); - _h_4j_costheta4->fill(fabs(cos(p4.theta())), weight); - _h_4j_costheta5->fill(fabs(cos(p5.theta())), weight); - _h_4j_costheta6->fill(fabs(cos(p6.theta())), weight); - _h_4j_cosomega34->fill(cos(p3.angle(p4)), weight); - _h_4j_cosomega35->fill(cos(p3.angle(p5)), weight); - _h_4j_cosomega36->fill(cos(p3.angle(p6)), weight); - _h_4j_cosomega45->fill(cos(p4.angle(p5)), weight); - _h_4j_cosomega46->fill(cos(p4.angle(p6)), weight); - _h_4j_cosomega56->fill(cos(p5.angle(p6)), weight); - _h_4j_mu34->fill(FourMomentum(p3+p4).mass()/sqrts, weight); - _h_4j_mu35->fill(FourMomentum(p3+p5).mass()/sqrts, weight); - _h_4j_mu36->fill(FourMomentum(p3+p6).mass()/sqrts, weight); - _h_4j_mu45->fill(FourMomentum(p4+p5).mass()/sqrts, weight); - _h_4j_mu46->fill(FourMomentum(p4+p6).mass()/sqrts, weight); - _h_4j_mu56->fill(FourMomentum(p5+p6).mass()/sqrts, weight); - _h_4j_theta_BZ->fill(acos(costheta_BZ)/degree, weight); - _h_4j_costheta_NR->fill(costheta_NR, weight); - - } - - - void D0_1996_S3214044::finalize() { - - normalize(_h_3j_x3, 1.0); - normalize(_h_3j_x5, 1.0); - normalize(_h_3j_costheta3, 1.0); - normalize(_h_3j_psi, 1.0); - normalize(_h_3j_mu34, 1.0); - normalize(_h_3j_mu35, 1.0); - normalize(_h_3j_mu45, 1.0); - normalize(_h_4j_x3, 1.0); - normalize(_h_4j_x4, 1.0); - normalize(_h_4j_x5, 1.0); - normalize(_h_4j_x6, 1.0); - normalize(_h_4j_costheta3, 1.0); - normalize(_h_4j_costheta4, 1.0); - normalize(_h_4j_costheta5, 1.0); - normalize(_h_4j_costheta6, 1.0); - normalize(_h_4j_cosomega34, 1.0); - normalize(_h_4j_cosomega35, 1.0); - normalize(_h_4j_cosomega36, 1.0); - normalize(_h_4j_cosomega45, 1.0); - normalize(_h_4j_cosomega46, 1.0); - normalize(_h_4j_cosomega56, 1.0); - normalize(_h_4j_mu34, 1.0); - normalize(_h_4j_mu35, 1.0); - normalize(_h_4j_mu36, 1.0); - normalize(_h_4j_mu45, 1.0); - normalize(_h_4j_mu46, 1.0); - normalize(_h_4j_mu56, 1.0); - normalize(_h_4j_theta_BZ, 1.0); - normalize(_h_4j_costheta_NR, 1.0); + + void _fourJetAnalysis(const Jets& jets, const double& weight) { + // >=4 jet events + FourMomentum jjjj(jets[0].momentum() + jets[1].momentum() + jets[2].momentum()+ jets[3].momentum()); + const double sqrts = jjjj.mass(); + if (sqrts < 200*GeV) return; + + LorentzTransform cms_boost(-jjjj.boostVector()); + vector<FourMomentum> jets_boosted; + foreach (Jet jet, jets) { + jets_boosted.push_back(cms_boost.transform(jet.momentum())); + } + sort(jets_boosted.begin(), jets_boosted.end(), FourMomentum::byEDescending()); + FourMomentum p3(jets_boosted[0]); + FourMomentum p4(jets_boosted[1]); + FourMomentum p5(jets_boosted[2]); + FourMomentum p6(jets_boosted[3]); + + Vector3 p3xp4 = p3.vector3().cross(p4.vector3()); + Vector3 p5xp6 = p5.vector3().cross(p6.vector3()); + const double costheta_BZ = p3xp4.dot(p5xp6)/p3xp4.mod()/p5xp6.mod(); + const double costheta_NR = (p3.vector3()-p4.vector3()).dot(p5.vector3()-p6.vector3())/ + (p3.vector3()-p4.vector3()).mod()/(p5.vector3()-p6.vector3()).mod(); + + _h_4j_x3->fill(2.0*p3.E()/sqrts, weight); + _h_4j_x4->fill(2.0*p4.E()/sqrts, weight); + _h_4j_x5->fill(2.0*p5.E()/sqrts, weight); + _h_4j_x6->fill(2.0*p6.E()/sqrts, weight); + _h_4j_costheta3->fill(fabs(cos(p3.theta())), weight); + _h_4j_costheta4->fill(fabs(cos(p4.theta())), weight); + _h_4j_costheta5->fill(fabs(cos(p5.theta())), weight); + _h_4j_costheta6->fill(fabs(cos(p6.theta())), weight); + _h_4j_cosomega34->fill(cos(p3.angle(p4)), weight); + _h_4j_cosomega35->fill(cos(p3.angle(p5)), weight); + _h_4j_cosomega36->fill(cos(p3.angle(p6)), weight); + _h_4j_cosomega45->fill(cos(p4.angle(p5)), weight); + _h_4j_cosomega46->fill(cos(p4.angle(p6)), weight); + _h_4j_cosomega56->fill(cos(p5.angle(p6)), weight); + _h_4j_mu34->fill(FourMomentum(p3+p4).mass()/sqrts, weight); + _h_4j_mu35->fill(FourMomentum(p3+p5).mass()/sqrts, weight); + _h_4j_mu36->fill(FourMomentum(p3+p6).mass()/sqrts, weight); + _h_4j_mu45->fill(FourMomentum(p4+p5).mass()/sqrts, weight); + _h_4j_mu46->fill(FourMomentum(p4+p6).mass()/sqrts, weight); + _h_4j_mu56->fill(FourMomentum(p5+p6).mass()/sqrts, weight); + _h_4j_theta_BZ->fill(acos(costheta_BZ)/degree, weight); + _h_4j_costheta_NR->fill(costheta_NR, weight); + + } + - } + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D *_h_3j_x3; + AIDA::IHistogram1D *_h_3j_x5; + AIDA::IHistogram1D *_h_3j_costheta3; + AIDA::IHistogram1D *_h_3j_psi; + AIDA::IHistogram1D *_h_3j_mu34; + AIDA::IHistogram1D *_h_3j_mu35; + AIDA::IHistogram1D *_h_3j_mu45; + + AIDA::IHistogram1D *_h_4j_x3; + AIDA::IHistogram1D *_h_4j_x4; + AIDA::IHistogram1D *_h_4j_x5; + AIDA::IHistogram1D *_h_4j_x6; + AIDA::IHistogram1D *_h_4j_costheta3; + AIDA::IHistogram1D *_h_4j_costheta4; + AIDA::IHistogram1D *_h_4j_costheta5; + AIDA::IHistogram1D *_h_4j_costheta6; + AIDA::IHistogram1D *_h_4j_cosomega34; + AIDA::IHistogram1D *_h_4j_cosomega35; + AIDA::IHistogram1D *_h_4j_cosomega36; + AIDA::IHistogram1D *_h_4j_cosomega45; + AIDA::IHistogram1D *_h_4j_cosomega46; + AIDA::IHistogram1D *_h_4j_cosomega56; + AIDA::IHistogram1D *_h_4j_mu34; + AIDA::IHistogram1D *_h_4j_mu35; + AIDA::IHistogram1D *_h_4j_mu36; + AIDA::IHistogram1D *_h_4j_mu45; + AIDA::IHistogram1D *_h_4j_mu46; + AIDA::IHistogram1D *_h_4j_mu56; + AIDA::IHistogram1D *_h_4j_theta_BZ; + AIDA::IHistogram1D *_h_4j_costheta_NR; + //@} + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_1996_S3214044> plugin_D0_1996_S3214044; - + } Modified: trunk/src/Analyses/D0_1996_S3324664.cc ============================================================================== --- trunk/src/Analyses/D0_1996_S3324664.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_1996_S3324664.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,95 +1,119 @@ // -*- C++ -*- -#include "Rivet/Rivet.hh" +#include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" -#include "Rivet/Analyses/D0_1996_S3324664.hh" +#include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/FinalState.hh" namespace Rivet { - D0_1996_S3324664::D0_1996_S3324664() : Analysis("D0_1996_S3324664") { - setBeams(PROTON, ANTIPROTON); - setNeedsCrossSection(false); + class D0_1996_S3324664 : public Analysis { + public: - const FinalState fs(-3.2, 3.2); - addProjection(fs, "FS"); - /// @todo Use correct jet algorithm - addProjection(FastJets(fs, FastJets::D0ILCONE, 0.7, 20.0*GeV), "ConeJets"); - } + /// @name Constructors etc. + //@{ + /// Constructor + D0_1996_S3324664() : Analysis("D0_1996_S3324664") { + setBeams(PROTON, ANTIPROTON); + setNeedsCrossSection(false); + + const FinalState fs(-3.2, 3.2); + addProjection(fs, "FS"); + /// @todo Use correct jet algorithm + addProjection(FastJets(fs, FastJets::D0ILCONE, 0.7, 20.0*GeV), "ConeJets"); + } + + + /// @name Analysis methods + //@{ - void D0_1996_S3324664::init() { - - _h_deta = bookHistogram1D(1, 1, 1); - _h_dphi.addHistogram(0.0, 2.0, bookHistogram1D(2, 1, 1)); - _h_dphi.addHistogram(2.0, 4.0, bookHistogram1D(2, 1, 2)); - _h_dphi.addHistogram(4.0, 6.0, bookHistogram1D(2, 1, 3)); - _h_cosdphi_deta = bookProfile1D(3, 1, 1); - } - + void init() { + _h_deta = bookHistogram1D(1, 1, 1); + _h_dphi.addHistogram(0.0, 2.0, bookHistogram1D(2, 1, 1)); + _h_dphi.addHistogram(2.0, 4.0, bookHistogram1D(2, 1, 2)); + _h_dphi.addHistogram(4.0, 6.0, bookHistogram1D(2, 1, 3)); + _h_cosdphi_deta = bookProfile1D(3, 1, 1); + } - void D0_1996_S3324664::analyze(const Event& event) { - const double weight = event.weight(); - - Jets jets; - foreach (const Jet& jet, applyProjection<FastJets>(event, "ConeJets").jets()) { - if (fabs(jet.momentum().eta())<3.0) { - jets.push_back(jet); + void analyze(const Event& event) { + const double weight = event.weight(); + + Jets jets; + foreach (const Jet& jet, applyProjection<FastJets>(event, "ConeJets").jets()) { + if (fabs(jet.momentum().eta()) < 3.0) { + jets.push_back(jet); + } + } + + if (jets.size() < 2) { + vetoEvent; } - } - - if (jets.size()<2) { - vetoEvent; - } - FourMomentum minjet=jets[0].momentum(); - FourMomentum maxjet=jets[1].momentum(); - double mineta = minjet.eta(); - double maxeta = maxjet.eta(); - - foreach(const Jet& jet, jets) { - double eta = jet.momentum().eta(); - if (eta < mineta) { - minjet = jet.momentum(); - mineta = eta; + FourMomentum minjet = jets[0].momentum(); + FourMomentum maxjet = jets[1].momentum(); + double mineta = minjet.eta(); + double maxeta = maxjet.eta(); + + foreach(const Jet& jet, jets) { + double eta = jet.momentum().eta(); + if (eta < mineta) { + minjet = jet.momentum(); + mineta = eta; + } + else if (eta > maxeta) { + maxjet = jet.momentum(); + maxeta = eta; + } } - else if (eta > maxeta) { - maxjet = jet.momentum(); - maxeta = eta; + + if (minjet.Et()<50*GeV && maxjet.Et()<50.0*GeV) { + vetoEvent; } + + double deta = maxjet.eta()-minjet.eta(); + double dphi = mapAngle0To2Pi(maxjet.phi()-minjet.phi()); + + _h_deta->fill(deta, weight); + _h_dphi.fill(deta, 1.0-dphi/M_PI, weight); + _h_cosdphi_deta->fill(deta, cos(M_PI-dphi), weight); + } - if (minjet.Et()<50*GeV && maxjet.Et()<50.0*GeV) { - vetoEvent; - } - double deta = maxjet.eta()-minjet.eta(); - double dphi = mapAngle0To2Pi(maxjet.phi()-minjet.phi()); + void finalize() { + // Normalised to #events + normalize(_h_deta, 8830.0); + + // I have no idea what this is normalised to... in the paper it says unity! + /// @todo Understand this! + foreach (IHistogram1D* histo, _h_dphi.getHistograms()) { + normalize(histo, 0.0798); + } + + } - _h_deta->fill(deta, weight); - _h_dphi.fill(deta, 1.0-dphi/M_PI, weight); - _h_cosdphi_deta->fill(deta, cos(M_PI-dphi), weight); + //@} - } + private: - void D0_1996_S3324664::finalize() { + /// @name Histograms + //@{ - normalize(_h_deta, 8830.0); // not normalised to cross section, but to #events - - // I have no idea what this is normalised to, in the paper it says unity!?! - foreach (IHistogram1D* histo, _h_dphi.getHistograms()) { - normalize(histo, 0.0798); - } - - } + AIDA::IHistogram1D *_h_deta; + BinnedHistogram<double> _h_dphi; + AIDA::IProfile1D *_h_cosdphi_deta; + //@} + }; // This global object acts as a hook for the plugin system AnalysisBuilder<D0_1996_S3324664> plugin_D0_1996_S3324664; + } Modified: trunk/src/Analyses/D0_2001_S4674421.cc ============================================================================== --- trunk/src/Analyses/D0_2001_S4674421.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2001_S4674421.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,8 +1,8 @@ // -*- C++ -*- -#include "Rivet/Tools/Logging.hh" +#include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" +#include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleIDMethods.hh" -#include "Rivet/Analyses/D0_2001_S4674421.hh" #include "Rivet/Projections/PVertex.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" @@ -10,12 +10,21 @@ namespace Rivet { + /// @brief D0 Run I differential W/Z boson cross-section analysis + /// @author Lars Sonnenschein + class D0_2001_S4674421 : public Analysis { + public: + + /// @name Constructors etc. + //@{ + + /// Constructor. // - @c _mwmz = ratio of \f$ mW/mZ \f$ used in the publication analysis // - @c _brwenu = ratio of \f$ BR(W->e,nu) \f$ used in the publication analysis // - @c _brzee = ratio of \f$ BR(Z->ee) \f$ used in the publication analysis // - @c _mZmin = lower Z mass cut used in the publication analysis // - @c _mZmax = upper Z mass cut used in the publication analysis - D0_2001_S4674421::D0_2001_S4674421() + D0_2001_S4674421() : Analysis("D0_2001_S4674421"), _mwmz(0.8820), _brwenu(0.1073), _brzee(0.033632), _mZmin(75.*GeV), _mZmax(105.*GeV) @@ -46,107 +55,138 @@ VetoedFinalState vfs(fs); vfs.vetoNeutrinos(); addProjection(vfs, "VFS"); - + } + + + /// @name Analysis methods + //@{ + + void init() { + _eventsFilledW = 0.0; + _eventsFilledZ = 0.0; + _h_dsigdpt_w = bookHistogram1D(1, 1, 1); + _h_dsigdpt_z = bookHistogram1D(1, 1, 2); + + vector<double> bins(23); + bins += 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 50, 60, 70, 80, 100, 120, 160, 200; + _h_dsigdpt_scaled_z = bookHistogram1D("d01-x01-y03", bins); + } + void analyze(const Event& event) { + const double weight = event.weight(); - void D0_2001_S4674421::init() { - _eventsFilledW = 0.0; - _eventsFilledZ = 0.0; - _h_dsigdpt_w = bookHistogram1D(1, 1, 1); - _h_dsigdpt_z = bookHistogram1D(1, 1, 2); - - vector<double> bins(23); - bins += 0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 25, 30, 35, 40, 50, 60, 70, 80, 100, 120, 160, 200; - _h_dsigdpt_scaled_z = bookHistogram1D("d01-x01-y03", bins); - } + const LeadingParticlesFinalState& eeFS = applyProjection<LeadingParticlesFinalState>(event, "eeFS"); + if (eeFS.particles().size() == 2) { + // If there is a Z candidate: + static size_t Zcount = 0; + // Fill Z pT distributions + const ParticleVector& Zdaughters = eeFS.particles(); + const FourMomentum pmom = Zdaughters[0].momentum() + Zdaughters[1].momentum(); + double mass = sqrt(pmom.invariant()); + if (mass/GeV > _mZmin && mass/GeV < _mZmax) { + ++Zcount; + _eventsFilledZ += weight; + getLog() << Log::DEBUG << "Z #" << Zcount << " pmom.pT() = " << pmom.pT()/GeV << " GeV" << endl; + _h_dsigdpt_z->fill(pmom.pT()/GeV, weight); + _h_dsigdpt_scaled_z->fill(pmom.pT()/GeV * _mwmz, weight); + } + } else { + // There is no Z -> ee candidate... so this must be a W event, right? + const LeadingParticlesFinalState& enuFS = applyProjection<LeadingParticlesFinalState>(event, "enuFS"); + const LeadingParticlesFinalState& enubFS = applyProjection<LeadingParticlesFinalState>(event, "enubFS"); + static size_t Wcount = 0; + + // Fill W pT distributions + ParticleVector Wdaughters; + if (enuFS.particles().size() == 2 && enubFS.isEmpty()) { + Wdaughters = enuFS.particles(); + } else if (enuFS.isEmpty() && enubFS.particles().size() == 2) { + Wdaughters = enubFS.particles(); + } + if (! Wdaughters.empty()) { + assert(Wdaughters.size() == 2); + const FourMomentum pmom = Wdaughters[0].momentum() + Wdaughters[1].momentum(); + ++Wcount; + _eventsFilledW += weight; + _h_dsigdpt_w->fill(pmom.pT()/GeV, weight); + } + } + } - void D0_2001_S4674421::analyze(const Event& event) { - const double weight = event.weight(); - - const LeadingParticlesFinalState& eeFS = applyProjection<LeadingParticlesFinalState>(event, "eeFS"); - if (eeFS.particles().size() == 2) { - // If there is a Z candidate: - static size_t Zcount = 0; - // Fill Z pT distributions - const ParticleVector& Zdaughters = eeFS.particles(); - const FourMomentum pmom = Zdaughters[0].momentum() + Zdaughters[1].momentum(); - double mass = sqrt(pmom.invariant()); - if (mass/GeV > _mZmin && mass/GeV < _mZmax) { - ++Zcount; - _eventsFilledZ += weight; - getLog() << Log::DEBUG << "Z #" << Zcount << " pmom.pT() = " << pmom.pT()/GeV << " GeV" << endl; - _h_dsigdpt_z->fill(pmom.pT()/GeV, weight); - _h_dsigdpt_scaled_z->fill(pmom.pT()/GeV * _mwmz, weight); - } - } else { - // There is no Z -> ee candidate... so this must be a W event, right? - const LeadingParticlesFinalState& enuFS = applyProjection<LeadingParticlesFinalState>(event, "enuFS"); - const LeadingParticlesFinalState& enubFS = applyProjection<LeadingParticlesFinalState>(event, "enubFS"); - static size_t Wcount = 0; - - // Fill W pT distributions - ParticleVector Wdaughters; - if (enuFS.particles().size() == 2 && enubFS.isEmpty()) { - Wdaughters = enuFS.particles(); - } else if (enuFS.isEmpty() && enubFS.particles().size() == 2) { - Wdaughters = enubFS.particles(); - } - if (! Wdaughters.empty()) { - assert(Wdaughters.size() == 2); - const FourMomentum pmom = Wdaughters[0].momentum() + Wdaughters[1].momentum(); - ++Wcount; - _eventsFilledW += weight; - _h_dsigdpt_w->fill(pmom.pT()/GeV, weight); + void finalize() { + // Get cross-section per event (i.e. per unit weight) from generator + const double xSecPerEvent = crossSection()/picobarn / sumOfWeights(); + + // Correct W pT distribution to W cross-section + const double xSecW = xSecPerEvent * _eventsFilledW; + + // Correct Z pT distribution to Z cross-section + const double xSecZ = xSecPerEvent * _eventsFilledZ; + + // Get W and Z pT integrals + const double wpt_integral = integral(_h_dsigdpt_w); + const double zpt_scaled_integral = integral(_h_dsigdpt_scaled_z); + + // Divide and scale ratio histos + AIDA::IDataPointSet* div = histogramFactory().divide(histoDir() + "/d02-x01-y01", *_h_dsigdpt_w, *_h_dsigdpt_scaled_z); + div->setTitle("$[\\mathrm{d}\\sigma/\\mathrm{d}p_\\perp(W)] / [\\mathrm{d}\\sigma/\\mathrm{d}(p_\\perp(Z) \\cdot M_W/M_Z)]$"); + if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_scaled_integral == 0) { + getLog() << Log::WARN << "Not filling ratio plot because input histos are empty" << endl; + } else { + // Scale factor converts event counts to cross-sections, and inverts the + // branching ratios since only one decay channel has been analysed for each boson. + const double scalefactor = (xSecW / wpt_integral) / (xSecZ / zpt_scaled_integral) * (_brzee / _brwenu); + for (int pt = 0; pt < div->size(); ++pt) { + assert(div->point(pt)->dimension() == 2); + AIDA::IMeasurement* m = div->point(pt)->coordinate(1); + m->setValue(m->value() * scalefactor); + m->setErrorPlus(m->errorPlus() * scalefactor); + m->setErrorMinus(m->errorPlus() * scalefactor); + } } + + // Normalize non-ratio histos + normalize(_h_dsigdpt_w, xSecW); + normalize(_h_dsigdpt_z, xSecZ); + normalize(_h_dsigdpt_scaled_z, xSecZ); + } - } - + //@} + + private: + + /// Analysis used ratio of mW/mZ + const double _mwmz; + + /// Ratio of \f$ BR(W->e,nu) \f$ used in the publication analysis + const double _brwenu; + + /// Ratio of \f$ \text{BR}( Z \to e^+ e^-) \f$ used in the publication analysis + const double _brzee; + + /// Invariant mass cuts for Z boson candidate (75 GeV < mZ < 105 GeV) + const double _mZmin, _mZmax; - void D0_2001_S4674421::finalize() { - // Get cross-section per event (i.e. per unit weight) from generator - const double xSecPerEvent = crossSection()/picobarn / sumOfWeights(); - - // Correct W pT distribution to W cross-section - const double xSecW = xSecPerEvent * _eventsFilledW; - - // Correct Z pT distribution to Z cross-section - const double xSecZ = xSecPerEvent * _eventsFilledZ; - - // Get W and Z pT integrals - const double wpt_integral = integral(_h_dsigdpt_w); - const double zpt_scaled_integral = integral(_h_dsigdpt_scaled_z); - - // Divide and scale ratio histos - AIDA::IDataPointSet* div = histogramFactory().divide(histoDir() + "/d02-x01-y01", *_h_dsigdpt_w, *_h_dsigdpt_scaled_z); - div->setTitle("$[\\mathrm{d}\\sigma/\\mathrm{d}p_\\perp(W)] / [\\mathrm{d}\\sigma/\\mathrm{d}(p_\\perp(Z) \\cdot M_W/M_Z)]$"); - if (xSecW == 0 || wpt_integral == 0 || xSecZ == 0 || zpt_scaled_integral == 0) { - getLog() << Log::WARN << "Not filling ratio plot because input histos are empty" << endl; - } else { - // Scale factor converts event counts to cross-sections, and inverts the - // branching ratios since only one decay channel has been analysed for each boson. - const double scalefactor = (xSecW / wpt_integral) / (xSecZ / zpt_scaled_integral) * (_brzee / _brwenu); - for (int pt = 0; pt < div->size(); ++pt) { - assert(div->point(pt)->dimension() == 2); - AIDA::IMeasurement* m = div->point(pt)->coordinate(1); - m->setValue(m->value() * scalefactor); - m->setErrorPlus(m->errorPlus() * scalefactor); - m->setErrorMinus(m->errorPlus() * scalefactor); - } - } - // Normalize non-ratio histos - normalize(_h_dsigdpt_w, xSecW); - normalize(_h_dsigdpt_z, xSecZ); - normalize(_h_dsigdpt_scaled_z, xSecZ); + // Event counters for cross section normalizations + double _eventsFilledW; + double _eventsFilledZ; + + //@{ + /// Histograms + AIDA::IHistogram1D* _h_dsigdpt_w; + AIDA::IHistogram1D* _h_dsigdpt_z; + AIDA::IHistogram1D* _h_dsigdpt_scaled_z; + //@} - } + }; Modified: trunk/src/Analyses/D0_2004_S5992206.cc ============================================================================== --- trunk/src/Analyses/D0_2004_S5992206.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2004_S5992206.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,7 +1,7 @@ // -*- C++ -*- -#include "Rivet/Tools/Logging.hh" +#include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" -#include "Rivet/Analyses/D0_2004_S5992206.hh" +#include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/PVertex.hh" #include "Rivet/Projections/TotalVisibleMomentum.hh" @@ -9,96 +9,133 @@ namespace Rivet { - // Constructor - D0_2004_S5992206::D0_2004_S5992206() - : Analysis("D0_2004_S5992206") - { - setBeams(PROTON, ANTIPROTON); - const FinalState fs(-3.0, 3.0); - addProjection(fs, "FS"); - addProjection(FastJets(fs, FastJets::D0ILCONE, 0.7, 6*GeV), "Jets"); - addProjection(TotalVisibleMomentum(fs), "CalMET"); - addProjection(PVertex(), "PV"); - - // Veto neutrinos, and muons with pT above 1.0 GeV - VetoedFinalState vfs(fs); - vfs - .addVetoPairId(NU_E) - .addVetoPairId(NU_MU) - .addVetoPairId(NU_TAU) - .addVetoDetail(MUON, 1.0, MAXDOUBLE); - addProjection(vfs, "VFS"); - } - - - // Book histograms - void D0_2004_S5992206::init() { - _histJetAzimuth_pTmax75_100 = bookHistogram1D(1, 2, 1); - _histJetAzimuth_pTmax100_130 = bookHistogram1D(2, 2, 1); - _histJetAzimuth_pTmax130_180 = bookHistogram1D(3, 2, 1); - _histJetAzimuth_pTmax180_ = bookHistogram1D(4, 2, 1); - } - - - // Do the analysis - void D0_2004_S5992206::analyze(const Event & event) { - - // Analyse and print some info - const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); - getLog() << Log::DEBUG << "Jet multiplicity before any pT cut = " << jetpro.size() << endl; - - // Find vertex and check that its z-component is < 50 cm from the nominal IP - const PVertex& pv = applyProjection<PVertex>(event, "PV"); - if (fabs(pv.position().z())/cm > 50.0) vetoEvent; - - const Jets jets = jetpro.jetsByPt(40.0*GeV); - if (jets.size() >= 2) { - getLog() << Log::DEBUG << "Jet multiplicity after pT > 40 GeV cut = " << jets.size() << endl; - } else { - vetoEvent; + /* @brief D0 Run II jet analysis + * @author Lars Sonnenschein + * + * Measurement of angular correlations in di-jet events. + * + * + * @par Run conditions + * + * @arg \f$ \sqrt{s} = \f$ 1960 GeV + * @arg Run with generic QCD events. + * @arg Several \f$ p_\perp^\text{min} \f$ cutoffs are probably required to fill the histograms: + * @arg \f$ p_\perp^\text{min} = \f$ 50, 75, 100, 150 GeV for the four pT ranges respecively + * + */ + class D0_2004_S5992206 : public Analysis { + + public: + + /// @name Constructors etc. + //@{ + + /// Constructor. + D0_2004_S5992206() : Analysis("D0_2004_S5992206") + { + setBeams(PROTON, ANTIPROTON); + const FinalState fs(-3.0, 3.0); + addProjection(fs, "FS"); + addProjection(FastJets(fs, FastJets::D0ILCONE, 0.7, 6*GeV), "Jets"); + addProjection(TotalVisibleMomentum(fs), "CalMET"); + addProjection(PVertex(), "PV"); + + // Veto neutrinos, and muons with pT above 1.0 GeV + VetoedFinalState vfs(fs); + vfs.vetoNeutrinos(); + vfs.addVetoDetail(MUON, 1.0, MAXDOUBLE); + addProjection(vfs, "VFS"); } - const double rap1 = jets[0].momentum().rapidity(); - const double rap2 = jets[1].momentum().rapidity(); - if (fabs(rap1) > 0.5 || fabs(rap2) > 0.5) { - vetoEvent; + + //@} + + + /// @name Analysis methods + //@{ + + /// Book histograms + void init() { + _histJetAzimuth_pTmax75_100 = bookHistogram1D(1, 2, 1); + _histJetAzimuth_pTmax100_130 = bookHistogram1D(2, 2, 1); + _histJetAzimuth_pTmax130_180 = bookHistogram1D(3, 2, 1); + _histJetAzimuth_pTmax180_ = bookHistogram1D(4, 2, 1); } - getLog() << Log::DEBUG << "Jet eta and pT requirements fulfilled" << endl; - const double pT1 = jets[0].momentum().pT(); - const TotalVisibleMomentum& caloMissEt = applyProjection<TotalVisibleMomentum>(event, "CalMET"); - getLog() << Log::DEBUG << "Missing Et = " << caloMissEt.momentum().pT()/GeV << endl; - if (caloMissEt.momentum().pT() > 0.7*pT1) { - vetoEvent; + + /// Do the analysis + void analyze(const Event & event) { + + // Analyse and print some info + const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); + getLog() << Log::DEBUG << "Jet multiplicity before any pT cut = " << jetpro.size() << endl; + + // Find vertex and check that its z-component is < 50 cm from the nominal IP + const PVertex& pv = applyProjection<PVertex>(event, "PV"); + if (fabs(pv.position().z())/cm > 50.0) vetoEvent; + + const Jets jets = jetpro.jetsByPt(40.0*GeV); + if (jets.size() >= 2) { + getLog() << Log::DEBUG << "Jet multiplicity after pT > 40 GeV cut = " << jets.size() << endl; + } else { + vetoEvent; + } + const double rap1 = jets[0].momentum().rapidity(); + const double rap2 = jets[1].momentum().rapidity(); + if (fabs(rap1) > 0.5 || fabs(rap2) > 0.5) { + vetoEvent; + } + getLog() << Log::DEBUG << "Jet eta and pT requirements fulfilled" << endl; + const double pT1 = jets[0].momentum().pT(); + + const TotalVisibleMomentum& caloMissEt = applyProjection<TotalVisibleMomentum>(event, "CalMET"); + getLog() << Log::DEBUG << "Missing Et = " << caloMissEt.momentum().pT()/GeV << endl; + if (caloMissEt.momentum().pT() > 0.7*pT1) { + vetoEvent; + } + + if (pT1/GeV >= 75.0) { + const double weight = event.weight(); + const double dphi = deltaPhi(jets[0].momentum().phi(), jets[1].momentum().phi()); + if (inRange(pT1/GeV, 75.0, 100.0)) { + _histJetAzimuth_pTmax75_100->fill(dphi, weight); + } else if (inRange(pT1/GeV, 100.0, 130.0)) { + _histJetAzimuth_pTmax100_130->fill(dphi, weight); + } else if (inRange(pT1/GeV, 130.0, 180.0)) { + _histJetAzimuth_pTmax130_180->fill(dphi, weight); + } else if (pT1/GeV > 180.0) { + _histJetAzimuth_pTmax180_->fill(dphi, weight); + } + } + } - if (pT1/GeV >= 75.0) { - const double weight = event.weight(); - const double dphi = deltaPhi(jets[0].momentum().phi(), jets[1].momentum().phi()); - if (inRange(pT1/GeV, 75.0, 100.0)) { - _histJetAzimuth_pTmax75_100->fill(dphi, weight); - } else if (inRange(pT1/GeV, 100.0, 130.0)) { - _histJetAzimuth_pTmax100_130->fill(dphi, weight); - } else if (inRange(pT1/GeV, 130.0, 180.0)) { - _histJetAzimuth_pTmax130_180->fill(dphi, weight); - } else if (pT1/GeV > 180.0) { - _histJetAzimuth_pTmax180_->fill(dphi, weight); - } + + // Finalize + void finalize() { + // Normalize histograms to unit area + normalize(_histJetAzimuth_pTmax75_100); + normalize(_histJetAzimuth_pTmax100_130); + normalize(_histJetAzimuth_pTmax130_180); + normalize(_histJetAzimuth_pTmax180_); } - - } + + //@} - // Finalize - void D0_2004_S5992206::finalize() { - // Normalize histograms to unit area - normalize(_histJetAzimuth_pTmax75_100); - normalize(_histJetAzimuth_pTmax100_130); - normalize(_histJetAzimuth_pTmax130_180); - normalize(_histJetAzimuth_pTmax180_); - } + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D* _histJetAzimuth_pTmax75_100; + AIDA::IHistogram1D* _histJetAzimuth_pTmax100_130; + AIDA::IHistogram1D* _histJetAzimuth_pTmax130_180; + AIDA::IHistogram1D* _histJetAzimuth_pTmax180_; + //@} + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2004_S5992206> plugin_D0_2004_S5992206; Modified: trunk/src/Analyses/D0_2006_S6438750.cc ============================================================================== --- trunk/src/Analyses/D0_2006_S6438750.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2006_S6438750.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,5 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2006_S6438750.hh" +#include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" @@ -9,100 +9,125 @@ namespace Rivet { - D0_2006_S6438750::D0_2006_S6438750() - : Analysis("D0_2006_S6438750") - { - setBeams(PROTON, ANTIPROTON); + /// @brief Inclusive isolated photon cross-section, differential in \f$ p_\perp(gamma) \f$. + /// @author Andy Buckley + /// @author Gavin Hesketh + class D0_2006_S6438750 : public Analysis { + + public: + + /// @name Constructors etc. + //@{ + + /// Default constructor. + D0_2006_S6438750() : Analysis("D0_2006_S6438750") + { + setBeams(PROTON, ANTIPROTON); + + /// @todo Use cross-section from generator + setNeedsCrossSection(true); + + // General FS for photon isolation + FinalState fs(-1.5, 1.5); + addProjection(fs, "AllFS"); + + // Get leading photon + LeadingParticlesFinalState photonfs(fs, -1.0, 1.0); + photonfs.addParticleId(PHOTON); + addProjection(photonfs, "LeadingPhoton"); + } - /// @todo Use cross-section from generator - setNeedsCrossSection(true); + //@} - // General FS for photon isolation - FinalState fs(-1.5, 1.5); - addProjection(fs, "AllFS"); - // Get leading photon - LeadingParticlesFinalState photonfs(fs, -1.0, 1.0); - photonfs.addParticleId(PHOTON); - addProjection(photonfs, "LeadingPhoton"); - } + /// @name Analysis methods + //@{ - - - // Book histograms - void D0_2006_S6438750::init() { - _h_pTgamma = bookHistogram1D(1, 1, 1); - } - - - - // Do the analysis - void D0_2006_S6438750::analyze(const Event& event) { - - // Get the photon - const FinalState& photonfs = applyProjection<FinalState>(event, "LeadingPhoton"); - if (photonfs.particles().size() != 1) { - getLog() << Log::DEBUG << "No photon found" << endl; - vetoEvent; - } - const FourMomentum photon = photonfs.particles().front().momentum(); - if (photon.pT()/GeV < 23) { - getLog() << Log::DEBUG << "Leading photon has pT < 23 GeV: " << photon.pT()/GeV << endl; - vetoEvent; + /// Book histograms + void init() { + _h_pTgamma = bookHistogram1D(1, 1, 1); } + - // Get all other particles - const FinalState& fs = applyProjection<FinalState>(event, "AllFS"); - if (fs.isEmpty()) { - vetoEvent; - } + /// Do the analysis + void analyze(const Event& event) { - // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy - const double egamma = photon.E(); - // Energy inside R = 0.2 - double econe_02 = 0.0; - // Energy between R = [0.2, 0.4] - double econe_02_04 = 0.0; - foreach (const Particle& p, fs.particles()) { - const double dr = deltaR(photon.pseudorapidity(), photon.azimuthalAngle(), - p.momentum().pseudorapidity(), p.momentum().azimuthalAngle()); - if (dr < 0.2) { - econe_02 += p.momentum().E(); - } else if (dr < 0.4) { - econe_02_04 += p.momentum().E(); + // Get the photon + const FinalState& photonfs = applyProjection<FinalState>(event, "LeadingPhoton"); + if (photonfs.particles().size() != 1) { + getLog() << Log::DEBUG << "No photon found" << endl; + vetoEvent; } + const FourMomentum photon = photonfs.particles().front().momentum(); + if (photon.pT()/GeV < 23) { + getLog() << Log::DEBUG << "Leading photon has pT < 23 GeV: " << photon.pT()/GeV << endl; + vetoEvent; + } + + // Get all other particles + const FinalState& fs = applyProjection<FinalState>(event, "AllFS"); + if (fs.isEmpty()) { + vetoEvent; + } + + // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy + const double egamma = photon.E(); + // Energy inside R = 0.2 + double econe_02 = 0.0; + // Energy between R = [0.2, 0.4] + double econe_02_04 = 0.0; + foreach (const Particle& p, fs.particles()) { + const double dr = deltaR(photon.pseudorapidity(), photon.azimuthalAngle(), + p.momentum().pseudorapidity(), p.momentum().azimuthalAngle()); + if (dr < 0.2) { + econe_02 += p.momentum().E(); + } else if (dr < 0.4) { + econe_02_04 += p.momentum().E(); + } + } + // Veto if outer hollow cone contains more than 10% of the energy of the inner cone + // or if the non-photon energy in the inner cone exceeds 5% of the photon energy. + if (econe_02_04/econe_02 > 0.1 || (econe_02-egamma)/egamma > 0.05) { + getLog() << Log::DEBUG << "Vetoing event because photon is insufficiently isolated" << endl; + vetoEvent; + } + + // Veto if leading jet is outside plotted rapidity regions + const double eta_gamma = fabs(photon.pseudorapidity()); + if (eta_gamma > 0.9) { + getLog() << Log::DEBUG << "Leading photon falls outside acceptance range; " + << "|eta_gamma| = " << eta_gamma << endl; + vetoEvent; + } + + // Fill histo + const double weight = event.weight(); + _h_pTgamma->fill(photon.pT(), weight); } - // Veto if outer hollow cone contains more than 10% of the energy of the inner cone - // or if the non-photon energy in the inner cone exceeds 5% of the photon energy. - if (econe_02_04/econe_02 > 0.1 || (econe_02-egamma)/egamma > 0.05) { - getLog() << Log::DEBUG << "Vetoing event because photon is insufficiently isolated" << endl; - vetoEvent; - } + + - // Veto if leading jet is outside plotted rapidity regions - const double eta_gamma = fabs(photon.pseudorapidity()); - if (eta_gamma > 0.9) { - getLog() << Log::DEBUG << "Leading photon falls outside acceptance range; " - << "|eta_gamma| = " << eta_gamma << endl; - vetoEvent; + // Finalize + void finalize() { + /// @todo Generator cross-section from Pythia gives ~7500, vs. expected 2988! + //normalize(_h_pTgamma, 2988.4869); + + const double lumi_gen = sumOfWeights()/crossSection(); + // Divide by effective lumi, plus rapidity bin width of 1.8 + scale(_h_pTgamma, 1/lumi_gen * 1/1.8); } - // Fill histo - const double weight = event.weight(); - _h_pTgamma->fill(photon.pT(), weight); - } + //@} + + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D* _h_pTgamma; + //@} - // Finalize - void D0_2006_S6438750::finalize() { - /// @todo Generator cross-section from Pythia gives ~7500, vs. expected 2988! - //normalize(_h_pTgamma, 2988.4869); - - const double lumi_gen = sumOfWeights()/crossSection(); - // Divide by effective lumi, plus rapidity bin width of 1.8 - scale(_h_pTgamma, 1/lumi_gen * 1/1.8); - } + }; Modified: trunk/src/Analyses/D0_2007_S7075677.cc ============================================================================== --- trunk/src/Analyses/D0_2007_S7075677.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2007_S7075677.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,64 +1,84 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2007_S7075677.hh" +#include "Rivet/Analysis.hh" +#include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/ZFinder.hh" -#include "Rivet/RivetAIDA.hh" namespace Rivet { - D0_2007_S7075677::D0_2007_S7075677() - : Analysis("D0_2007_S7075677") - { - // Run II Z rapidity - setBeams(PROTON, ANTIPROTON); + /// @brief Measurement of D0 Run II Z pT diff cross-section shape + /// @author Andy Buckley + /// @author Gavin Hesketh + /// @author Frank Siegert + class D0_2007_S7075677 : public Analysis { + + public: + + /// Default constructor. + D0_2007_S7075677() : Analysis("D0_2007_S7075677") + { + // Run II Z rapidity + setBeams(PROTON, ANTIPROTON); + + + /// @todo Ask Gavin Hesketh about his first implemention without eta cuts. + vector<pair<double, double> > etaRanges; + // Remove eta cuts for the moment, because it seems like they have been + // corrected for. + // etaRanges.push_back(make_pair(-3.2, -1.5)); + // etaRanges.push_back(make_pair(-0.9, 0.9)); + // etaRanges.push_back(make_pair(1.5, 3.2)); + ZFinder zfinder(etaRanges, 15.0*GeV, ELECTRON, 71.0*GeV, 111.0*GeV, 0.2); + addProjection(zfinder, "ZFinder"); + } - std::vector<std::pair<double, double> > etaRanges; - // remove eta cuts for the moment, because it seems, like they have been - // corrected for. - // todo: ask gavin hesketh about it, he first implemented this - // analysis without eta cuts. - // etaRanges.push_back(make_pair(-3.2, -1.5)); - // etaRanges.push_back(make_pair(-0.9, 0.9)); - // etaRanges.push_back(make_pair(1.5, 3.2)); - ZFinder zfinder(etaRanges, 15.0*GeV, ELECTRON, 71.0*GeV, 111.0*GeV, 0.2); - addProjection(zfinder, "ZFinder"); - } + /// @name Analysis methods + //@{ - - // Book histograms - void D0_2007_S7075677::init() { - _h_yZ = bookHistogram1D(1, 1, 1); - } - + /// Book histograms + void init() { + _h_yZ = bookHistogram1D(1, 1, 1); + } - // Do the analysis - void D0_2007_S7075677::analyze(const Event & e) { - double weight = e.weight(); - - const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size() == 1) { - const ParticleVector& el(zfinder.constituentsFinalState().particles()); - if (el[0].momentum().pT() > 25.0*GeV || el[1].momentum().pT() > 25.0*GeV) { - double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); - _h_yZ->fill(yZ, weight); + /// Do the analysis + void analyze(const Event & e) { + const double weight = e.weight(); + + const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); + if (zfinder.particles().size() == 1) { + const ParticleVector& el(zfinder.constituentsFinalState().particles()); + if (el[0].momentum().pT() > 25.0*GeV || el[1].momentum().pT() > 25.0*GeV) { + double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); + _h_yZ->fill(yZ, weight); + } + } + else { + getLog() << Log::DEBUG << "No unique lepton pair found." << endl; } } - else { - getLog() << Log::DEBUG << "no unique lepton pair found." << endl; + + + // Finalize + void finalize() { + // Data seems to have been normalized for the avg of the two sides + // (+ve & -ve rapidity) rather than the sum, hence the 0.5: + normalize(_h_yZ, 0.5); } - } + //@} + + + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D * _h_yZ; + //@} - // Finalize - void D0_2007_S7075677::finalize() { - // Data seems to have been normalized for the avg of the two sides - // (+ve & -ve rapidity) rather than the sum, hence the 0.5: - normalize(_h_yZ, 0.5); - } + }; Modified: trunk/src/Analyses/D0_2008_S6879055.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S6879055.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2008_S6879055.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,6 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2008_S6879055.hh" +#include "Rivet/Analysis.hh" +#include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" @@ -7,148 +8,168 @@ #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/PVertex.hh" #include "Rivet/Projections/FastJets.hh" -#include "Rivet/RivetAIDA.hh" namespace Rivet { - D0_2008_S6879055::D0_2008_S6879055() - : Analysis("D0_2008_S6879055") - { - setBeams(PROTON, ANTIPROTON); - - // Basic final state - FinalState fs(-5.0, 5.0); - addProjection(fs, "FS"); - - // Leading electrons in tracking acceptance - LeadingParticlesFinalState lpfs(fs, -1.1, 1.1, 25*GeV); - lpfs.addParticleId(ELECTRON).addParticleId(POSITRON); - addProjection(lpfs, "LeadingElectronsFS"); - - // Invariant mass selection around Z pole - InvMassFinalState electronsFromZ(lpfs, make_pair(ELECTRON, POSITRON), 75*GeV, 105*GeV); - addProjection(electronsFromZ,"ElectronsFromZ"); - - // Vetoed FS for jets - VetoedFinalState vfs(fs); - // Add particle/antiparticle vetoing - vfs.vetoNeutrinos(); - // Veto the electrons from Z decay - vfs.addVetoOnThisFinalState(electronsFromZ); - addProjection(vfs, "JetFS"); - - // Jet finder - FastJets jets(vfs, FastJets::D0ILCONE, 0.5, 20.0*GeV); - addProjection(jets, "Jets"); - - // Vertex - PVertex vertex; - addProjection(vertex, "PrimaryVertex"); - } - - - - // Book histograms - void D0_2008_S6879055::init() { - _crossSectionRatio = bookHistogram1D(1, 1, 1); - _pTjet1 = bookHistogram1D(2, 1, 1); - _pTjet2 = bookHistogram1D(3, 1, 1); - _pTjet3 = bookHistogram1D(4, 1, 1); - } - - - - // Do the analysis - void D0_2008_S6879055::analyze(const Event& event) { - const double weight = event.weight(); - - // Skip if the event is empty - const FinalState& fs = applyProjection<FinalState>(event, "FS"); - if (fs.isEmpty()) { - vetoEvent; - } - - // Check that the primary vertex is within 60 cm in z from (0,0,0) - const PVertex& vertex = applyProjection<PVertex>(event, "PrimaryVertex"); - getLog() << Log::DEBUG << "Primary vertex is at " << vertex.position()/cm << " cm" << endl; - if (fabs(vertex.position().z())/cm > 60) { - getLog() << Log::DEBUG << "Vertex z-position " << vertex.position().z()/cm << " is outside cuts" << endl; - vetoEvent; - } - - // Find the Z candidates - const InvMassFinalState& invmassfs = applyProjection<InvMassFinalState>(event, "ElectronsFromZ"); - // If there is no Z candidate in the FinalState, skip the event - if (invmassfs.particles().size() != 2) { - getLog() << Log::DEBUG << "No Z candidate found" << endl; - vetoEvent; - } - - // Now build the list of jets on a FS without the electrons from Z - // Additional cuts on jets: |eta| < 2.5 and dR(j,leading electron) > 0.4 - const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); - const Jets jets = jetpro.jetsByPt(20.0*GeV); - vector<FourMomentum> finaljet_list; - foreach (const Jet& j, jets) { - const double jeta = j.momentum().pseudorapidity(); - const double jphi = j.momentum().azimuthalAngle(); - if (fabs(jeta) > 2.5) continue; - - FourMomentum e0 = invmassfs.particles()[0].momentum(); - FourMomentum e1 = invmassfs.particles()[1].momentum(); - const double e0eta = e0.pseudorapidity(); - const double e0phi = e0.azimuthalAngle(); - if (deltaR(e0eta, e0phi, jeta, jphi) < 0.4) continue; - - const double e1eta = e1.pseudorapidity(); - const double e1phi = e1.azimuthalAngle(); - if (deltaR(e1eta, e1phi, jeta, jphi) < 0.4) continue; - - // If we pass all cuts... - finaljet_list.push_back(j.momentum()); - } - getLog() << Log::DEBUG << "Num jets passing = " << finaljet_list.size() << endl; - - // For normalisation of crossSection data (includes events with no jets passing cuts) - _crossSectionRatio->fill(0, weight); - - // Fill jet pT and multiplicities - if (finaljet_list.size() >= 1) { - _crossSectionRatio->fill(1, weight); - _pTjet1->fill(finaljet_list[0].pT(), weight); - } - if (finaljet_list.size() >= 2) { - _crossSectionRatio->fill(2, weight); - _pTjet2->fill(finaljet_list[1].pT(), weight); - } - if (finaljet_list.size() >= 3) { - _crossSectionRatio->fill(3, weight); - _pTjet3->fill(finaljet_list[2].pT(), weight); - } - if (finaljet_list.size() >= 4) { - _crossSectionRatio->fill(4, weight); - } - } - - - - // Finalize - void D0_2008_S6879055::finalize() { - // Now divide by the inclusive result - _crossSectionRatio->scale(1.0/_crossSectionRatio->binHeight(0)); - - // Normalise jet pT's to integral of data - // there is no other way to do this, because these quantities are not - // detector corrected - normalize(_pTjet1, 10439.0); - normalize(_pTjet2, 1461.5); - normalize(_pTjet3, 217.0); - } - + /// @brief Measurement of the ratio sigma(Z/gamma* + n jets)/sigma(Z/gamma*) + class D0_2008_S6879055 : public Analysis { + public: + + /// Default constructor. + D0_2008_S6879055() : Analysis("D0_2008_S6879055") + { + setBeams(PROTON, ANTIPROTON); + + // Basic final state + FinalState fs(-5.0, 5.0); + addProjection(fs, "FS"); + + // Leading electrons in tracking acceptance + LeadingParticlesFinalState lpfs(fs, -1.1, 1.1, 25*GeV); + lpfs.addParticleId(ELECTRON).addParticleId(POSITRON); + addProjection(lpfs, "LeadingElectronsFS"); + + // Invariant mass selection around Z pole + InvMassFinalState electronsFromZ(lpfs, make_pair(ELECTRON, POSITRON), 75*GeV, 105*GeV); + addProjection(electronsFromZ,"ElectronsFromZ"); + + // Vetoed FS for jets + VetoedFinalState vfs(fs); + // Add particle/antiparticle vetoing + vfs.vetoNeutrinos(); + // Veto the electrons from Z decay + vfs.addVetoOnThisFinalState(electronsFromZ); + addProjection(vfs, "JetFS"); + + // Jet finder + FastJets jets(vfs, FastJets::D0ILCONE, 0.5, 20.0*GeV); + addProjection(jets, "Jets"); + + // Vertex + PVertex vertex; + addProjection(vertex, "PrimaryVertex"); + } + + + /// @name Analysis methods + //@{ + + // Book histograms + void init() { + _crossSectionRatio = bookHistogram1D(1, 1, 1); + _pTjet1 = bookHistogram1D(2, 1, 1); + _pTjet2 = bookHistogram1D(3, 1, 1); + _pTjet3 = bookHistogram1D(4, 1, 1); + } + + + + /// Do the analysis + void analyze(const Event& event) { + const double weight = event.weight(); + + // Skip if the event is empty + const FinalState& fs = applyProjection<FinalState>(event, "FS"); + if (fs.isEmpty()) { + vetoEvent; + } + + // Check that the primary vertex is within 60 cm in z from (0,0,0) + const PVertex& vertex = applyProjection<PVertex>(event, "PrimaryVertex"); + getLog() << Log::DEBUG << "Primary vertex is at " << vertex.position()/cm << " cm" << endl; + if (fabs(vertex.position().z())/cm > 60) { + getLog() << Log::DEBUG << "Vertex z-position " << vertex.position().z()/cm << " is outside cuts" << endl; + vetoEvent; + } + + // Find the Z candidates + const InvMassFinalState& invmassfs = applyProjection<InvMassFinalState>(event, "ElectronsFromZ"); + // If there is no Z candidate in the FinalState, skip the event + if (invmassfs.particles().size() != 2) { + getLog() << Log::DEBUG << "No Z candidate found" << endl; + vetoEvent; + } + + // Now build the list of jets on a FS without the electrons from Z + // Additional cuts on jets: |eta| < 2.5 and dR(j,leading electron) > 0.4 + const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); + const Jets jets = jetpro.jetsByPt(20.0*GeV); + vector<FourMomentum> finaljet_list; + foreach (const Jet& j, jets) { + const double jeta = j.momentum().pseudorapidity(); + const double jphi = j.momentum().azimuthalAngle(); + if (fabs(jeta) > 2.5) continue; + + FourMomentum e0 = invmassfs.particles()[0].momentum(); + FourMomentum e1 = invmassfs.particles()[1].momentum(); + const double e0eta = e0.pseudorapidity(); + const double e0phi = e0.azimuthalAngle(); + if (deltaR(e0eta, e0phi, jeta, jphi) < 0.4) continue; + + const double e1eta = e1.pseudorapidity(); + const double e1phi = e1.azimuthalAngle(); + if (deltaR(e1eta, e1phi, jeta, jphi) < 0.4) continue; + + // If we pass all cuts... + finaljet_list.push_back(j.momentum()); + } + getLog() << Log::DEBUG << "Num jets passing = " << finaljet_list.size() << endl; + + // For normalisation of crossSection data (includes events with no jets passing cuts) + _crossSectionRatio->fill(0, weight); + + // Fill jet pT and multiplicities + if (finaljet_list.size() >= 1) { + _crossSectionRatio->fill(1, weight); + _pTjet1->fill(finaljet_list[0].pT(), weight); + } + if (finaljet_list.size() >= 2) { + _crossSectionRatio->fill(2, weight); + _pTjet2->fill(finaljet_list[1].pT(), weight); + } + if (finaljet_list.size() >= 3) { + _crossSectionRatio->fill(3, weight); + _pTjet3->fill(finaljet_list[2].pT(), weight); + } + if (finaljet_list.size() >= 4) { + _crossSectionRatio->fill(4, weight); + } + } + + + + /// Finalize + void finalize() { + // Now divide by the inclusive result + _crossSectionRatio->scale(1.0/_crossSectionRatio->binHeight(0)); + + // Normalise jet pT's to integral of data + // there is no other way to do this, because these quantities are not + // detector corrected + normalize(_pTjet1, 10439.0); + normalize(_pTjet2, 1461.5); + normalize(_pTjet3, 217.0); + } + + //@} + + + private: + + /// @name Histograms + //@{ + AIDA::IHistogram1D * _crossSectionRatio; + AIDA::IHistogram1D * _pTjet1; + AIDA::IHistogram1D * _pTjet2; + AIDA::IHistogram1D * _pTjet3; + //@} + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2008_S6879055> plugin_D0_2008_S6879055; - + } Modified: trunk/src/Analyses/D0_2008_S7554427.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7554427.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2008_S7554427.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,63 +1,87 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2008_S7554427.hh" +#include "Rivet/Analysis.hh" +#include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" -#include "Rivet/RivetAIDA.hh" namespace Rivet { - D0_2008_S7554427::D0_2008_S7554427() - : Analysis("D0_2008_S7554427") - { - // Run II Z pT - setBeams(PROTON, ANTIPROTON); - - ZFinder zfinder(-MaxRapidity, MaxRapidity, 0.0*GeV, ELECTRON, - 40.0*GeV, 200.0*GeV, 0.2); - addProjection(zfinder, "ZFinder"); - } - - + /// @brief Measurement of D0 Run II Z pT differential cross-section shape + /// @author Andy Buckley + /// @author Gavin Hesketh + /// @author Frank Siegert + class D0_2008_S7554427 : public Analysis { + + public: + + /// Default constructor. + D0_2008_S7554427() + : Analysis("D0_2008_S7554427") + { + // Run II Z pT + setBeams(PROTON, ANTIPROTON); + + ZFinder zfinder(-MaxRapidity, MaxRapidity, 0.0*GeV, ELECTRON, + 40.0*GeV, 200.0*GeV, 0.2); + addProjection(zfinder, "ZFinder"); + } + + + /// @name Analysis methods + //@{ - // Book histograms - void D0_2008_S7554427::init() { - _h_ZpT = bookHistogram1D(1, 1, 1); - _h_forward_ZpT = bookHistogram1D(3, 1, 1); - } + /// Book histograms + void init() { + _h_ZpT = bookHistogram1D(1, 1, 1); + _h_forward_ZpT = bookHistogram1D(3, 1, 1); + } - // Do the analysis - void D0_2008_S7554427::analyze(const Event & e) { - double weight = e.weight(); - - const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size() == 1) { - double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); - double pTZ = zfinder.particles()[0].momentum().pT(); - _h_ZpT->fill(pTZ, weight); - if (yZ > 2.0) { - _h_forward_ZpT->fill(pTZ, weight); + /// Do the analysis + void analyze(const Event & e) { + const double weight = e.weight(); + + const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); + if (zfinder.particles().size() == 1) { + double yZ = fabs(zfinder.particles()[0].momentum().rapidity()); + double pTZ = zfinder.particles()[0].momentum().pT(); + _h_ZpT->fill(pTZ, weight); + if (yZ > 2.0) { + _h_forward_ZpT->fill(pTZ, weight); + } + } + else { + getLog() << Log::DEBUG << "no unique lepton pair found." << endl; } + } - else { - getLog() << Log::DEBUG << "no unique lepton pair found." << endl; + + + + // Finalize + void finalize() { + normalize(_h_ZpT); + normalize(_h_forward_ZpT); } - - } - + + //@} - // Finalize - void D0_2008_S7554427::finalize() { - normalize(_h_ZpT); - normalize(_h_forward_ZpT); - } + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D * _h_ZpT; + AIDA::IHistogram1D * _h_forward_ZpT; + //@} + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2008_S7554427> plugin_D0_2008_S7554427; Modified: trunk/src/Analyses/D0_2008_S7662670.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7662670.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2008_S7662670.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,105 +1,136 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2008_S7662670.hh" +#include "Rivet/Analysis.hh" +#include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "Rivet/Projections/FastJets.hh" -#include "Rivet/RivetAIDA.hh" namespace Rivet { - D0_2008_S7662670::D0_2008_S7662670() - : Analysis("D0_2008_S7662670") - { - setBeams(PROTON, ANTIPROTON); - setNeedsCrossSection(true); - - //full final state - FinalState fs(-5.0, 5.0); - addProjection(fs, "FS"); - - FastJets jetpro(fs, FastJets::D0ILCONE, 0.7, 6*GeV); - addProjection(jetpro, "Jets"); - } - - - - // Book histograms - void D0_2008_S7662670::init() - { - _h_dsigdptdy_y00_04 = bookHistogram1D(1, 1, 1); - _h_dsigdptdy_y04_08 = bookHistogram1D(2, 1, 1); - _h_dsigdptdy_y08_12 = bookHistogram1D(3, 1, 1); - _h_dsigdptdy_y12_16 = bookHistogram1D(4, 1, 1); - _h_dsigdptdy_y16_20 = bookHistogram1D(5, 1, 1); - _h_dsigdptdy_y20_24 = bookHistogram1D(6, 1, 1); - } - + /// @brief Measurement of D0 differential jet cross sections + /// @author Andy Buckley + /// @author Gavin Hesketh + class D0_2008_S7662670 : public Analysis { + + public: + + /// @name Constructors etc. + //@{ + + /// Constructor + D0_2008_S7662670() + : Analysis("D0_2008_S7662670") + { + setBeams(PROTON, ANTIPROTON); + setNeedsCrossSection(true); + + // Full final state + FinalState fs(-5.0, 5.0); + addProjection(fs, "FS"); + + FastJets jetpro(fs, FastJets::D0ILCONE, 0.7, 6*GeV); + addProjection(jetpro, "Jets"); + } + + //@} - // Do the analysis - void D0_2008_S7662670::analyze(const Event & event) { - const double weight = event.weight(); - - // Skip if the event is empty - const FinalState& fs = applyProjection<FinalState>(event, "FS"); - if (fs.isEmpty()) { - getLog() << Log::DEBUG << "Empty event!" << endl; - vetoEvent; - } + /// @name Analysis methods + //@{ - // Find the jets - const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); - // If there are no jets, skip the event - if (jetpro.jets().size() == 0) { - getLog() << Log::DEBUG << "No jets found" << endl; - vetoEvent; + /// Book histograms + void init() + { + _h_dsigdptdy_y00_04 = bookHistogram1D(1, 1, 1); + _h_dsigdptdy_y04_08 = bookHistogram1D(2, 1, 1); + _h_dsigdptdy_y08_12 = bookHistogram1D(3, 1, 1); + _h_dsigdptdy_y12_16 = bookHistogram1D(4, 1, 1); + _h_dsigdptdy_y16_20 = bookHistogram1D(5, 1, 1); + _h_dsigdptdy_y20_24 = bookHistogram1D(6, 1, 1); } + + + + /// Do the analysis + void analyze(const Event& event) { + const double weight = event.weight(); + + // Skip if the event is empty + const FinalState& fs = applyProjection<FinalState>(event, "FS"); + if (fs.isEmpty()) { + getLog() << Log::DEBUG << "Empty event!" << endl; + vetoEvent; + } + + // Find the jets + const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets"); + // If there are no jets, skip the event + if (jetpro.jets().size() == 0) { + getLog() << Log::DEBUG << "No jets found" << endl; + vetoEvent; + } - // Fill histo for each jet - foreach (const Jet& j, jetpro.jets()) { - const double pt = j.momentum().pT(); - const double y = fabs(j.momentum().rapidity()); - if (pt/GeV > 50) { - getLog() << Log::TRACE << "Filling histos: pT = " << pt/GeV - << ", |y| = " << y << endl; - if (y < 0.4) { - _h_dsigdptdy_y00_04->fill(pt/GeV, weight); - } else if (y < 0.8) { - _h_dsigdptdy_y04_08->fill(pt/GeV, weight); - } else if (y < 1.2) { - _h_dsigdptdy_y08_12->fill(pt/GeV, weight); - } else if (y < 1.6) { - _h_dsigdptdy_y12_16->fill(pt/GeV, weight); - } else if (y < 2.0) { - _h_dsigdptdy_y16_20->fill(pt/GeV, weight); - } else if (y < 2.4) { - _h_dsigdptdy_y20_24->fill(pt/GeV, weight); + // Fill histo for each jet + foreach (const Jet& j, jetpro.jets()) { + const double pt = j.momentum().pT(); + const double y = fabs(j.momentum().rapidity()); + if (pt/GeV > 50) { + getLog() << Log::TRACE << "Filling histos: pT = " << pt/GeV + << ", |y| = " << y << endl; + if (y < 0.4) { + _h_dsigdptdy_y00_04->fill(pt/GeV, weight); + } else if (y < 0.8) { + _h_dsigdptdy_y04_08->fill(pt/GeV, weight); + } else if (y < 1.2) { + _h_dsigdptdy_y08_12->fill(pt/GeV, weight); + } else if (y < 1.6) { + _h_dsigdptdy_y12_16->fill(pt/GeV, weight); + } else if (y < 2.0) { + _h_dsigdptdy_y16_20->fill(pt/GeV, weight); + } else if (y < 2.4) { + _h_dsigdptdy_y20_24->fill(pt/GeV, weight); + } } } + } + - } - - - // Finalize - void D0_2008_S7662670::finalize() { - /// Scale by L_eff = sig_MC * L_exp / num_MC - const double lumi_mc = sumOfWeights() / crossSection(); - const double scalefactor = 1 / lumi_mc; - scale(_h_dsigdptdy_y00_04, scalefactor); - scale(_h_dsigdptdy_y04_08, scalefactor); - scale(_h_dsigdptdy_y08_12, scalefactor); - scale(_h_dsigdptdy_y12_16, scalefactor); - scale(_h_dsigdptdy_y16_20, scalefactor); - scale(_h_dsigdptdy_y20_24, scalefactor); - } + /// Finalize + void finalize() { + /// Scale by L_eff = sig_MC * L_exp / num_MC + const double lumi_mc = sumOfWeights() / crossSection(); + const double scalefactor = 1 / lumi_mc; + scale(_h_dsigdptdy_y00_04, scalefactor); + scale(_h_dsigdptdy_y04_08, scalefactor); + scale(_h_dsigdptdy_y08_12, scalefactor); + scale(_h_dsigdptdy_y12_16, scalefactor); + scale(_h_dsigdptdy_y16_20, scalefactor); + scale(_h_dsigdptdy_y20_24, scalefactor); + } + //@} + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D* _h_dsigdptdy_y00_04; + AIDA::IHistogram1D* _h_dsigdptdy_y04_08; + AIDA::IHistogram1D* _h_dsigdptdy_y08_12; + AIDA::IHistogram1D* _h_dsigdptdy_y12_16; + AIDA::IHistogram1D* _h_dsigdptdy_y16_20; + AIDA::IHistogram1D* _h_dsigdptdy_y20_24; + //@} + + }; + + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2008_S7662670> plugin_D0_2008_S7662670; - + } Modified: trunk/src/Analyses/D0_2008_S7719523.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7719523.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2008_S7719523.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,5 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2008_S7719523.hh" +#include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/LeadingParticlesFinalState.hh" @@ -10,165 +10,195 @@ namespace Rivet { - D0_2008_S7719523::D0_2008_S7719523() - : Analysis("D0_2008_S7719523") - { - setBeams(PROTON, ANTIPROTON); - setNeedsCrossSection(true); - - // General FS - FinalState fs(-5.0, 5.0); - addProjection(fs, "FS"); - - // Get leading photon - LeadingParticlesFinalState photonfs(fs, -1.0, 1.0); - photonfs.addParticleId(PHOTON); - addProjection(photonfs, "LeadingPhoton"); - - // FS for jets excludes the leading photon - VetoedFinalState vfs(fs); - vfs.addVetoOnThisFinalState(photonfs); - addProjection(vfs, "JetFS"); - } - - - - // Book histograms - void D0_2008_S7719523::init() { - /// @todo Cross-section units in label - - _h_central_same_cross_section = bookHistogram1D(1, 1, 1); - _h_central_opp_cross_section = bookHistogram1D(2, 1, 1); - _h_forward_same_cross_section = bookHistogram1D(3, 1, 1); - _h_forward_opp_cross_section = bookHistogram1D(4, 1, 1); - } - - + /// @brief Measurement of isolated gamma + jet + X differential cross-sections + /// Inclusive isolated gamma + jet cross-sections, differential in pT(gamma), for + /// various photon and jet rapidity bins. + /// + /// @author Andy Buckley + /// @author Gavin Hesketh + class D0_2008_S7719523 : public Analysis { + + public: + + /// @name Constructors etc. + //@{ + + /// Constructor + D0_2008_S7719523() + : Analysis("D0_2008_S7719523") + { + setBeams(PROTON, ANTIPROTON); + setNeedsCrossSection(true); + + // General FS + FinalState fs(-5.0, 5.0); + addProjection(fs, "FS"); + + // Get leading photon + LeadingParticlesFinalState photonfs(fs, -1.0, 1.0); + photonfs.addParticleId(PHOTON); + addProjection(photonfs, "LeadingPhoton"); + + // FS for jets excludes the leading photon + VetoedFinalState vfs(fs); + vfs.addVetoOnThisFinalState(photonfs); + addProjection(vfs, "JetFS"); + } + + //@} - // Do the analysis - void D0_2008_S7719523::analyze(const Event& event) { - const double weight = event.weight(); - - // Get the photon - const FinalState& photonfs = applyProjection<FinalState>(event, "LeadingPhoton"); - if (photonfs.particles().size() != 1) { - getLog() << Log::DEBUG << "No photon found" << endl; - vetoEvent; - } - const FourMomentum photon = photonfs.particles().front().momentum(); - if (photon.pT()/GeV < 30) { - getLog() << Log::DEBUG << "Leading photon has pT < 30 GeV: " << photon.pT()/GeV << endl; - vetoEvent; - } - // Get all charged particles - const FinalState& fs = applyProjection<FinalState>(event, "JetFS"); - if (fs.isEmpty()) { - vetoEvent; + /// @name Analysis methods + //@{ + + /// Book histograms + void init() { + _h_central_same_cross_section = bookHistogram1D(1, 1, 1); + _h_central_opp_cross_section = bookHistogram1D(2, 1, 1); + _h_forward_same_cross_section = bookHistogram1D(3, 1, 1); + _h_forward_opp_cross_section = bookHistogram1D(4, 1, 1); } + + - // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy - const double egamma = photon.E(); - double econe = 0.0; - foreach (const Particle& p, fs.particles()) { - if (deltaR(photon, p.momentum()) < 0.4) { - econe += p.momentum().E(); - // Veto as soon as E_cone gets larger - if (econe/egamma > 0.07) { - getLog() << Log::DEBUG << "Vetoing event because photon is insufficiently isolated" << endl; - vetoEvent; + /// Do the analysis + void analyze(const Event& event) { + const double weight = event.weight(); + + // Get the photon + const FinalState& photonfs = applyProjection<FinalState>(event, "LeadingPhoton"); + if (photonfs.particles().size() != 1) { + getLog() << Log::DEBUG << "No photon found" << endl; + vetoEvent; + } + const FourMomentum photon = photonfs.particles().front().momentum(); + if (photon.pT()/GeV < 30) { + getLog() << Log::DEBUG << "Leading photon has pT < 30 GeV: " << photon.pT()/GeV << endl; + vetoEvent; + } + + // Get all charged particles + const FinalState& fs = applyProjection<FinalState>(event, "JetFS"); + if (fs.isEmpty()) { + vetoEvent; + } + + // Isolate photon by ensuring that a 0.4 cone around it contains less than 7% of the photon's energy + const double egamma = photon.E(); + double econe = 0.0; + foreach (const Particle& p, fs.particles()) { + if (deltaR(photon, p.momentum()) < 0.4) { + econe += p.momentum().E(); + // Veto as soon as E_cone gets larger + if (econe/egamma > 0.07) { + getLog() << Log::DEBUG << "Vetoing event because photon is insufficiently isolated" << endl; + vetoEvent; + } } } - } - - - /// @todo Allow proj creation w/o FS as ctor arg, so that calc can be used more easily. - FastJets jetpro(fs, FastJets::D0ILCONE, 0.7); //< @todo This fs arg makes no sense! - jetpro.calc(fs.particles()); - Jets isolated_jets; - foreach (const Jet& j, jetpro.jets()) { - const FourMomentum pjet = j.momentum(); - const double dr = deltaR(photon.pseudorapidity(), photon.azimuthalAngle(), - pjet.pseudorapidity(), pjet.azimuthalAngle()); - if (dr > 0.7 && pjet.pT()/GeV > 15) { - isolated_jets.push_back(j); + + + /// @todo Allow proj creation w/o FS as ctor arg, so that calc can be used more easily. + FastJets jetpro(fs, FastJets::D0ILCONE, 0.7); //< @todo This fs arg makes no sense! + jetpro.calc(fs.particles()); + Jets isolated_jets; + foreach (const Jet& j, jetpro.jets()) { + const FourMomentum pjet = j.momentum(); + const double dr = deltaR(photon.pseudorapidity(), photon.azimuthalAngle(), + pjet.pseudorapidity(), pjet.azimuthalAngle()); + if (dr > 0.7 && pjet.pT()/GeV > 15) { + isolated_jets.push_back(j); + } + } + + getLog() << Log::DEBUG << "Num jets after isolation and pT cuts = " + << isolated_jets.size() << endl; + if (isolated_jets.empty()) { + getLog() << Log::DEBUG << "No jets pass cuts" << endl; + vetoEvent; } + + // Sort by pT and get leading jet + sort(isolated_jets.begin(), isolated_jets.end(), cmpJetsByPt); + const FourMomentum leadingJet = isolated_jets.front().momentum(); + int photon_jet_sign = sign( leadingJet.rapidity() * photon.rapidity() ); + + // Veto if leading jet is outside plotted rapidity regions + const double abs_y1 = fabs(leadingJet.rapidity()); + if (inRange(abs_y1, 0.8, 1.5) || abs_y1 > 2.5) { + getLog() << Log::DEBUG << "Leading jet falls outside acceptance range; |y1| = " + << abs_y1 << endl; + vetoEvent; + } + + // Fill histos + if (fabs(leadingJet.rapidity()) < 0.8) { + if (photon_jet_sign >= 1) { + _h_central_same_cross_section->fill(photon.pT(), weight); + } else { + _h_central_opp_cross_section->fill(photon.pT(), weight); + } + } else if (inRange( fabs(leadingJet.rapidity()), 1.5, 2.5)) { + if (photon_jet_sign >= 1) { + _h_forward_same_cross_section->fill(photon.pT(), weight); + } else { + _h_forward_opp_cross_section->fill(photon.pT(), weight); + } + } + } - getLog() << Log::DEBUG << "Num jets after isolation and pT cuts = " - << isolated_jets.size() << endl; - if (isolated_jets.empty()) { - getLog() << Log::DEBUG << "No jets pass cuts" << endl; - vetoEvent; - } - - // Sort by pT and get leading jet - sort(isolated_jets.begin(), isolated_jets.end(), cmpJetsByPt); - const FourMomentum leadingJet = isolated_jets.front().momentum(); - int photon_jet_sign = sign( leadingJet.rapidity() * photon.rapidity() ); - - // Veto if leading jet is outside plotted rapidity regions - const double abs_y1 = fabs(leadingJet.rapidity()); - if (inRange(abs_y1, 0.8, 1.5) || abs_y1 > 2.5) { - getLog() << Log::DEBUG << "Leading jet falls outside acceptance range; |y1| = " - << abs_y1 << endl; - vetoEvent; - } - - // Fill histos - if (fabs(leadingJet.rapidity()) < 0.8) { - if (photon_jet_sign >= 1) { - _h_central_same_cross_section->fill(photon.pT(), weight); - } else { - _h_central_opp_cross_section->fill(photon.pT(), weight); - } - } else if (inRange( fabs(leadingJet.rapidity()), 1.5, 2.5)) { - if (photon_jet_sign >= 1) { - _h_forward_same_cross_section->fill(photon.pT(), weight); - } else { - _h_forward_opp_cross_section->fill(photon.pT(), weight); - } + + + /// Finalize + void finalize() { + const double lumi_gen = sumOfWeights()/crossSection(); + const double dy_photon = 2.0; + const double dy_jet_central = 1.6; + const double dy_jet_forward = 2.0; + + // Cross-section ratios (6 plots) + // Central/central and forward/forward ratios + AIDA::IHistogramFactory& hf = histogramFactory(); + const string dir = histoDir(); + + hf.divide(dir + "/d05-x01-y01", *_h_central_opp_cross_section, *_h_central_same_cross_section); + hf.divide(dir + "/d08-x01-y01", *_h_forward_opp_cross_section, *_h_forward_same_cross_section); + + // Central/forward ratio combinations + hf.divide(dir + "/d06-x01-y01", *_h_central_same_cross_section, + *_h_forward_same_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); + hf.divide(dir + "/d07-x01-y01", *_h_central_opp_cross_section, + *_h_forward_same_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); + hf.divide(dir + "/d09-x01-y01", *_h_central_same_cross_section, + *_h_forward_opp_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); + hf.divide(dir + "/d10-x01-y01", *_h_central_opp_cross_section, + *_h_forward_opp_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); + + // Use generator cross section for remaining histograms + scale(_h_central_same_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central); + scale(_h_central_opp_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central); + scale(_h_forward_same_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward); + scale(_h_forward_opp_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward); } + + //@} - } - - - - // Finalize - void D0_2008_S7719523::finalize() { - const double lumi_gen = sumOfWeights()/crossSection(); - const double dy_photon = 2.0; - const double dy_jet_central = 1.6; - const double dy_jet_forward = 2.0; - - // Cross-section ratios (6 plots) - // Central/central and forward/forward ratios - AIDA::IHistogramFactory& hf = histogramFactory(); - const string dir = histoDir(); - - hf.divide(dir + "/d05-x01-y01", *_h_central_opp_cross_section, *_h_central_same_cross_section); - hf.divide(dir + "/d08-x01-y01", *_h_forward_opp_cross_section, *_h_forward_same_cross_section); - - // Central/forward ratio combinations - hf.divide(dir + "/d06-x01-y01", *_h_central_same_cross_section, - *_h_forward_same_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); - hf.divide(dir + "/d07-x01-y01", *_h_central_opp_cross_section, - *_h_forward_same_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); - hf.divide(dir + "/d09-x01-y01", *_h_central_same_cross_section, - *_h_forward_opp_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); - hf.divide(dir + "/d10-x01-y01", *_h_central_opp_cross_section, - *_h_forward_opp_cross_section)->scale(dy_jet_forward/dy_jet_central, 1); - - // Use generator cross section for remaining histograms - scale(_h_central_same_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central); - scale(_h_central_opp_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_central); - scale(_h_forward_same_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward); - scale(_h_forward_opp_cross_section, 1.0/lumi_gen * 1.0/dy_photon * 1.0/dy_jet_forward); - } + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D* _h_central_same_cross_section; + AIDA::IHistogram1D* _h_central_opp_cross_section; + AIDA::IHistogram1D* _h_forward_same_cross_section; + AIDA::IHistogram1D* _h_forward_opp_cross_section; + //@} + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2008_S7719523> plugin_D0_2008_S7719523; - + } Modified: trunk/src/Analyses/D0_2008_S7837160.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7837160.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2008_S7837160.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,5 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2008_S7837160.hh" +#include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Tools/ParticleIDMethods.hh" #include "Rivet/Projections/FinalState.hh" @@ -10,168 +10,193 @@ namespace Rivet { - D0_2008_S7837160::D0_2008_S7837160() - : Analysis("D0_2008_S7837160") - { - // Run II W charge asymmetry - setBeams(PROTON, ANTIPROTON); - - // Leading electrons - FinalState fs(-5.0, 5.0); - - LeadingParticlesFinalState efs(fs); - efs.addParticleId(ELECTRON).addParticleId(POSITRON); - addProjection(efs, "WDecayE"); - - LeadingParticlesFinalState nufs(fs); - nufs.addParticleId(NU_E).addParticleId(NU_EBAR); - addProjection(nufs, "WDecayNu"); - - // Final state w/o electron - VetoedFinalState vfs(fs); - /// @todo A better way would be to have a "only photons FS". Add this projection. - /// @todo Make it easier to exclude all neutrinos - vfs.addVetoOnThisFinalState(efs); - vfs.addVetoPairId(NU_E).addVetoPairId(NU_MU).addVetoPairId(NU_TAU); - addProjection(vfs, "NoElectronFS"); - } - - - - // Book histograms - void D0_2008_S7837160::init() { - _h_dsigplus_deta_25_35 = bookHistogram1D("dsigplus_deta_25_35", 10, 0.0, 3.2); - _h_dsigminus_deta_25_35 = bookHistogram1D("dsigminus_deta_25_35", 10, 0.0, 3.2); - _h_dsigplus_deta_35 = bookHistogram1D("dsigplus_deta_35", 10, 0.0, 3.2); - _h_dsigminus_deta_35 = bookHistogram1D("dsigminus_deta_35", 10, 0.0, 3.2); - _h_dsigplus_deta_25 = bookHistogram1D("dsigplus_deta_25", 10, 0.0, 3.2); - _h_dsigminus_deta_25 = bookHistogram1D("dsigminus_deta_25", 10, 0.0, 3.2); - } - - - // Do the analysis - void D0_2008_S7837160::analyze(const Event & event) { - const double weight = event.weight(); - - // Find the W decay products - const FinalState& efs = applyProjection<FinalState>(event, "WDecayE"); - const FinalState& nufs = applyProjection<FinalState>(event, "WDecayNu"); - - // If there is no e/nu_e pair in the FinalState, skip the event - if (efs.particles().size() < 1 || nufs.particles().size() < 1) { - getLog() << Log::DEBUG << "No e/nu_e pair found " << endl; - vetoEvent; - } - - // Identify leading nu and electron - ParticleVector es = efs.particles(); - sort(es.begin(), es.end(), cmpParticleByEt); - Particle leading_e = es[0]; - // - ParticleVector nus = nufs.particles(); - sort(nus.begin(), nus.end(), cmpParticleByEt); - Particle leading_nu = nus[0]; - - // Require that the neutrino has Et > 25 GeV - const FourMomentum nu = leading_nu.momentum(); - if (nu.Et()/GeV < 25) { - getLog() << Log::DEBUG << "Neutrino fails Et cut" << endl; - vetoEvent; - } - - // Get "raw" electron 4-momentum and add back in photons that could have radiated from the electron - FourMomentum e = leading_e.momentum(); - const ParticleVector allparts = applyProjection<FinalState>(event, "NoElectronFS").particles(); - const double HALO_RADIUS = 0.2; - foreach (const Particle& p, allparts) { - if (p.pdgId() == PHOTON) { - const double pho_eta = p.momentum().pseudorapidity(); - const double pho_phi = p.momentum().azimuthalAngle(); - if (deltaR(e.pseudorapidity(), e.azimuthalAngle(), pho_eta, pho_phi) < HALO_RADIUS) { - e += p.momentum(); + /// @brief Measurement of W charge asymmetry from D0 Run II + /// @author Andy Buckley + /// @author Gavin Hesketh + class D0_2008_S7837160 : public Analysis { + + public: + + /// Default constructor. + D0_2008_S7837160() + : Analysis("D0_2008_S7837160") + { + // Run II W charge asymmetry + setBeams(PROTON, ANTIPROTON); + + // Leading electrons + FinalState fs(-5.0, 5.0); + + LeadingParticlesFinalState efs(fs); + efs.addParticleId(ELECTRON).addParticleId(POSITRON); + addProjection(efs, "WDecayE"); + + LeadingParticlesFinalState nufs(fs); + nufs.addParticleId(NU_E).addParticleId(NU_EBAR); + addProjection(nufs, "WDecayNu"); + + // Final state w/o electron + VetoedFinalState vfs(fs); + /// @todo A better way would be to have a "only photons FS". Add this projection. + vfs.addVetoOnThisFinalState(efs); + vfs.vetoNeutrinos(); + addProjection(vfs, "NoElectronFS"); + } + + + /// @name Analysis methods + //@{ + + // Book histograms + void init() { + _h_dsigplus_deta_25_35 = bookHistogram1D("dsigplus_deta_25_35", 10, 0.0, 3.2); + _h_dsigminus_deta_25_35 = bookHistogram1D("dsigminus_deta_25_35", 10, 0.0, 3.2); + _h_dsigplus_deta_35 = bookHistogram1D("dsigplus_deta_35", 10, 0.0, 3.2); + _h_dsigminus_deta_35 = bookHistogram1D("dsigminus_deta_35", 10, 0.0, 3.2); + _h_dsigplus_deta_25 = bookHistogram1D("dsigplus_deta_25", 10, 0.0, 3.2); + _h_dsigminus_deta_25 = bookHistogram1D("dsigminus_deta_25", 10, 0.0, 3.2); + } + + + /// Do the analysis + void analyze(const Event & event) { + const double weight = event.weight(); + + // Find the W decay products + const FinalState& efs = applyProjection<FinalState>(event, "WDecayE"); + const FinalState& nufs = applyProjection<FinalState>(event, "WDecayNu"); + + // If there is no e/nu_e pair in the FinalState, skip the event + if (efs.particles().size() < 1 || nufs.particles().size() < 1) { + getLog() << Log::DEBUG << "No e/nu_e pair found " << endl; + vetoEvent; + } + + // Identify leading nu and electron + ParticleVector es = efs.particles(); + sort(es.begin(), es.end(), cmpParticleByEt); + Particle leading_e = es[0]; + // + ParticleVector nus = nufs.particles(); + sort(nus.begin(), nus.end(), cmpParticleByEt); + Particle leading_nu = nus[0]; + + // Require that the neutrino has Et > 25 GeV + const FourMomentum nu = leading_nu.momentum(); + if (nu.Et() < 25*GeV) { + getLog() << Log::DEBUG << "Neutrino fails Et cut" << endl; + vetoEvent; + } + + // Get "raw" electron 4-momentum and add back in photons that could have radiated from the electron + FourMomentum e = leading_e.momentum(); + /// @todo Use ClusteredPhotons photon summing projection + const ParticleVector allparts = applyProjection<FinalState>(event, "NoElectronFS").particles(); + const double HALO_RADIUS = 0.2; + foreach (const Particle& p, allparts) { + if (p.pdgId() == PHOTON) { + const double pho_eta = p.momentum().pseudorapidity(); + const double pho_phi = p.momentum().azimuthalAngle(); + if (deltaR(e.pseudorapidity(), e.azimuthalAngle(), pho_eta, pho_phi) < HALO_RADIUS) { + e += p.momentum(); + } } } - } - - // Require that the electron has Et > 25 GeV - if (e.Et()/GeV < 25) { - getLog() << Log::DEBUG << "Electron fails Et cut" << endl; - vetoEvent; - } - - - const double eta_e = fabs(e.pseudorapidity()); - const double et_e = e.Et(); - const int chg_e = PID::threeCharge(leading_e.pdgId()); - if (et_e/GeV < 35) { - // 25 < ET < 35 - if (chg_e < 0) { - _h_dsigminus_deta_25_35->fill(eta_e, weight); + + // Require that the electron has Et > 25 GeV + if (e.Et() < 25*GeV) { + getLog() << Log::DEBUG << "Electron fails Et cut" << endl; + vetoEvent; + } + + + const double eta_e = fabs(e.pseudorapidity()); + const double et_e = e.Et(); + const int chg_e = PID::threeCharge(leading_e.pdgId()); + if (et_e < 35*GeV) { + // 25 < ET < 35 + if (chg_e < 0) { + _h_dsigminus_deta_25_35->fill(eta_e, weight); + } else { + _h_dsigplus_deta_25_35->fill(eta_e, weight); + } } else { - _h_dsigplus_deta_25_35->fill(eta_e, weight); + // ET > 35 + if (chg_e < 0) { + _h_dsigminus_deta_35->fill(eta_e, weight); + } else { + _h_dsigplus_deta_35->fill(eta_e, weight); + } } - } else { - // ET > 35 + // Inclusive: ET > 25 if (chg_e < 0) { - _h_dsigminus_deta_35->fill(eta_e, weight); + _h_dsigminus_deta_25->fill(eta_e, weight); } else { - _h_dsigplus_deta_35->fill(eta_e, weight); + _h_dsigplus_deta_25->fill(eta_e, weight); } } - // Inclusive: ET > 25 - if (chg_e < 0) { - _h_dsigminus_deta_25->fill(eta_e, weight); - } else { - _h_dsigplus_deta_25->fill(eta_e, weight); - } - } - - - // Finalize - void D0_2008_S7837160::finalize() { - // Construct asymmetry: (dsig+/deta - dsig-/deta) / (dsig+/deta + dsig-/deta) for each Et region - AIDA::IHistogramFactory& hf = histogramFactory(); - - const string basetitle = "W charge asymmetry for "; - const string xlabel = "$|\\eta|$ of leading electron"; - const string ylabel = "A = " - "$(\\frac{\\mathrm{d}{\\sigma^+}}{\\mathrm{d}{|\\eta|}} - \\frac{\\mathrm{d}{\\sigma^-}}{\\mathrm{d}{|\\eta|}}) / " - "(\\frac{\\mathrm{d}{\\sigma^+}}{\\mathrm{d}{|\\eta|}} + \\frac{\\mathrm{d}{\\sigma^-}}{\\mathrm{d}{|\\eta|}})$"; - - IHistogram1D* num25_35 = hf.subtract("/num25_35", *_h_dsigplus_deta_25_35, *_h_dsigminus_deta_25_35); - IHistogram1D* denom25_35 = hf.add("/denom25_35", *_h_dsigplus_deta_25_35, *_h_dsigminus_deta_25_35); - assert(num25_35 && denom25_35); - IDataPointSet* tot25_35 = hf.divide(histoDir() + "/d01-x01-y01", *num25_35, *denom25_35); - tot25_35->setTitle(basetitle + "$25 < E_\\perp < 35$ GeV"); - tot25_35->setXTitle(xlabel); - tot25_35->setYTitle(ylabel); - hf.destroy(num25_35); - hf.destroy(denom25_35); - // - IHistogram1D* num35 = hf.subtract("/num35", *_h_dsigplus_deta_35, *_h_dsigminus_deta_35); - IHistogram1D* denom35 = hf.add("/denom35", *_h_dsigplus_deta_35, *_h_dsigminus_deta_35); - assert(num35 && denom35); - IDataPointSet* tot35 = hf.divide(histoDir() + "/d02-x01-y01", *num35, *denom35); - tot35->setTitle(basetitle + "$E_\\perp > 35$ GeV"); - tot35->setXTitle(xlabel); - tot35->setYTitle(ylabel); - hf.destroy(num35); - hf.destroy(denom35); - // - IHistogram1D* num25 = hf.subtract("/num25", *_h_dsigplus_deta_25, *_h_dsigminus_deta_25); - IHistogram1D* denom25 = hf.add("/denom25", *_h_dsigplus_deta_25, *_h_dsigminus_deta_25); - assert(num25 && denom25); - IDataPointSet* tot25 = hf.divide(histoDir() + "/d03-x01-y01", *num25, *denom25); - tot25->setTitle(basetitle + "$E_\\perp > 35$ GeV"); - tot25->setXTitle(xlabel); - tot25->setYTitle(ylabel); - hf.destroy(num25); - hf.destroy(denom25); - } - - - + + + /// Finalize + void finalize() { + // Construct asymmetry: (dsig+/deta - dsig-/deta) / (dsig+/deta + dsig-/deta) for each Et region + AIDA::IHistogramFactory& hf = histogramFactory(); + + const string basetitle = "W charge asymmetry for "; + const string xlabel = "$|\\eta|$ of leading electron"; + const string ylabel = "A = " + "$(\\frac{\\mathrm{d}{\\sigma^+}}{\\mathrm{d}{|\\eta|}} - \\frac{\\mathrm{d}{\\sigma^-}}{\\mathrm{d}{|\\eta|}}) / " + "(\\frac{\\mathrm{d}{\\sigma^+}}{\\mathrm{d}{|\\eta|}} + \\frac{\\mathrm{d}{\\sigma^-}}{\\mathrm{d}{|\\eta|}})$"; + + IHistogram1D* num25_35 = hf.subtract("/num25_35", *_h_dsigplus_deta_25_35, *_h_dsigminus_deta_25_35); + IHistogram1D* denom25_35 = hf.add("/denom25_35", *_h_dsigplus_deta_25_35, *_h_dsigminus_deta_25_35); + assert(num25_35 && denom25_35); + IDataPointSet* tot25_35 = hf.divide(histoDir() + "/d01-x01-y01", *num25_35, *denom25_35); + tot25_35->setTitle(basetitle + "$25 < E_\\perp < 35$ GeV"); + tot25_35->setXTitle(xlabel); + tot25_35->setYTitle(ylabel); + hf.destroy(num25_35); + hf.destroy(denom25_35); + // + IHistogram1D* num35 = hf.subtract("/num35", *_h_dsigplus_deta_35, *_h_dsigminus_deta_35); + IHistogram1D* denom35 = hf.add("/denom35", *_h_dsigplus_deta_35, *_h_dsigminus_deta_35); + assert(num35 && denom35); + IDataPointSet* tot35 = hf.divide(histoDir() + "/d02-x01-y01", *num35, *denom35); + tot35->setTitle(basetitle + "$E_\\perp > 35$ GeV"); + tot35->setXTitle(xlabel); + tot35->setYTitle(ylabel); + hf.destroy(num35); + hf.destroy(denom35); + // + IHistogram1D* num25 = hf.subtract("/num25", *_h_dsigplus_deta_25, *_h_dsigminus_deta_25); + IHistogram1D* denom25 = hf.add("/denom25", *_h_dsigplus_deta_25, *_h_dsigminus_deta_25); + assert(num25 && denom25); + IDataPointSet* tot25 = hf.divide(histoDir() + "/d03-x01-y01", *num25, *denom25); + tot25->setTitle(basetitle + "$E_\\perp > 35$ GeV"); + tot25->setXTitle(xlabel); + tot25->setYTitle(ylabel); + hf.destroy(num25); + hf.destroy(denom25); + } + + //@} + + + private: + + /// @name Histograms + //@{ + /// @todo Move into histo manager + AIDA::IHistogram1D *_h_dsigplus_deta_25_35, *_h_dsigminus_deta_25_35; + AIDA::IHistogram1D *_h_dsigplus_deta_35, *_h_dsigminus_deta_35; + AIDA::IHistogram1D *_h_dsigplus_deta_25, *_h_dsigminus_deta_25; + //@} + + }; + + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2008_S7837160> plugin_D0_2008_S7837160; - + } Modified: trunk/src/Analyses/D0_2008_S7863608.cc ============================================================================== --- trunk/src/Analyses/D0_2008_S7863608.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2008_S7863608.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,5 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2008_S7863608.hh" +#include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" @@ -8,95 +8,122 @@ namespace Rivet { - D0_2008_S7863608::D0_2008_S7863608() - : Analysis("D0_2008_S7863608") - { - setBeams(PROTON, ANTIPROTON); - setNeedsCrossSection(true); - - ZFinder zfinder(-1.7, 1.7, 15.0*GeV, MUON, 65.0*GeV, 115.0*GeV, 0.2); - addProjection(zfinder, "ZFinder"); + /// @brief Measurement differential Z/gamma* + jet +X cross sections + /// @author Gavin Hesketh, Andy Buckley, Frank Siegert + class D0_2008_S7863608 : public Analysis { + + public: + + /// @name Construction + //@{ + /// Constructor + D0_2008_S7863608() + : Analysis("D0_2008_S7863608") + { + setBeams(PROTON, ANTIPROTON); + setNeedsCrossSection(true); + + ZFinder zfinder(-1.7, 1.7, 15.0*GeV, MUON, 65.0*GeV, 115.0*GeV, 0.2); + addProjection(zfinder, "ZFinder"); + + FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5, 20.0*GeV); + addProjection(conefinder, "ConeFinder"); + } - FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5, 20.0*GeV); - addProjection(conefinder, "ConeFinder"); - } - + //@} - // Book histograms - void D0_2008_S7863608::init() { - - _h_jet_pT_cross_section = bookHistogram1D(1, 1, 1); - _h_jet_y_cross_section = bookHistogram1D(2, 1, 1); - _h_Z_pT_cross_section = bookHistogram1D(3, 1, 1); - _h_Z_y_cross_section = bookHistogram1D(4, 1, 1); - _h_total_cross_section = bookHistogram1D(5, 1, 1); + /// @name Analysis methods + //@{ + + /// Book histograms + void init() { + _h_jet_pT_cross_section = bookHistogram1D(1, 1, 1); + _h_jet_y_cross_section = bookHistogram1D(2, 1, 1); + _h_Z_pT_cross_section = bookHistogram1D(3, 1, 1); + _h_Z_y_cross_section = bookHistogram1D(4, 1, 1); + _h_total_cross_section = bookHistogram1D(5, 1, 1); + } + - } - - - - // Do the analysis - void D0_2008_S7863608::analyze(const Event & e) { - double weight = e.weight(); - const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size()==1) { - const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder"); - const Jets& jets = jetpro.jetsByPt(20.0*GeV); - Jets jets_cut; - foreach (const Jet& j, jets) { - if (fabs(j.momentum().pseudorapidity()) < 2.8) { - jets_cut.push_back(j); - } - } - - // Return if there are no jets: - if(jets_cut.size()<1) { - getLog() << Log::DEBUG << "Skipping event " << e.genEvent().event_number() - << " because no jets pass cuts " << endl; - vetoEvent; - } + // Do the analysis + void analyze(const Event& e) { + const double weight = e.weight(); - // cut on Delta R between jet and muons - foreach (const Jet& j, jets_cut) { - foreach (const Particle& mu, zfinder.constituentsFinalState().particles()) { - if (deltaR(mu.momentum().pseudorapidity(), mu.momentum().azimuthalAngle(), - j.momentum().pseudorapidity(), j.momentum().azimuthalAngle()) < 0.5) { - vetoEvent; + const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); + if (zfinder.particles().size()==1) { + const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder"); + const Jets& jets = jetpro.jetsByPt(20.0*GeV); + Jets jets_cut; + foreach (const Jet& j, jets) { + if (fabs(j.momentum().pseudorapidity()) < 2.8) { + jets_cut.push_back(j); } } + + // Return if there are no jets: + if(jets_cut.size()<1) { + getLog() << Log::DEBUG << "Skipping event " << e.genEvent().event_number() + << " because no jets pass cuts " << endl; + vetoEvent; + } + + // cut on Delta R between jet and muons + foreach (const Jet& j, jets_cut) { + foreach (const Particle& mu, zfinder.constituentsFinalState().particles()) { + if (deltaR(mu.momentum().pseudorapidity(), mu.momentum().azimuthalAngle(), + j.momentum().pseudorapidity(), j.momentum().azimuthalAngle()) < 0.5) { + vetoEvent; + } + } + } + + const FourMomentum Zmom = zfinder.particles()[0].momentum(); + + // In jet pT + _h_jet_pT_cross_section->fill( jets_cut[0].momentum().pT(), weight); + _h_jet_y_cross_section->fill( fabs(jets_cut[0].momentum().rapidity()), weight); + + // In Z pT + _h_Z_pT_cross_section->fill(Zmom.pT(), weight); + _h_Z_y_cross_section->fill(fabs(Zmom.rapidity()), weight); + + _h_total_cross_section->fill(1960.0, weight); } - - const FourMomentum Zmom = zfinder.particles()[0].momentum(); - - // In jet pT - _h_jet_pT_cross_section->fill( jets_cut[0].momentum().pT(), weight); - _h_jet_y_cross_section->fill( fabs(jets_cut[0].momentum().rapidity()), weight); - - // In Z pT - _h_Z_pT_cross_section->fill(Zmom.pT(), weight); - _h_Z_y_cross_section->fill(fabs(Zmom.rapidity()), weight); - - _h_total_cross_section->fill(1960.0, weight); } - } - + + + + /// Finalize + void finalize() { + const double invlumi = crossSection()/sumOfWeights(); + scale(_h_total_cross_section, invlumi); + scale(_h_jet_pT_cross_section, invlumi); + scale(_h_jet_y_cross_section, invlumi); + scale(_h_Z_pT_cross_section, invlumi); + scale(_h_Z_y_cross_section, invlumi); + } + + //@} - // Finalize - void D0_2008_S7863608::finalize() { - const double invlumi = crossSection()/sumOfWeights(); - scale(_h_total_cross_section, invlumi); - scale(_h_jet_pT_cross_section, invlumi); - scale(_h_jet_y_cross_section, invlumi); - scale(_h_Z_pT_cross_section, invlumi); - scale(_h_Z_y_cross_section, invlumi); - } + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D * _h_jet_pT_cross_section; + AIDA::IHistogram1D * _h_jet_y_cross_section; + AIDA::IHistogram1D * _h_Z_pT_cross_section; + AIDA::IHistogram1D * _h_Z_y_cross_section; + AIDA::IHistogram1D * _h_total_cross_section; + //@} + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2008_S7863608> plugin_D0_2008_S7863608; - + } Modified: trunk/src/Analyses/D0_2009_S8202443.cc ============================================================================== --- trunk/src/Analyses/D0_2009_S8202443.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2009_S8202443.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,5 +1,5 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2009_S8202443.hh" +#include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" @@ -9,119 +9,151 @@ namespace Rivet { - D0_2009_S8202443::D0_2009_S8202443() - : Analysis("D0_2009_S8202443"), - _sum_of_weights(0.0), _sum_of_weights_constrained(0.0) - { - setBeams(PROTON, ANTIPROTON); - - // Leptons in constrained tracking acceptance - std::vector<std::pair<double, double> > etaRanges; - etaRanges.push_back(make_pair(-2.5, -1.5)); - etaRanges.push_back(make_pair(-1.1, 1.1)); - etaRanges.push_back(make_pair(1.5, 2.5)); - ZFinder zfinder_constrained(etaRanges, 25.0*GeV, ELECTRON, - 65.0*GeV, 115.0*GeV, 0.2); - addProjection(zfinder_constrained, "ZFinderConstrained"); - FastJets conefinder_constrained(zfinder_constrained.remainingFinalState(), - FastJets::D0ILCONE, 0.5, 20.0*GeV); - addProjection(conefinder_constrained, "ConeFinderConstrained"); - - // Unconstrained leptons - ZFinder zfinder(FinalState(), ELECTRON, 65.0*GeV, 115.0*GeV, 0.2); - addProjection(zfinder, "ZFinder"); - FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5, 20.0*GeV); - addProjection(conefinder, "ConeFinder"); - } - - - // Book histograms - void D0_2009_S8202443::init() { - _h_jet1_pT_constrained = bookHistogram1D(1, 1, 1); - _h_jet2_pT_constrained = bookHistogram1D(3, 1, 1); - _h_jet3_pT_constrained = bookHistogram1D(5, 1, 1); - _h_jet1_pT = bookHistogram1D(2, 1, 1); - _h_jet2_pT = bookHistogram1D(4, 1, 1); - _h_jet3_pT = bookHistogram1D(6, 1, 1); - } + class D0_2009_S8202443 : public Analysis { + public: - - // Do the analysis - void D0_2009_S8202443::analyze(const Event & e) { - double weight = e.weight(); - - // unconstrained electrons first - const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); - if (zfinder.particles().size()==1) { - _sum_of_weights += weight; - const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder"); - const Jets& jets = jetpro.jetsByPt(20.0*GeV); - Jets jets_cut; - foreach (const Jet& j, jets) { - if (fabs(j.momentum().pseudorapidity()) < 2.5) { - jets_cut.push_back(j); - } - } + /// @name Construction + //@{ + /// Constructor + D0_2009_S8202443() + : Analysis("D0_2009_S8202443"), + _sum_of_weights(0.0), _sum_of_weights_constrained(0.0) + { + setBeams(PROTON, ANTIPROTON); - if (jets_cut.size()>0) { - _h_jet1_pT->fill(jets_cut[0].momentum().pT()/GeV, weight); - } - if (jets_cut.size()>1) { - _h_jet2_pT->fill(jets_cut[1].momentum().pT()/GeV, weight); - } - if (jets_cut.size()>2) { - _h_jet3_pT->fill(jets_cut[2].momentum().pT()/GeV, weight); - } - } - else { - getLog() << Log::DEBUG << "no unique lepton pair found." << endl; + // Leptons in constrained tracking acceptance + vector<pair<double, double> > etaRanges; + etaRanges.push_back(make_pair(-2.5, -1.5)); + etaRanges.push_back(make_pair(-1.1, 1.1)); + etaRanges.push_back(make_pair(1.5, 2.5)); + ZFinder zfinder_constrained(etaRanges, 25.0*GeV, ELECTRON, + 65.0*GeV, 115.0*GeV, 0.2); + addProjection(zfinder_constrained, "ZFinderConstrained"); + FastJets conefinder_constrained(zfinder_constrained.remainingFinalState(), + FastJets::D0ILCONE, 0.5, 20.0*GeV); + addProjection(conefinder_constrained, "ConeFinderConstrained"); + + // Unconstrained leptons + ZFinder zfinder(FinalState(), ELECTRON, 65.0*GeV, 115.0*GeV, 0.2); + addProjection(zfinder, "ZFinder"); + FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5, 20.0*GeV); + addProjection(conefinder, "ConeFinder"); + } + + //@} + + + /// @name Analysis methods + //@{ + + /// Book histograms + void init() { + _h_jet1_pT_constrained = bookHistogram1D(1, 1, 1); + _h_jet2_pT_constrained = bookHistogram1D(3, 1, 1); + _h_jet3_pT_constrained = bookHistogram1D(5, 1, 1); + _h_jet1_pT = bookHistogram1D(2, 1, 1); + _h_jet2_pT = bookHistogram1D(4, 1, 1); + _h_jet3_pT = bookHistogram1D(6, 1, 1); } - - - // constrained electrons - const ZFinder& zfinder_constrained = applyProjection<ZFinder>(e, "ZFinderConstrained"); - if (zfinder_constrained.particles().size()==1) { - _sum_of_weights_constrained += weight; - const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinderConstrained"); - const Jets& jets = jetpro.jetsByPt(20.0*GeV); - Jets jets_cut; - foreach (const Jet& j, jets) { - if (fabs(j.momentum().pseudorapidity()) < 2.5) { - jets_cut.push_back(j); + + + + // Do the analysis + void analyze(const Event& e) { + double weight = e.weight(); + + // unconstrained electrons first + const ZFinder& zfinder = applyProjection<ZFinder>(e, "ZFinder"); + if (zfinder.particles().size()==1) { + _sum_of_weights += weight; + const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinder"); + const Jets& jets = jetpro.jetsByPt(20.0*GeV); + Jets jets_cut; + foreach (const Jet& j, jets) { + if (fabs(j.momentum().pseudorapidity()) < 2.5) { + jets_cut.push_back(j); + } + } + + if (jets_cut.size()>0) { + _h_jet1_pT->fill(jets_cut[0].momentum().pT()/GeV, weight); + } + if (jets_cut.size()>1) { + _h_jet2_pT->fill(jets_cut[1].momentum().pT()/GeV, weight); + } + if (jets_cut.size()>2) { + _h_jet3_pT->fill(jets_cut[2].momentum().pT()/GeV, weight); } } - - if (jets_cut.size()>0) { - _h_jet1_pT_constrained->fill(jets_cut[0].momentum().pT()/GeV, weight); + else { + getLog() << Log::DEBUG << "no unique lepton pair found." << endl; } - if (jets_cut.size()>1) { - _h_jet2_pT_constrained->fill(jets_cut[1].momentum().pT()/GeV, weight); + + + // constrained electrons + const ZFinder& zfinder_constrained = applyProjection<ZFinder>(e, "ZFinderConstrained"); + if (zfinder_constrained.particles().size()==1) { + _sum_of_weights_constrained += weight; + const JetAlg& jetpro = applyProjection<JetAlg>(e, "ConeFinderConstrained"); + const Jets& jets = jetpro.jetsByPt(20.0*GeV); + Jets jets_cut; + foreach (const Jet& j, jets) { + if (fabs(j.momentum().pseudorapidity()) < 2.5) { + jets_cut.push_back(j); + } + } + + if (jets_cut.size()>0) { + _h_jet1_pT_constrained->fill(jets_cut[0].momentum().pT()/GeV, weight); + } + if (jets_cut.size()>1) { + _h_jet2_pT_constrained->fill(jets_cut[1].momentum().pT()/GeV, weight); + } + if (jets_cut.size()>2) { + _h_jet3_pT_constrained->fill(jets_cut[2].momentum().pT()/GeV, weight); + } } - if (jets_cut.size()>2) { - _h_jet3_pT_constrained->fill(jets_cut[2].momentum().pT()/GeV, weight); + else { + getLog() << Log::DEBUG << "no unique lepton pair found." << endl; + vetoEvent; } } - else { - getLog() << Log::DEBUG << "no unique lepton pair found." << endl; - vetoEvent; + + + + // Finalize + void finalize() { + scale(_h_jet1_pT, 1.0/_sum_of_weights); + scale(_h_jet2_pT, 1.0/_sum_of_weights); + scale(_h_jet3_pT, 1.0/_sum_of_weights); + scale(_h_jet1_pT_constrained, 1.0/_sum_of_weights_constrained); + scale(_h_jet2_pT_constrained, 1.0/_sum_of_weights_constrained); + scale(_h_jet3_pT_constrained, 1.0/_sum_of_weights_constrained); } - } + + //@} + private: - // Finalize - void D0_2009_S8202443::finalize() { - scale(_h_jet1_pT, 1.0/_sum_of_weights); - scale(_h_jet2_pT, 1.0/_sum_of_weights); - scale(_h_jet3_pT, 1.0/_sum_of_weights); - scale(_h_jet1_pT_constrained, 1.0/_sum_of_weights_constrained); - scale(_h_jet2_pT_constrained, 1.0/_sum_of_weights_constrained); - scale(_h_jet3_pT_constrained, 1.0/_sum_of_weights_constrained); - } + /// @name Histograms + //@{ + AIDA::IHistogram1D * _h_jet1_pT; + AIDA::IHistogram1D * _h_jet2_pT; + AIDA::IHistogram1D * _h_jet3_pT; + AIDA::IHistogram1D * _h_jet1_pT_constrained; + AIDA::IHistogram1D * _h_jet2_pT_constrained; + AIDA::IHistogram1D * _h_jet3_pT_constrained; + //@} + + double _sum_of_weights, _sum_of_weights_constrained; + }; + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2009_S8202443> plugin_D0_2009_S8202443; - + } Modified: trunk/src/Analyses/D0_2009_S8320160.cc ============================================================================== --- trunk/src/Analyses/D0_2009_S8320160.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2009_S8320160.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,6 +1,7 @@ // -*- C++ -*- -#include "Rivet/Analyses/D0_2009_S8320160.hh" +#include "Rivet/Analysis.hh" #include "Rivet/Tools/Logging.hh" +#include "Rivet/Tools/BinnedHistogram.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/FastJets.hh" @@ -9,65 +10,89 @@ namespace Rivet { - D0_2009_S8320160::D0_2009_S8320160() - : Analysis("D0_2009_S8320160") - { - setBeams(PROTON, ANTIPROTON); - - FinalState fs; - FastJets conefinder(fs, FastJets::D0ILCONE, 0.7); - addProjection(conefinder, "ConeFinder"); - } - + class D0_2009_S8320160 : public Analysis { - // Book histograms - void D0_2009_S8320160::init() { - _h_chi_dijet.addHistogram(250., 300., bookHistogram1D(1, 1, 1)); - _h_chi_dijet.addHistogram(300., 400., bookHistogram1D(2, 1, 1)); - _h_chi_dijet.addHistogram(400., 500., bookHistogram1D(3, 1, 1)); - _h_chi_dijet.addHistogram(500., 600., bookHistogram1D(4, 1, 1)); - _h_chi_dijet.addHistogram(600., 700., bookHistogram1D(5, 1, 1)); - _h_chi_dijet.addHistogram(700., 800., bookHistogram1D(6, 1, 1)); - _h_chi_dijet.addHistogram(800., 900., bookHistogram1D(7, 1, 1)); - _h_chi_dijet.addHistogram(900., 1000., bookHistogram1D(8, 1, 1)); - _h_chi_dijet.addHistogram(1000., 1100., bookHistogram1D(9, 1, 1)); - _h_chi_dijet.addHistogram(1100., 1960, bookHistogram1D(10, 1, 1)); - } + public: + /// @name Construction + //@{ + /// Constructor + D0_2009_S8320160() + : Analysis("D0_2009_S8320160") + { + setBeams(PROTON, ANTIPROTON); + + FinalState fs; + FastJets conefinder(fs, FastJets::D0ILCONE, 0.7); + addProjection(conefinder, "ConeFinder"); + } + + //@} - // Do the analysis - void D0_2009_S8320160::analyze(const Event & e) { - double weight = e.weight(); - const Jets& jets = applyProjection<JetAlg>(e, "ConeFinder").jetsByPt(); + /// @name Analysis methods + //@{ - if(jets.size()<2) vetoEvent; + // Book histograms + void init() { + _h_chi_dijet.addHistogram(250., 300., bookHistogram1D(1, 1, 1)); + _h_chi_dijet.addHistogram(300., 400., bookHistogram1D(2, 1, 1)); + _h_chi_dijet.addHistogram(400., 500., bookHistogram1D(3, 1, 1)); + _h_chi_dijet.addHistogram(500., 600., bookHistogram1D(4, 1, 1)); + _h_chi_dijet.addHistogram(600., 700., bookHistogram1D(5, 1, 1)); + _h_chi_dijet.addHistogram(700., 800., bookHistogram1D(6, 1, 1)); + _h_chi_dijet.addHistogram(800., 900., bookHistogram1D(7, 1, 1)); + _h_chi_dijet.addHistogram(900., 1000., bookHistogram1D(8, 1, 1)); + _h_chi_dijet.addHistogram(1000., 1100., bookHistogram1D(9, 1, 1)); + _h_chi_dijet.addHistogram(1100., 1960, bookHistogram1D(10, 1, 1)); + } - FourMomentum j0(jets[0].momentum()); - FourMomentum j1(jets[1].momentum()); - double y0 = j0.rapidity(); - double y1 = j1.rapidity(); - if (fabs(y0+y1)>2) vetoEvent; - double mjj = FourMomentum(j0+j1).mass(); - double chi = exp(fabs(y0-y1)); - _h_chi_dijet.fill(mjj, chi, weight); - } - - - - // Finalize - void D0_2009_S8320160::finalize() { - foreach (AIDA::IHistogram1D* hist, _h_chi_dijet.getHistograms()) { - normalize(hist); + /// Do the analysis + void analyze(const Event & e) { + const double weight = e.weight(); + + const Jets& jets = applyProjection<JetAlg>(e, "ConeFinder").jetsByPt(); + if (jets.size() < 2) vetoEvent; + + FourMomentum j0(jets[0].momentum()); + FourMomentum j1(jets[1].momentum()); + double y0 = j0.rapidity(); + double y1 = j1.rapidity(); + + if (fabs(y0+y1)>2) vetoEvent; + + double mjj = FourMomentum(j0+j1).mass(); + double chi = exp(fabs(y0-y1)); + _h_chi_dijet.fill(mjj, chi, weight); + } + + + + /// Finalize + void finalize() { + foreach (AIDA::IHistogram1D* hist, _h_chi_dijet.getHistograms()) { + normalize(hist); + } } - } + //@} + + + private: + + /// @name Histograms + //@{ + BinnedHistogram<double> _h_chi_dijet; + //@} + + }; + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2009_S8320160> plugin_D0_2009_S8320160; - + } Modified: trunk/src/Analyses/D0_2009_S8349509.cc ============================================================================== --- trunk/src/Analyses/D0_2009_S8349509.cc Mon Aug 31 08:25:18 2009 (r1784) +++ trunk/src/Analyses/D0_2009_S8349509.cc Mon Aug 31 09:50:19 2009 (r1785) @@ -1,8 +1,7 @@ // -*- C++ -*- -#include "Rivet/Rivet.hh" +#include "Rivet/Analysis.hh" #include "Rivet/RivetAIDA.hh" #include "Rivet/Tools/Logging.hh" -#include "Rivet/Analyses/D0_2009_S8349509.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" @@ -10,101 +9,139 @@ namespace Rivet { - D0_2009_S8349509::D0_2009_S8349509() : Analysis("D0_2009_S8349509") { - setBeams(PROTON, ANTIPROTON); + class D0_2009_S8349509 : public Analysis { + public: - ZFinder zfinder(-1.7, 1.7, 15.0*GeV, MUON, 65.0*GeV, 115.0*GeV, 0.2); - addProjection(zfinder, "ZFinder"); + /// @name Constructors etc. + //@{ - FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5, 20.0*GeV); - addProjection(conefinder, "ConeFinder"); - } - - - void D0_2009_S8349509::init() { - - _h_dphi_jet_Z25 = bookHistogram1D(1, 1, 1); - _h_dphi_jet_Z45 = bookHistogram1D(2, 1, 1); - - _h_dy_jet_Z25 = bookHistogram1D(3, 1, 1); - _h_dy_jet_Z45 = bookHistogram1D(4, 1, 1); - - _h_yboost_jet_Z25 = bookHistogram1D(5, 1, 1); - _h_yboost_jet_Z45 = bookHistogram1D(6, 1, 1); + /// Constructor + D0_2009_S8349509() : + Analysis("D0_2009_S8349509") + { + setBeams(PROTON, ANTIPROTON); + + ZFinder zfinder(-1.7, 1.7, 15.0*GeV, MUON, 65.0*GeV, 115.0*GeV, 0.2); + addProjection(zfinder, "ZFinder"); + + FastJets conefinder(zfinder.remainingFinalState(), FastJets::D0ILCONE, 0.5, 20.0*GeV); + addProjection(conefinder, "ConeFinder"); + } - _inclusive_Z_sumofweights = 0.0; - } - + //@} - void D0_2009_S8349509::analyze(const Event& event) { - double weight = event.weight(); - - const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder"); - if (zfinder.particles().size()==1) { - // count inclusive sum of weights for histogram normalisation - _inclusive_Z_sumofweights += weight; + /// @name Analysis methods + //@{ + + /// Book histograms + void init() { - Jets jets; - foreach (const Jet& j, applyProjection<JetAlg>(event, "ConeFinder").jetsByPt()) { - if (fabs(j.momentum().pseudorapidity()) < 2.8) { - jets.push_back(j); - break; - } - } + _h_dphi_jet_Z25 = bookHistogram1D(1, 1, 1); + _h_dphi_jet_Z45 = bookHistogram1D(2, 1, 1); - // Return if there are no jets: - if(jets.size()<1) { - getLog() << Log::DEBUG << "Skipping event " << event.genEvent().event_number() - << " because no jets pass cuts " << endl; - vetoEvent; - } + _h_dy_jet_Z25 = bookHistogram1D(3, 1, 1); + _h_dy_jet_Z45 = bookHistogram1D(4, 1, 1); + + _h_yboost_jet_Z25 = bookHistogram1D(5, 1, 1); + _h_yboost_jet_Z45 = bookHistogram1D(6, 1, 1); + + _inclusive_Z_sumofweights = 0.0; + } + + + void analyze(const Event& event) { + const double weight = event.weight(); - // cut on Delta R between jet and muons - foreach (const Jet& j, jets) { - foreach (const Particle& mu, zfinder.constituentsFinalState().particles()) { - if (deltaR(mu.momentum(), j.momentum()) < 0.5) { - vetoEvent; + const ZFinder& zfinder = applyProjection<ZFinder>(event, "ZFinder"); + if (zfinder.particles().size()==1) { + // count inclusive sum of weights for histogram normalisation + _inclusive_Z_sumofweights += weight; + + Jets jets; + foreach (const Jet& j, applyProjection<JetAlg>(event, "ConeFinder").jetsByPt()) { + if (fabs(j.momentum().pseudorapidity()) < 2.8) { + jets.push_back(j); + break; } } + + // Return if there are no jets: + if (jets.size() < 1) { + getLog() << Log::DEBUG << "Skipping event " << event.genEvent().event_number() + << " because no jets pass cuts " << endl; + vetoEvent; + } + + // Cut on Delta R between jet and muons + foreach (const Jet& j, jets) { + foreach (const Particle& mu, zfinder.constituentsFinalState().particles()) { + if (deltaR(mu.momentum(), j.momentum()) < 0.5) { + vetoEvent; + } + } + } + + const FourMomentum Zmom = zfinder.particles()[0].momentum(); + const FourMomentum jetmom = jets[0].momentum(); + double yZ = Zmom.rapidity(); + double yjet = jetmom.rapidity(); + double dphi = deltaPhi(Zmom.phi(), jetmom.phi()); + double dy = fabs(yZ-yjet); + double yboost = fabs(yZ+yjet)/2.0; + + if (Zmom.pT() > 25.0*GeV) { + _h_dphi_jet_Z25->fill(dphi,weight); + _h_dy_jet_Z25->fill(dy, weight); + _h_yboost_jet_Z25->fill(yboost, weight); + } + if (Zmom.pT() > 45.0*GeV) { + _h_dphi_jet_Z45->fill(dphi,weight); + _h_dy_jet_Z45->fill(dy, weight); + _h_yboost_jet_Z45->fill(yboost, weight); + } } - const FourMomentum Zmom = zfinder.particles()[0].momentum(); - const FourMomentum jetmom = jets[0].momentum(); - double yZ = Zmom.rapidity(); - double yjet = jetmom.rapidity(); - double dphi = deltaPhi(Zmom.phi(), jetmom.phi()); - double dy = fabs(yZ-yjet); - double yboost = fabs(yZ+yjet)/2.0; - - if (Zmom.pT()>25.0*GeV) { - _h_dphi_jet_Z25->fill(dphi,weight); - _h_dy_jet_Z25->fill(dy, weight); - _h_yboost_jet_Z25->fill(yboost, weight); - } - if (Zmom.pT()>45.0*GeV) { - _h_dphi_jet_Z45->fill(dphi,weight); - _h_dy_jet_Z45->fill(dy, weight); - _h_yboost_jet_Z45->fill(yboost, weight); - } } + + + void finalize() { + if (_inclusive_Z_sumofweights == 0.0) return; + scale(_h_dphi_jet_Z25, 1.0/_inclusive_Z_sumofweights); + scale(_h_dphi_jet_Z45, 1.0/_inclusive_Z_sumofweights); + scale(_h_dy_jet_Z25, 1.0/_inclusive_Z_sumofweights); + scale(_h_dy_jet_Z45, 1.0/_inclusive_Z_sumofweights); + scale(_h_yboost_jet_Z25, 1.0/_inclusive_Z_sumofweights); + scale(_h_yboost_jet_Z45, 1.0/_inclusive_Z_sumofweights); + } + + //@} - } + private: + // Data members like post-cuts event weight counters go here - void D0_2009_S8349509::finalize() { - if (_inclusive_Z_sumofweights==0.0) return; - scale(_h_dphi_jet_Z25, 1.0/_inclusive_Z_sumofweights); - scale(_h_dphi_jet_Z45, 1.0/_inclusive_Z_sumofweights); - scale(_h_dy_jet_Z25, 1.0/_inclusive_Z_sumofweights); - scale(_h_dy_jet_Z45, 1.0/_inclusive_Z_sumofweights); - scale(_h_yboost_jet_Z25, 1.0/_inclusive_Z_sumofweights); - scale(_h_yboost_jet_Z45, 1.0/_inclusive_Z_sumofweights); - } + private: + /// @name Histograms + //@{ + AIDA::IHistogram1D *_h_dphi_jet_Z25; + AIDA::IHistogram1D *_h_dphi_jet_Z45; + AIDA::IHistogram1D *_h_dy_jet_Z25; + AIDA::IHistogram1D *_h_dy_jet_Z45; + + AIDA::IHistogram1D *_h_yboost_jet_Z25; + AIDA::IHistogram1D *_h_yboost_jet_Z45; + //@} + + double _inclusive_Z_sumofweights; + }; + + + // This global object acts as a hook for the plugin system AnalysisBuilder<D0_2009_S8349509> plugin_D0_2009_S8349509; - + }
More information about the Rivet-svn mailing list |