ATLAS_2011_I945498.cc
Go to the documentation of this file.
00001 // -*- C++ -*- 00002 #include "Rivet/Analysis.hh" 00003 00004 #include "Rivet/Projections/ZFinder.hh" 00005 #include "Rivet/Projections/FastJets.hh" 00006 #include "Rivet/Projections/FinalState.hh" 00007 #include "Rivet/Projections/VetoedFinalState.hh" 00008 #include "Rivet/Projections/IdentifiedFinalState.hh" 00009 #include "Rivet/Projections/LeadingParticlesFinalState.hh" 00010 00011 00012 namespace Rivet { 00013 00014 00015 /// ATLAS Z+jets in pp at 7 TeV 00016 class ATLAS_2011_I945498 : public Analysis { 00017 public: 00018 00019 /// Constructor 00020 ATLAS_2011_I945498() 00021 : Analysis("ATLAS_2011_I945498") 00022 { } 00023 00024 00025 /// Book histograms and initialise projections before the run 00026 void init() { 00027 00028 // Variable initialisation 00029 _isZeeSample = false; 00030 _isZmmSample = false; 00031 for (size_t chn = 0; chn < 3; ++chn) { 00032 weights_nj0[chn] = 0; 00033 weights_nj1[chn] = 0; 00034 weights_nj2[chn] = 0; 00035 weights_nj3[chn] = 0; 00036 weights_nj4[chn] = 0; 00037 } 00038 00039 // Set up projections 00040 FinalState fs; 00041 ZFinder zfinder_mu(fs, EtaIn(-2.4, 2.4) & (Cuts::pT >= 20*GeV), PID::MUON, 66*GeV, 116*GeV, 0.1, ZFinder::CLUSTERNODECAY); 00042 addProjection(zfinder_mu, "ZFinder_mu"); 00043 00044 Cut cuts = ( EtaIn(-2.47, -1.52) 00045 | EtaIn(-1.37, 1.37) 00046 | EtaIn( 1.52, 2.47) ) & (Cuts::pT >= 20.0*GeV); 00047 00048 ZFinder zfinder_el(fs, cuts, PID::ELECTRON, 66*GeV, 116*GeV, 0.1, ZFinder::CLUSTERNODECAY); 00049 addProjection(zfinder_el, "ZFinder_el"); 00050 00051 Cut cuts25_20 = EtaIn(-2.5,2.5) & (Cuts::pT >= 20.0*GeV); 00052 // For combined cross-sections (combined phase space + dressed level) 00053 ZFinder zfinder_comb_mu(fs, cuts25_20, PID::MUON, 66.0*GeV, 116.0*GeV, 0.1, ZFinder::CLUSTERNODECAY); 00054 addProjection(zfinder_comb_mu, "ZFinder_comb_mu"); 00055 ZFinder zfinder_comb_el(fs, cuts25_20, PID::ELECTRON, 66.0*GeV, 116.0*GeV, 0.1, ZFinder::CLUSTERNODECAY); 00056 addProjection(zfinder_comb_el, "ZFinder_comb_el"); 00057 00058 // Define veto FS in order to prevent Z-decay products entering the jet algorithm 00059 VetoedFinalState remfs; 00060 remfs.addVetoOnThisFinalState(zfinder_el); 00061 remfs.addVetoOnThisFinalState(zfinder_mu); 00062 VetoedFinalState remfs_comb; 00063 remfs_comb.addVetoOnThisFinalState(zfinder_comb_el); 00064 remfs_comb.addVetoOnThisFinalState(zfinder_comb_mu); 00065 00066 FastJets jets(remfs, FastJets::ANTIKT, 0.4); 00067 jets.useInvisibles(); 00068 addProjection(jets, "jets"); 00069 FastJets jets_comb(remfs_comb, FastJets::ANTIKT, 0.4); 00070 jets_comb.useInvisibles(); 00071 addProjection(jets_comb, "jets_comb"); 00072 00073 // 0=el, 1=mu, 2=comb 00074 for (size_t chn = 0; chn < 3; ++chn) { 00075 _h_njet_incl[chn] = bookHisto1D(1, 1, chn+1); 00076 _h_njet_ratio[chn] = bookScatter2D(2, 1, chn+1); 00077 _h_ptjet[chn] = bookHisto1D(3, 1, chn+1); 00078 _h_ptlead[chn] = bookHisto1D(4, 1, chn+1); 00079 _h_ptseclead[chn] = bookHisto1D(5, 1, chn+1); 00080 _h_yjet[chn] = bookHisto1D(6, 1, chn+1); 00081 _h_ylead[chn] = bookHisto1D(7, 1, chn+1); 00082 _h_yseclead[chn] = bookHisto1D(8, 1, chn+1); 00083 _h_mass[chn] = bookHisto1D(9, 1, chn+1); 00084 _h_deltay[chn] = bookHisto1D(10, 1, chn+1); 00085 _h_deltaphi[chn] = bookHisto1D(11, 1, chn+1); 00086 _h_deltaR[chn] = bookHisto1D(12, 1, chn+1); 00087 } 00088 } 00089 00090 00091 // Jet selection criteria universal for electron and muon channel 00092 /// @todo Replace with a Cut passed to jetsByPt 00093 Jets selectJets(const ZFinder* zf, const FastJets* allJets) { 00094 const FourMomentum l1 = zf->constituents()[0].momentum(); 00095 const FourMomentum l2 = zf->constituents()[1].momentum(); 00096 Jets jets; 00097 foreach (const Jet& jet, allJets->jetsByPt(30*GeV)) { 00098 const FourMomentum jmom = jet.momentum(); 00099 if (fabs(jmom.rapidity()) < 4.4 && 00100 deltaR(l1, jmom) > 0.5 && deltaR(l2, jmom) > 0.5) { 00101 jets.push_back(jet); 00102 } 00103 } 00104 return jets; 00105 } 00106 00107 00108 /// Perform the per-event analysis 00109 void analyze(const Event& event) { 00110 const double weight = event.weight(); 00111 00112 vector<const ZFinder*> zfs; 00113 zfs.push_back(& (applyProjection<ZFinder>(event, "ZFinder_el"))); 00114 zfs.push_back(& (applyProjection<ZFinder>(event, "ZFinder_mu"))); 00115 zfs.push_back(& (applyProjection<ZFinder>(event, "ZFinder_comb_el"))); 00116 zfs.push_back(& (applyProjection<ZFinder>(event, "ZFinder_comb_mu"))); 00117 00118 vector<const FastJets*> fjs; 00119 fjs.push_back(& (applyProjection<FastJets>(event, "jets"))); 00120 fjs.push_back(& (applyProjection<FastJets>(event, "jets_comb"))); 00121 00122 // Determine what kind of MC sample this is 00123 const bool isZee = (zfs[0]->bosons().size() == 1) || (zfs[2]->bosons().size() == 1); 00124 const bool isZmm = (zfs[1]->bosons().size() == 1) || (zfs[3]->bosons().size() == 1); 00125 if (isZee) _isZeeSample = true; 00126 if (isZmm) _isZmmSample = true; 00127 00128 // Require exactly one electronic or muonic Z-decay in the event 00129 bool isZeemm = ( (zfs[0]->bosons().size() == 1 && zfs[1]->bosons().size() != 1) || 00130 (zfs[1]->bosons().size() == 1 && zfs[0]->bosons().size() != 1) ); 00131 bool isZcomb = ( (zfs[2]->bosons().size() == 1 && zfs[3]->bosons().size() != 1) || 00132 (zfs[3]->bosons().size() == 1 && zfs[2]->bosons().size() != 1) ); 00133 if (!isZeemm && !isZcomb) vetoEvent; 00134 00135 vector<int> zfIDs; 00136 vector<int> fjIDs; 00137 if (isZeemm) { 00138 int chn = zfs[0]->bosons().size() == 1 ? 0 : 1; 00139 zfIDs.push_back(chn); 00140 fjIDs.push_back(0); 00141 } 00142 if (isZcomb) { 00143 int chn = zfs[2]->bosons().size() == 1 ? 2 : 3; 00144 zfIDs.push_back(chn); 00145 fjIDs.push_back(1); 00146 } 00147 00148 for (size_t izf = 0; izf < zfIDs.size(); ++izf) { 00149 int zfID = zfIDs[izf]; 00150 int fjID = fjIDs[izf]; 00151 00152 int chn = zfID; 00153 if (zfID == 2 || zfID == 3) chn = 2; 00154 00155 Jets jets = selectJets(zfs[zfID], fjs[fjID]); 00156 00157 switch (jets.size()) { 00158 case 0: 00159 weights_nj0[chn] += weight; 00160 break; 00161 case 1: 00162 weights_nj0[chn] += weight; 00163 weights_nj1[chn] += weight; 00164 break; 00165 case 2: 00166 weights_nj0[chn] += weight; 00167 weights_nj1[chn] += weight; 00168 weights_nj2[chn] += weight; 00169 break; 00170 case 3: 00171 weights_nj0[chn] += weight; 00172 weights_nj1[chn] += weight; 00173 weights_nj2[chn] += weight; 00174 weights_nj3[chn] += weight; 00175 break; 00176 default: // >= 4 00177 weights_nj0[chn] += weight; 00178 weights_nj1[chn] += weight; 00179 weights_nj2[chn] += weight; 00180 weights_nj3[chn] += weight; 00181 weights_nj4[chn] += weight; 00182 } 00183 00184 // Require at least one jet 00185 if (jets.empty()) continue; 00186 00187 // Fill jet multiplicities 00188 for (size_t ijet = 1; ijet <= jets.size(); ++ijet) { 00189 _h_njet_incl[chn]->fill(ijet, weight); 00190 } 00191 00192 // Loop over selected jets, fill inclusive jet distributions 00193 for (size_t ijet = 0; ijet < jets.size(); ++ijet) { 00194 _h_ptjet[chn]->fill(jets[ijet].momentum().pT()/GeV, weight); 00195 _h_yjet [chn]->fill(fabs(jets[ijet].momentum().rapidity()), weight); 00196 } 00197 00198 // Leading jet histos 00199 const double ptlead = jets[0].momentum().pT()/GeV; 00200 const double yabslead = fabs(jets[0].momentum().rapidity()); 00201 _h_ptlead[chn]->fill(ptlead, weight); 00202 _h_ylead [chn]->fill(yabslead, weight); 00203 00204 if (jets.size() >= 2) { 00205 // Second jet histos 00206 const double pt2ndlead = jets[1].momentum().pT()/GeV; 00207 const double yabs2ndlead = fabs(jets[1].momentum().rapidity()); 00208 _h_ptseclead[chn] ->fill(pt2ndlead, weight); 00209 _h_yseclead [chn] ->fill(yabs2ndlead, weight); 00210 00211 // Dijet histos 00212 const double deltaphi = fabs(deltaPhi(jets[1], jets[0])); 00213 const double deltarap = fabs(jets[0].momentum().rapidity() - jets[1].momentum().rapidity()) ; 00214 const double deltar = fabs(deltaR(jets[0], jets[1], RAPIDITY)); 00215 const double mass = (jets[0].momentum() + jets[1].momentum()).mass(); 00216 _h_mass [chn] ->fill(mass/GeV, weight); 00217 _h_deltay [chn] ->fill(deltarap, weight); 00218 _h_deltaphi[chn] ->fill(deltaphi, weight); 00219 _h_deltaR [chn] ->fill(deltar, weight); 00220 } 00221 } 00222 } 00223 00224 00225 /// @name Ratio calculator util functions 00226 //@{ 00227 00228 /// Calculate the ratio, being careful about div-by-zero 00229 double ratio(double a, double b) { 00230 return (b != 0) ? a/b : 0; 00231 } 00232 00233 /// Calculate the ratio error, being careful about div-by-zero 00234 double ratio_err(double a, double b) { 00235 return (b != 0) ? sqrt(a/sqr(b) + sqr(a)/(b*b*b)) : 0; 00236 } 00237 00238 //@} 00239 00240 00241 void finalize() { 00242 // Fill ratio histograms 00243 for (size_t chn = 0; chn < 3; ++chn) { 00244 _h_njet_ratio[chn]->addPoint(1, ratio(weights_nj1[chn], weights_nj0[chn]), 0, ratio_err(weights_nj1[chn], weights_nj0[chn])); 00245 _h_njet_ratio[chn]->addPoint(2, ratio(weights_nj2[chn], weights_nj1[chn]), 0, ratio_err(weights_nj2[chn], weights_nj1[chn])); 00246 _h_njet_ratio[chn]->addPoint(3, ratio(weights_nj3[chn], weights_nj2[chn]), 0, ratio_err(weights_nj3[chn], weights_nj2[chn])); 00247 _h_njet_ratio[chn]->addPoint(4, ratio(weights_nj4[chn], weights_nj3[chn]), 0, ratio_err(weights_nj4[chn], weights_nj3[chn])); 00248 } 00249 00250 // Scale other histos 00251 for (size_t chn = 0; chn < 3; ++chn) { 00252 // For ee and mumu channels: normalize to Njet inclusive cross-section 00253 double xs = (chn == 2) ? crossSectionPerEvent()/picobarn : 1 / weights_nj0[chn]; 00254 // For inclusive MC sample(ee/mmu channels together) we want the single-lepton-flavor xsec 00255 if (_isZeeSample && _isZmmSample) xs /= 2; 00256 00257 // Special case histogram: always not normalized 00258 scale(_h_njet_incl[chn], (chn < 2) ? crossSectionPerEvent()/picobarn : xs); 00259 00260 scale(_h_ptjet[chn] , xs); 00261 scale(_h_ptlead[chn] , xs); 00262 scale(_h_ptseclead[chn], xs); 00263 scale(_h_yjet[chn] , xs); 00264 scale(_h_ylead[chn] , xs); 00265 scale(_h_yseclead[chn] , xs); 00266 scale(_h_deltaphi[chn] , xs); 00267 scale(_h_deltay[chn] , xs); 00268 scale(_h_deltaR[chn] , xs); 00269 scale(_h_mass[chn] , xs); 00270 } 00271 00272 } 00273 00274 //@} 00275 00276 00277 private: 00278 00279 bool _isZeeSample; 00280 bool _isZmmSample; 00281 00282 double weights_nj0[3]; 00283 double weights_nj1[3]; 00284 double weights_nj2[3]; 00285 double weights_nj3[3]; 00286 double weights_nj4[3]; 00287 00288 Scatter2DPtr _h_njet_ratio[3]; 00289 Histo1DPtr _h_njet_incl[3]; 00290 Histo1DPtr _h_ptjet[3]; 00291 Histo1DPtr _h_ptlead[3]; 00292 Histo1DPtr _h_ptseclead[3]; 00293 Histo1DPtr _h_yjet[3]; 00294 Histo1DPtr _h_ylead[3]; 00295 Histo1DPtr _h_yseclead[3]; 00296 Histo1DPtr _h_deltaphi[3]; 00297 Histo1DPtr _h_deltay[3]; 00298 Histo1DPtr _h_deltaR[3]; 00299 Histo1DPtr _h_mass[3]; 00300 00301 }; 00302 00303 00304 DECLARE_RIVET_PLUGIN(ATLAS_2011_I945498); 00305 00306 00307 } Generated on Thu Feb 6 2014 17:38:41 for The Rivet MC analysis system by 1.7.6.1 |