ATLAS_2012_I943401.cc
Go to the documentation of this file.
00001 // -*- C++ -*- 00002 #include "Rivet/Analysis.hh" 00003 #include "Rivet/Tools/BinnedHistogram.hh" 00004 #include "Rivet/Projections/FinalState.hh" 00005 #include "Rivet/Projections/ChargedFinalState.hh" 00006 #include "Rivet/Projections/VisibleFinalState.hh" 00007 #include "Rivet/Projections/IdentifiedFinalState.hh" 00008 #include "Rivet/Projections/FastJets.hh" 00009 #include "Rivet/Projections/VetoedFinalState.hh" 00010 00011 namespace Rivet { 00012 00013 00014 class ATLAS_2012_I943401 : public Analysis { 00015 public: 00016 00017 /// @name Constructors etc. 00018 //@{ 00019 00020 /// Constructor 00021 00022 ATLAS_2012_I943401() 00023 : Analysis("ATLAS_2012_I943401") 00024 { } 00025 00026 //@} 00027 00028 00029 public: 00030 00031 /// @name Analysis methods 00032 //@{ 00033 00034 /// Book histograms and initialise projections before the run 00035 void init() { 00036 00037 // projection to find the electrons 00038 IdentifiedFinalState elecs(EtaIn(-2.47, 2.47) & (Cuts::pT >= 20.0*GeV)); 00039 elecs.acceptIdPair(PID::ELECTRON); 00040 addProjection(elecs, "elecs"); 00041 00042 // projection to find the muons 00043 IdentifiedFinalState muons(EtaIn(-2.4, 2.4) & (Cuts::pT >= 10.0*GeV)); 00044 muons.acceptIdPair(PID::MUON); 00045 addProjection(muons, "muons"); 00046 00047 // jet finder 00048 VetoedFinalState vfs; 00049 vfs.addVetoPairId(PID::MUON); 00050 addProjection(FastJets(vfs, FastJets::ANTIKT, 0.4), 00051 "AntiKtJets04"); 00052 00053 // all tracks (to do deltaR with leptons) 00054 addProjection(ChargedFinalState(-3.0,3.0,0.5*GeV),"cfs"); 00055 00056 // for pTmiss 00057 addProjection(VisibleFinalState(EtaIn(-4.5, 4.5)),"vfs"); 00058 00059 // book histograms 00060 00061 // counts in signal regions 00062 _count_OS_SR1 = bookHisto1D("count_OS_SR1", 1, 0., 1.); 00063 _count_OS_SR2 = bookHisto1D("count_OS_SR2", 1, 0., 1.); 00064 _count_OS_SR3 = bookHisto1D("count_OS_SR3", 1, 0., 1.); 00065 _count_SS_SR1 = bookHisto1D("count_SS_SR1", 1, 0., 1.); 00066 _count_SS_SR2 = bookHisto1D("count_SS_SR2", 1, 0., 1.); 00067 _count_FS_SR1 = bookHisto1D("count_FS_SR1", 1, 0., 1.); 00068 _count_FS_SR2 = bookHisto1D("count_FS_SR2", 1, 0., 1.); 00069 _count_FS_SR3 = bookHisto1D("count_FS_SR3", 1, 0., 1.); 00070 00071 // histograms from paper 00072 00073 _hist_mll_SS_D = bookHisto1D( 1,1,1); 00074 _hist_mll_SS_B = bookHisto1D( 1,1,2); 00075 _hist_eTmiss_SS_D = bookHisto1D( 2,1,1); 00076 _hist_eTmiss_SS_B = bookHisto1D( 2,1,2); 00077 _hist_mll_SS_2Jet_D = bookHisto1D( 3,1,1); 00078 _hist_mll_SS_2Jet_B = bookHisto1D( 3,1,2); 00079 _hist_njet_SS_D = bookHisto1D( 5,1,1); 00080 _hist_njet_SS_B = bookHisto1D( 5,1,2); 00081 _hist_pT_j1_SS_D = bookHisto1D( 6,1,1); 00082 _hist_pT_j1_SS_B = bookHisto1D( 6,1,2); 00083 _hist_pT_j2_SS_D = bookHisto1D( 7,1,1); 00084 _hist_pT_j2_SS_B = bookHisto1D( 7,1,2); 00085 _hist_pT_l1_SS_D = bookHisto1D( 8,1,1); 00086 _hist_pT_l1_SS_B = bookHisto1D( 8,1,2); 00087 _hist_pT_l2_SS_D = bookHisto1D( 9,1,1); 00088 _hist_pT_l2_SS_B = bookHisto1D( 9,1,2); 00089 _hist_mll_OS_D = bookHisto1D(10,1,1); 00090 _hist_mll_OS_B = bookHisto1D(10,1,2); 00091 _hist_eTmiss_OS_D = bookHisto1D(11,1,1); 00092 _hist_eTmiss_OS_B = bookHisto1D(11,1,2); 00093 _hist_eTmiss_3Jet_OS_D = bookHisto1D(12,1,1); 00094 _hist_eTmiss_3Jet_OS_B = bookHisto1D(12,1,2); 00095 _hist_eTmiss_4Jet_OS_D = bookHisto1D(13,1,1); 00096 _hist_eTmiss_4Jet_OS_B = bookHisto1D(13,1,2); 00097 _hist_njet_OS_D = bookHisto1D(14,1,1); 00098 _hist_njet_OS_B = bookHisto1D(14,1,2); 00099 _hist_pT_j1_OS_D = bookHisto1D(15,1,1); 00100 _hist_pT_j1_OS_B = bookHisto1D(15,1,2); 00101 _hist_pT_j2_OS_D = bookHisto1D(16,1,1); 00102 _hist_pT_j2_OS_B = bookHisto1D(16,1,2); 00103 _hist_pT_l1_OS_D = bookHisto1D(17,1,1); 00104 _hist_pT_l1_OS_B = bookHisto1D(17,1,2); 00105 _hist_pT_l2_OS_D = bookHisto1D(18,1,1); 00106 _hist_pT_l2_OS_B = bookHisto1D(18,1,2); 00107 //???? 00108 // <dataPointSet name="d04-x01-y01" dimension="2" path="/REF/ATLAS_2011_I943401" title="EVENTS/10 GEV" > 00109 // <dataPointSet name="d04-x01-y02" dimension="2" path="/REF/ATLAS_2011_I943401" title="EVENTS/10 GEV" > 00110 } 00111 00112 /// Perform the event analysis 00113 void analyze(const Event& event) { 00114 // event weight 00115 const double weight = event.weight(); 00116 00117 // get the jet candidates 00118 Jets cand_jets; 00119 foreach (const Jet& jet, 00120 applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV) ) { 00121 if ( fabs( jet.eta() ) < 2.8 ) { 00122 cand_jets.push_back(jet); 00123 } 00124 } 00125 00126 // electron candidates 00127 Particles cand_e = 00128 applyProjection<IdentifiedFinalState>(event, "elecs").particlesByPt(); 00129 00130 // Discard jets that overlap with electrons 00131 Jets recon_jets; 00132 foreach ( const Jet& jet, cand_jets ) { 00133 bool away_from_e = true; 00134 foreach ( const Particle & e, cand_e ) { 00135 if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) { 00136 away_from_e = false; 00137 break; 00138 } 00139 } 00140 if ( away_from_e ) recon_jets.push_back( jet ); 00141 } 00142 // get the charged tracks for isolation 00143 Particles chg_tracks = 00144 applyProjection<ChargedFinalState>(event, "cfs").particles(); 00145 00146 // Reconstructed electrons 00147 Particles recon_e; 00148 foreach ( const Particle & e, cand_e ) { 00149 // check not near a jet 00150 bool e_near_jet = false; 00151 foreach ( const Jet& jet, recon_jets ) { 00152 if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) { 00153 e_near_jet = true; 00154 break; 00155 } 00156 } 00157 if ( e_near_jet ) continue; 00158 // check the isolation 00159 double pTinCone = -e.pT(); 00160 foreach ( const Particle & track, chg_tracks ) { 00161 if ( deltaR(e.momentum(),track.momentum()) < 0.2 ) 00162 pTinCone += track.pT(); 00163 } 00164 if ( pTinCone < 0.1*e.momentum().perp() ) 00165 recon_e.push_back(e); 00166 } 00167 00168 // Reconstructed Muons 00169 Particles recon_mu; 00170 Particles cand_mu = 00171 applyProjection<IdentifiedFinalState>(event,"muons").particlesByPt(); 00172 foreach ( const Particle & mu, cand_mu ) { 00173 // check not near a jet 00174 bool mu_near_jet = false; 00175 foreach ( const Jet& jet, recon_jets ) { 00176 if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) { 00177 mu_near_jet = true; 00178 break; 00179 } 00180 } 00181 if ( mu_near_jet ) continue; 00182 // isolation 00183 double pTinCone = -mu.pT(); 00184 foreach ( const Particle & track, chg_tracks ) { 00185 if ( deltaR(mu.momentum(),track.momentum()) < 0.2 ) 00186 pTinCone += track.pT(); 00187 } 00188 if ( pTinCone < 1.8*GeV ) 00189 recon_mu.push_back(mu); 00190 } 00191 00192 // pTmiss 00193 Particles vfs_particles 00194 = applyProjection<VisibleFinalState>(event, "vfs").particles(); 00195 FourMomentum pTmiss; 00196 foreach ( const Particle & p, vfs_particles ) { 00197 pTmiss -= p.momentum(); 00198 } 00199 double eTmiss = pTmiss.pT(); 00200 00201 // ATLAS calo problem 00202 if(rand()/static_cast<double>(RAND_MAX)<=0.42) { 00203 foreach ( const Particle & e, recon_e ) { 00204 double eta = e.eta(); 00205 double phi = e.momentum().azimuthalAngle(MINUSPI_PLUSPI); 00206 if(eta>-0.1&&eta<1.5&&phi>-0.9&&phi<-0.5) 00207 vetoEvent; 00208 } 00209 foreach ( const Jet & jet, recon_jets ) { 00210 double eta = jet.rapidity(); 00211 double phi = jet.momentum().azimuthalAngle(MINUSPI_PLUSPI); 00212 if(jet.momentum().perp()>40 && eta>-0.1&&eta<1.5&&phi>-0.9&&phi<-0.5) 00213 vetoEvent; 00214 } 00215 } 00216 00217 // Exactly two leptons for each event 00218 if ( recon_mu.size() + recon_e.size() != 2) 00219 vetoEvent; 00220 // two electrons highest pT > 25 00221 Particles recon_leptons; 00222 if(recon_e.size()==2&&recon_e[0].momentum().perp()>25.) { 00223 recon_leptons = recon_e; 00224 } 00225 // two muons highest pT > 20 00226 else if(recon_mu.size()==2&&recon_mu[0].momentum().perp()>20.) { 00227 recon_leptons = recon_mu; 00228 } 00229 else if(recon_e.size()==1 && recon_mu.size()==1 && 00230 (recon_e[0].momentum().perp()>25. ||recon_mu[0].momentum().perp()>20. )) { 00231 if(recon_mu[0].momentum().perp()<recon_e[0].momentum().perp()) { 00232 recon_leptons.push_back(recon_e [0]); 00233 recon_leptons.push_back(recon_mu[0]); 00234 } 00235 else { 00236 recon_leptons.push_back(recon_mu[0]); 00237 recon_leptons.push_back(recon_e [0]); 00238 } 00239 } 00240 // fails trigger 00241 else 00242 vetoEvent; 00243 00244 double mll = (recon_leptons[0].momentum()+recon_leptons[1].momentum()).mass(); 00245 // lepton pair mass > 12. 00246 if(mll < 12.) vetoEvent; 00247 00248 // same sign or opposite sign event 00249 int sign = recon_leptons[0].pdgId()*recon_leptons[1].pdgId(); 00250 00251 // same sign leptons 00252 if(sign>0) { 00253 _hist_mll_SS_D ->fill(mll ,weight); 00254 _hist_mll_SS_B ->fill(mll ,weight); 00255 _hist_eTmiss_SS_D->fill(eTmiss,weight); 00256 _hist_eTmiss_SS_B->fill(eTmiss,weight); 00257 if(recon_jets.size()>=2) { 00258 _hist_mll_SS_2Jet_D ->fill(mll ,weight); 00259 _hist_mll_SS_2Jet_B ->fill(mll ,weight); 00260 } 00261 _hist_njet_SS_D ->fill(recon_jets.size(),weight); 00262 _hist_njet_SS_B ->fill(recon_jets.size(),weight); 00263 if(!recon_jets.empty()) { 00264 _hist_pT_j1_SS_D->fill(recon_jets[0].momentum().perp(),weight); 00265 _hist_pT_j1_SS_B->fill(recon_jets[0].momentum().perp(),weight); 00266 } 00267 if(recon_jets.size()>2) { 00268 _hist_pT_j2_SS_D->fill(recon_jets[1].momentum().perp(),weight); 00269 _hist_pT_j2_SS_B->fill(recon_jets[1].momentum().perp(),weight); 00270 } 00271 _hist_pT_l1_SS_D->fill(recon_leptons[0].momentum().perp(),weight); 00272 _hist_pT_l1_SS_B->fill(recon_leptons[0].momentum().perp(),weight); 00273 _hist_pT_l2_SS_D->fill(recon_leptons[1].momentum().perp(),weight); 00274 _hist_pT_l2_SS_B->fill(recon_leptons[1].momentum().perp(),weight); 00275 // SS-SR1 00276 if(eTmiss>100.) { 00277 _count_SS_SR1->fill(0.5,weight); 00278 } 00279 // SS-SR2 00280 if(eTmiss>80. && recon_jets.size()>=2 && 00281 recon_jets[1].momentum().perp()>50.) { 00282 _count_SS_SR2->fill(0.5,weight); 00283 } 00284 } 00285 // opposite sign 00286 else { 00287 _hist_mll_OS_D->fill(mll ,weight); 00288 _hist_mll_OS_B->fill(mll ,weight); 00289 _hist_eTmiss_OS_D->fill(eTmiss,weight); 00290 _hist_eTmiss_OS_B->fill(eTmiss,weight); 00291 if(recon_jets.size()>=3){ 00292 _hist_eTmiss_3Jet_OS_D->fill(eTmiss,weight); 00293 _hist_eTmiss_3Jet_OS_B->fill(eTmiss,weight); 00294 } 00295 if(recon_jets.size()>=4){ 00296 _hist_eTmiss_4Jet_OS_D->fill(eTmiss,weight); 00297 _hist_eTmiss_4Jet_OS_B->fill(eTmiss,weight); 00298 } 00299 _hist_njet_OS_D->fill(recon_jets.size(),weight); 00300 _hist_njet_OS_B->fill(recon_jets.size(),weight); 00301 if(!recon_jets.empty()) { 00302 _hist_pT_j1_OS_D->fill(recon_jets[0].momentum().perp(),weight); 00303 _hist_pT_j1_OS_B->fill(recon_jets[0].momentum().perp(),weight); 00304 } 00305 if(recon_jets.size()>2) { 00306 _hist_pT_j2_OS_D->fill(recon_jets[1].momentum().perp(),weight); 00307 _hist_pT_j2_OS_B->fill(recon_jets[1].momentum().perp(),weight); 00308 } 00309 _hist_pT_l1_OS_D->fill(recon_leptons[0].momentum().perp(),weight); 00310 _hist_pT_l1_OS_B->fill(recon_leptons[0].momentum().perp(),weight); 00311 _hist_pT_l2_OS_D->fill(recon_leptons[1].momentum().perp(),weight); 00312 _hist_pT_l2_OS_B->fill(recon_leptons[1].momentum().perp(),weight); 00313 // different signal regions 00314 // OS-SR1 00315 if(eTmiss>250.) { 00316 _count_OS_SR1->fill(0.5,weight); 00317 } 00318 // OS-SR2 00319 if(eTmiss>220. && recon_jets.size()>=3 && 00320 recon_jets[0].momentum().perp()>80. && 00321 recon_jets[2].momentum().perp()>40.) { 00322 _count_OS_SR2->fill(0.5,weight); 00323 } 00324 // OS-SR3 00325 if(eTmiss>100. && recon_jets.size()>=4 && 00326 recon_jets[0].momentum().perp()>100. && 00327 recon_jets[3].momentum().perp()>70.) { 00328 _count_OS_SR3->fill(0.5,weight); 00329 } 00330 // same flavour analysis 00331 static const double beta = 0.75; 00332 static const double tau_e = 0.96; 00333 static const double tau_mu = 0.816; 00334 double fs_weight = weight; 00335 if(abs(recon_leptons[0].pdgId())==PID::ELECTRON && abs(recon_leptons[1].pdgId())==PID::ELECTRON) { 00336 fs_weight /= beta*(1.-sqr(1.-tau_e)); 00337 } 00338 else if(abs(recon_leptons[0].pdgId())==PID::MUON && abs(recon_leptons[1].pdgId())==PID::MUON) { 00339 fs_weight *= beta/(1.-sqr(1.-tau_mu)); 00340 } 00341 else { 00342 fs_weight /= -(1.-(1.-tau_e)*(1.-tau_mu)); 00343 } 00344 // FS-SR1 00345 if(eTmiss>80.&& (mll<80.||mll>100.)) { 00346 _count_FS_SR1->fill(0.5,fs_weight); 00347 } 00348 // FS-SR2 00349 if(eTmiss>80.&&recon_jets.size()>=2) { 00350 _count_FS_SR2->fill(0.5,fs_weight); 00351 } 00352 // FS-SR3 00353 if(eTmiss>250.) { 00354 _count_FS_SR3->fill(0.5,fs_weight); 00355 } 00356 } 00357 } 00358 00359 //@} 00360 00361 00362 void finalize() { 00363 00364 double norm = crossSection()/femtobarn*1.04/sumOfWeights(); 00365 // event counts 00366 scale(_count_OS_SR1,norm); 00367 scale(_count_OS_SR2,norm); 00368 scale(_count_OS_SR3,norm); 00369 scale(_count_SS_SR1,norm); 00370 scale(_count_SS_SR2,norm); 00371 scale(_count_FS_SR1,norm); 00372 scale(_count_FS_SR2,norm); 00373 scale(_count_FS_SR3,norm); 00374 // histograms 00375 scale(_hist_mll_SS_D ,norm*20.); 00376 scale(_hist_mll_SS_B ,norm*20.); 00377 scale(_hist_eTmiss_SS_D ,norm*20.); 00378 scale(_hist_eTmiss_SS_B ,norm*20.); 00379 scale(_hist_mll_SS_2Jet_D,norm*50.); 00380 scale(_hist_mll_SS_2Jet_B,norm*50.); 00381 scale(_hist_njet_SS_D ,norm ); 00382 scale(_hist_njet_SS_B ,norm ); 00383 scale(_hist_pT_j1_SS_D ,norm*20.); 00384 scale(_hist_pT_j1_SS_B ,norm*20.); 00385 scale(_hist_pT_j2_SS_D ,norm*20.); 00386 scale(_hist_pT_j2_SS_B ,norm*20.); 00387 scale(_hist_pT_l1_SS_D ,norm*5. ); 00388 scale(_hist_pT_l1_SS_B ,norm*5. ); 00389 scale(_hist_pT_l2_SS_D ,norm*5. ); 00390 scale(_hist_pT_l2_SS_B ,norm*5. ); 00391 00392 scale(_hist_mll_OS_D ,norm*10.); 00393 scale(_hist_mll_OS_B ,norm*10.); 00394 scale(_hist_eTmiss_OS_D ,norm*10.); 00395 scale(_hist_eTmiss_OS_B ,norm*10.); 00396 scale(_hist_eTmiss_3Jet_OS_D,norm*10.); 00397 scale(_hist_eTmiss_3Jet_OS_B,norm*10.); 00398 scale(_hist_eTmiss_4Jet_OS_D,norm*10.); 00399 scale(_hist_eTmiss_4Jet_OS_B,norm*10.); 00400 scale(_hist_njet_OS_D ,norm ); 00401 scale(_hist_njet_OS_B ,norm ); 00402 scale(_hist_pT_j1_OS_D ,norm*20.); 00403 scale(_hist_pT_j1_OS_B ,norm*20.); 00404 scale(_hist_pT_j2_OS_D ,norm*20.); 00405 scale(_hist_pT_j2_OS_B ,norm*20.); 00406 scale(_hist_pT_l1_OS_D ,norm*20.); 00407 scale(_hist_pT_l1_OS_B ,norm*20.); 00408 scale(_hist_pT_l2_OS_D ,norm*20.); 00409 scale(_hist_pT_l2_OS_B ,norm*20.); 00410 } 00411 00412 private: 00413 00414 /// @name Histograms 00415 //@{ 00416 Histo1DPtr _count_OS_SR1; 00417 Histo1DPtr _count_OS_SR2; 00418 Histo1DPtr _count_OS_SR3; 00419 Histo1DPtr _count_SS_SR1; 00420 Histo1DPtr _count_SS_SR2; 00421 Histo1DPtr _count_FS_SR1; 00422 Histo1DPtr _count_FS_SR2; 00423 Histo1DPtr _count_FS_SR3; 00424 00425 Histo1DPtr _hist_mll_SS_D; 00426 Histo1DPtr _hist_mll_SS_B; 00427 Histo1DPtr _hist_eTmiss_SS_D; 00428 Histo1DPtr _hist_eTmiss_SS_B; 00429 Histo1DPtr _hist_mll_SS_2Jet_D; 00430 Histo1DPtr _hist_mll_SS_2Jet_B; 00431 Histo1DPtr _hist_njet_SS_D; 00432 Histo1DPtr _hist_njet_SS_B; 00433 Histo1DPtr _hist_pT_j1_SS_D; 00434 Histo1DPtr _hist_pT_j1_SS_B; 00435 Histo1DPtr _hist_pT_j2_SS_D; 00436 Histo1DPtr _hist_pT_j2_SS_B; 00437 Histo1DPtr _hist_pT_l1_SS_D; 00438 Histo1DPtr _hist_pT_l1_SS_B; 00439 Histo1DPtr _hist_pT_l2_SS_D; 00440 Histo1DPtr _hist_pT_l2_SS_B; 00441 00442 Histo1DPtr _hist_mll_OS_D; 00443 Histo1DPtr _hist_mll_OS_B; 00444 Histo1DPtr _hist_eTmiss_OS_D; 00445 Histo1DPtr _hist_eTmiss_OS_B; 00446 Histo1DPtr _hist_eTmiss_3Jet_OS_D; 00447 Histo1DPtr _hist_eTmiss_3Jet_OS_B; 00448 Histo1DPtr _hist_eTmiss_4Jet_OS_D; 00449 Histo1DPtr _hist_eTmiss_4Jet_OS_B; 00450 Histo1DPtr _hist_njet_OS_D ; 00451 Histo1DPtr _hist_njet_OS_B ; 00452 Histo1DPtr _hist_pT_j1_OS_D; 00453 Histo1DPtr _hist_pT_j1_OS_B; 00454 Histo1DPtr _hist_pT_j2_OS_D; 00455 Histo1DPtr _hist_pT_j2_OS_B; 00456 Histo1DPtr _hist_pT_l1_OS_D; 00457 Histo1DPtr _hist_pT_l1_OS_B; 00458 Histo1DPtr _hist_pT_l2_OS_D; 00459 Histo1DPtr _hist_pT_l2_OS_B; 00460 //@} 00461 }; 00462 00463 // The hook for the plugin system 00464 DECLARE_RIVET_PLUGIN(ATLAS_2012_I943401); 00465 00466 } Generated on Thu Feb 6 2014 17:38:42 for The Rivet MC analysis system by 
1.7.6.1
|