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