00001
00002 #include "Rivet/Analysis.hh"
00003 #include "Rivet/RivetAIDA.hh"
00004 #include "Rivet/Tools/Logging.hh"
00005 #include "Rivet/Projections/FastJets.hh"
00006 #include "Rivet/Projections/VetoedFinalState.hh"
00007 #include "Rivet/Projections/VisibleFinalState.hh"
00008 #include "Rivet/Projections/MissingMomentum.hh"
00009
00010 namespace Rivet {
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 class D0_2004_S5992206 : public Analysis {
00027
00028 public:
00029
00030
00031
00032
00033
00034 D0_2004_S5992206() : Analysis("D0_2004_S5992206")
00035 {
00036 setBeams(PROTON, ANTIPROTON);
00037 }
00038
00039
00040
00041
00042
00043
00044
00045 void init() {
00046
00047 const FinalState fs(-3.0, 3.0);
00048 addProjection(fs, "FS");
00049
00050 VetoedFinalState vfs(fs);
00051 vfs.vetoNeutrinos();
00052 vfs.addVetoPairDetail(MUON, 1.0*GeV, MAXDOUBLE);
00053 addProjection(vfs, "VFS");
00054 addProjection(FastJets(vfs, FastJets::D0ILCONE, 0.7), "Jets");
00055 addProjection(MissingMomentum(vfs), "CalMET");
00056
00057
00058 _histJetAzimuth_pTmax75_100 = bookHistogram1D(1, 2, 1);
00059 _histJetAzimuth_pTmax100_130 = bookHistogram1D(2, 2, 1);
00060 _histJetAzimuth_pTmax130_180 = bookHistogram1D(3, 2, 1);
00061 _histJetAzimuth_pTmax180_ = bookHistogram1D(4, 2, 1);
00062 }
00063
00064
00065
00066 void analyze(const Event& event) {
00067
00068
00069 const JetAlg& jetpro = applyProjection<JetAlg>(event, "Jets");
00070 getLog() << Log::DEBUG << "Jet multiplicity before any pT cut = " << jetpro.size() << endl;
00071
00072 const Jets jets = jetpro.jetsByPt(40.0*GeV);
00073 if (jets.size() >= 2) {
00074 getLog() << Log::DEBUG << "Jet multiplicity after pT > 40 GeV cut = " << jets.size() << endl;
00075 } else {
00076 vetoEvent;
00077 }
00078 const double rap1 = jets[0].momentum().rapidity();
00079 const double rap2 = jets[1].momentum().rapidity();
00080 if (fabs(rap1) > 0.5 || fabs(rap2) > 0.5) {
00081 vetoEvent;
00082 }
00083 getLog() << Log::DEBUG << "Jet eta and pT requirements fulfilled" << endl;
00084 const double pT1 = jets[0].momentum().pT();
00085
00086 const MissingMomentum& caloMissEt = applyProjection<MissingMomentum>(event, "CalMET");
00087 getLog() << Log::DEBUG << "Missing vector Et = " << caloMissEt.vectorET()/GeV << " GeV" << endl;
00088 if (caloMissEt.vectorET() > 0.7*pT1) {
00089 MSG_DEBUG("Vetoing event with too much missing Et: "
00090 << caloMissEt.vectorET()/GeV << " GeV > "
00091 << 0.7*pT1/GeV << " GeV");
00092 vetoEvent;
00093 }
00094
00095 if (pT1/GeV >= 75.0) {
00096 const double weight = event.weight();
00097 const double dphi = deltaPhi(jets[0].momentum().phi(), jets[1].momentum().phi());
00098 if (inRange(pT1/GeV, 75.0, 100.0)) {
00099 _histJetAzimuth_pTmax75_100->fill(dphi, weight);
00100 } else if (inRange(pT1/GeV, 100.0, 130.0)) {
00101 _histJetAzimuth_pTmax100_130->fill(dphi, weight);
00102 } else if (inRange(pT1/GeV, 130.0, 180.0)) {
00103 _histJetAzimuth_pTmax130_180->fill(dphi, weight);
00104 } else if (pT1/GeV > 180.0) {
00105 _histJetAzimuth_pTmax180_->fill(dphi, weight);
00106 }
00107 }
00108
00109 }
00110
00111
00112
00113 void finalize() {
00114
00115 normalize(_histJetAzimuth_pTmax75_100);
00116 normalize(_histJetAzimuth_pTmax100_130);
00117 normalize(_histJetAzimuth_pTmax130_180);
00118 normalize(_histJetAzimuth_pTmax180_);
00119 }
00120
00121
00122
00123
00124 private:
00125
00126
00127
00128 AIDA::IHistogram1D* _histJetAzimuth_pTmax75_100;
00129 AIDA::IHistogram1D* _histJetAzimuth_pTmax100_130;
00130 AIDA::IHistogram1D* _histJetAzimuth_pTmax130_180;
00131 AIDA::IHistogram1D* _histJetAzimuth_pTmax180_;
00132
00133
00134 };
00135
00136
00137
00138
00139 AnalysisBuilder<D0_2004_S5992206> plugin_D0_2004_S5992206;
00140
00141 }