ATLAS_2012_CONF_2012_109.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 /// @author Peter Richardson 00014 class ATLAS_2012_CONF_2012_109 : public Analysis { 00015 public: 00016 00017 /// @name Constructors etc. 00018 //@{ 00019 00020 /// Constructor 00021 ATLAS_2012_CONF_2012_109() 00022 : Analysis("ATLAS_2012_CONF_2012_109") 00023 { } 00024 00025 //@} 00026 00027 00028 public: 00029 00030 /// @name Analysis methods 00031 //@{ 00032 00033 /// Book histograms and initialise projections before the run 00034 void init() { 00035 00036 // Projection to find the electrons 00037 IdentifiedFinalState elecs(EtaIn(-2.47, 2.47) 00038 & (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) 00044 & (Cuts::pT >= 10.0*GeV)); 00045 muons.acceptIdPair(PID::MUON); 00046 addProjection(muons, "muons"); 00047 00048 // Jet finder 00049 VetoedFinalState vfs; 00050 vfs.addVetoPairId(PID::MUON); 00051 addProjection(FastJets(vfs, FastJets::ANTIKT, 0.4), "AntiKtJets04"); 00052 00053 // All tracks (to do deltaR with leptons) 00054 addProjection(ChargedFinalState(-3.0,3.0),"cfs"); 00055 00056 // Used for pTmiss (N.B. the real 'vfs' extends beyond 4.5 to |eta| = 4.9) 00057 addProjection(VisibleFinalState(-4.5,4.5),"vfs"); 00058 00059 // Book histograms 00060 _count_A_tight = bookHisto1D("count_A_tight" , 1, 0., 1.); 00061 _count_A_medium = bookHisto1D("count_A_medium" , 1, 0., 1.); 00062 _count_A_loose = bookHisto1D("count_A_loose" , 1, 0., 1.); 00063 _count_B_tight = bookHisto1D("count_B_tight" , 1, 0., 1.); 00064 _count_B_medium = bookHisto1D("count_B_medium" , 1, 0., 1.); 00065 _count_C_tight = bookHisto1D("count_C_tight" , 1, 0., 1.); 00066 _count_C_medium = bookHisto1D("count_C_medium" , 1, 0., 1.); 00067 _count_C_loose = bookHisto1D("count_C_loose" , 1, 0., 1.); 00068 _count_D_tight = bookHisto1D("count_D_tight" , 1, 0., 1.); 00069 _count_E_tight = bookHisto1D("count_E_tight" , 1, 0., 1.); 00070 _count_E_medium = bookHisto1D("count_E_medium" , 1, 0., 1.); 00071 _count_E_loose = bookHisto1D("count_E_loose" , 1, 0., 1.); 00072 00073 _hist_meff_A_medium = bookHisto1D("meff_A_medium" , 40, 0., 4000.); 00074 _hist_meff_A_tight = bookHisto1D("meff_A_tight" , 40, 0., 4000.); 00075 _hist_meff_B_medium = bookHisto1D("meff_B_medium" , 40, 0., 4000.); 00076 _hist_meff_B_tight = bookHisto1D("meff_B_tight" , 40, 0., 4000.); 00077 _hist_meff_C_medium = bookHisto1D("meff_C_medium" , 40, 0., 4000.); 00078 _hist_meff_C_tight = bookHisto1D("meff_C_tight" , 40, 0., 4000.); 00079 _hist_meff_D = bookHisto1D("meff_D" , 40, 0., 4000.); 00080 _hist_meff_E_loose = bookHisto1D("meff_E_loose" , 40, 0., 4000.); 00081 _hist_meff_E_medium = bookHisto1D("meff_E_medium" , 40, 0., 4000.); 00082 _hist_meff_E_tight = bookHisto1D("meff_E_tight" , 40, 0., 4000.); 00083 00084 } 00085 00086 00087 /// Perform the per-event analysis 00088 void analyze(const Event& event) { 00089 const double weight = event.weight(); 00090 00091 Jets cand_jets; 00092 const Jets jets = applyProjection<FastJets>(event, "AntiKtJets04").jetsByPt(20.0*GeV); 00093 foreach (const Jet& jet, jets) { 00094 if ( fabs( jet.eta() ) < 4.9 ) { 00095 cand_jets.push_back(jet); 00096 } 00097 } 00098 00099 const Particles cand_e = applyProjection<IdentifiedFinalState>(event, "elecs").particlesByPt(); 00100 00101 // Muon isolation not mentioned in hep-exp 1109.6572 but assumed to still be applicable 00102 Particles cand_mu; 00103 const Particles chg_tracks = applyProjection<ChargedFinalState>(event, "cfs").particles(); 00104 const Particles muons = applyProjection<IdentifiedFinalState>(event, "muons").particlesByPt(); 00105 foreach (const Particle& mu, muons) { 00106 double pTinCone = -mu.pT(); 00107 foreach (const Particle& track, chg_tracks) { 00108 if ( deltaR(mu.momentum(),track.momentum()) <= 0.2 ) { 00109 pTinCone += track.pT(); 00110 } 00111 } 00112 if ( pTinCone < 1.8*GeV ) cand_mu.push_back(mu); 00113 } 00114 00115 // Resolve jet-lepton overlap for jets with |eta| < 2.8 00116 Jets recon_jets; 00117 foreach ( const Jet& jet, cand_jets ) { 00118 if ( fabs( jet.eta() ) >= 2.8 ) continue; 00119 bool away_from_e = true; 00120 foreach ( const Particle & e, cand_e ) { 00121 if ( deltaR(e.momentum(),jet.momentum()) <= 0.2 ) { 00122 away_from_e = false; 00123 break; 00124 } 00125 } 00126 if ( away_from_e ) recon_jets.push_back( jet ); 00127 } 00128 00129 Particles recon_e, recon_mu; 00130 00131 foreach ( const Particle & e, cand_e ) { 00132 bool away = true; 00133 foreach ( const Jet& jet, recon_jets ) { 00134 if ( deltaR(e.momentum(),jet.momentum()) < 0.4 ) { 00135 away = false; 00136 break; 00137 } 00138 } 00139 if ( away ) recon_e.push_back( e ); 00140 } 00141 00142 foreach ( const Particle & mu, cand_mu ) { 00143 bool away = true; 00144 foreach ( const Jet& jet, recon_jets ) { 00145 if ( deltaR(mu.momentum(),jet.momentum()) < 0.4 ) { 00146 away = false; 00147 break; 00148 } 00149 } 00150 if ( away ) recon_mu.push_back( mu ); 00151 } 00152 00153 // pTmiss 00154 // Based on all candidate electrons, muons and jets, plus everything else with |eta| < 4.5 00155 // i.e. everything in our projection "vfs" plus the jets with |eta| > 4.5 00156 Particles vfs_particles = applyProjection<VisibleFinalState>(event, "vfs").particles(); 00157 FourMomentum pTmiss; 00158 foreach ( const Particle & p, vfs_particles ) { 00159 pTmiss -= p.momentum(); 00160 } 00161 foreach ( const Jet& jet, cand_jets ) { 00162 if ( fabs( jet.eta() ) > 4.5 ) pTmiss -= jet.momentum(); 00163 } 00164 double eTmiss = pTmiss.pT(); 00165 00166 // no electron pT> 20 or muons pT>10 00167 if ( !recon_mu.empty() || !recon_e.empty() ) { 00168 MSG_DEBUG("Charged leptons left after selection"); 00169 vetoEvent; 00170 } 00171 00172 if ( eTmiss <= 160 * GeV ) { 00173 MSG_DEBUG("Not enough eTmiss: " << eTmiss << " < 130"); 00174 vetoEvent; 00175 } 00176 00177 // check the hardest two jets 00178 if ( recon_jets.size()<2 || 00179 recon_jets[0].pT() <= 130.0 * GeV || 00180 recon_jets[0].pT() <= 60.0 * GeV ) { 00181 MSG_DEBUG("No hard leading jet in " << recon_jets.size() << " jets"); 00182 vetoEvent; 00183 } 00184 00185 // check the charged and EM fractions of the hard jets to avoid photons 00186 for (unsigned int ix = 0; ix < 2; ++ix) { 00187 // jets over 100 GeV 00188 if (recon_jets[ix].pT() < 100*GeV || 00189 recon_jets[ix].eta() > 2.) continue; //< @todo Should be |eta|? 00190 double fch(0.), fem(0.), eTotal(0.); 00191 foreach(const Particle & part, recon_jets[ix].particles()) { 00192 long id = abs(part.pdgId()); 00193 if(PID::threeCharge(id)!=0) 00194 fch += part.momentum().E(); 00195 if (id == PID::PHOTON || id == PID::ELECTRON || id == PID::PI0) 00196 fem += part.momentum().E(); 00197 } 00198 fch /= eTotal; 00199 fem /= eTotal; 00200 // remove events with hard photon 00201 if (fch < 0.02 || (fch < 0.05 && fem > 0.09)) vetoEvent; 00202 } 00203 00204 // ==================== observables ==================== 00205 00206 int Njets = 0; 00207 double min_dPhi_All = 999.999; //< @todo Use std::numeric_limits! 00208 double min_dPhi_2 = 999.999; //< @todo Use std::numeric_limits! 00209 double min_dPhi_3 = 999.999; //< @todo Use std::numeric_limits! 00210 double pTmiss_phi = pTmiss.phi(); 00211 00212 foreach ( const Jet& jet, recon_jets ) { 00213 if ( jet.pT() < 40*GeV ) continue; 00214 double dPhi = deltaPhi( pTmiss_phi, jet.momentum().phi()); 00215 if ( Njets < 2 ) min_dPhi_2 = min( min_dPhi_2, dPhi ); 00216 if ( Njets < 3 ) min_dPhi_3 = min( min_dPhi_3, dPhi ); 00217 min_dPhi_All = min( min_dPhi_All, dPhi ); 00218 ++Njets; 00219 } 00220 00221 // inclusive meff 00222 double m_eff_inc = eTmiss; 00223 foreach ( const Jet& jet, recon_jets ) { 00224 double perp = jet.pT(); 00225 if(perp>40.) m_eff_inc += perp; 00226 } 00227 00228 // region A 00229 double m_eff_Nj = eTmiss + recon_jets[0].pT() + recon_jets[1].pT(); 00230 if( min_dPhi_2 > 0.4 && eTmiss/m_eff_Nj > 0.3 ) { 00231 _hist_meff_A_tight ->fill(m_eff_inc,weight); 00232 if(eTmiss/m_eff_Nj > 0.4) 00233 _hist_meff_A_medium->fill(m_eff_inc,weight); 00234 if(m_eff_inc>1900.) 00235 _count_A_tight ->fill(0.5,weight); 00236 if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.4) 00237 _count_A_medium->fill(0.5,weight); 00238 if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.4) 00239 _count_A_loose ->fill(0.5,weight); 00240 } 00241 00242 // for rest of regions 3 jets pT> 60 needed 00243 if(recon_jets.size()<3 || recon_jets[2].momentum().perp()<60.) 00244 vetoEvent; 00245 00246 // region B 00247 m_eff_Nj += recon_jets[2].momentum().perp(); 00248 if( min_dPhi_3 > 0.4 && eTmiss/m_eff_Nj > 0.25 ) { 00249 _hist_meff_B_tight->fill(m_eff_inc,weight); 00250 if(eTmiss/m_eff_Nj > 0.3) 00251 _hist_meff_B_medium->fill(m_eff_inc,weight); 00252 if(m_eff_inc>1900.) 00253 _count_B_tight ->fill(0.5,weight); 00254 if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.3) 00255 _count_B_medium->fill(0.5,weight); 00256 } 00257 00258 // for rest of regions 4 jets pT> 60 needed 00259 if(recon_jets.size()<4 || recon_jets[3].momentum().perp()<60.) 00260 vetoEvent; 00261 00262 // region C 00263 m_eff_Nj += recon_jets[3].momentum().perp(); 00264 if( min_dPhi_3 > 0.4 && min_dPhi_All > 0.2 && eTmiss/m_eff_Nj > 0.25 ) { 00265 _hist_meff_C_tight->fill(m_eff_inc,weight); 00266 if( eTmiss/m_eff_Nj > 0.3 ) 00267 _hist_meff_C_medium->fill(m_eff_inc,weight); 00268 if(m_eff_inc>1900.) 00269 _count_C_tight ->fill(0.5,weight); 00270 if(m_eff_inc>1300. && eTmiss/m_eff_Nj > 0.3) 00271 _count_C_medium->fill(0.5,weight); 00272 if(m_eff_inc>1000. && eTmiss/m_eff_Nj > 0.3) 00273 _count_C_loose ->fill(0.5,weight); 00274 } 00275 00276 // for rest of regions 5 jets pT> 40 needed 00277 if(recon_jets.size()<5 || recon_jets[4].momentum().perp()<40.) 00278 vetoEvent; 00279 00280 // region D 00281 m_eff_Nj += recon_jets[4].momentum().perp(); 00282 if( min_dPhi_3 > 0.4 && min_dPhi_All > 0.2 && eTmiss/m_eff_Nj > 0.15 ) { 00283 _hist_meff_D->fill(m_eff_inc,weight); 00284 if(m_eff_inc>1700.) _count_D_tight ->fill(0.5,weight); 00285 } 00286 00287 // for rest of regions 6 jets pT> 40 needed 00288 if(recon_jets.size()<6 || recon_jets[5].momentum().perp()<40.) 00289 vetoEvent; 00290 00291 // region E 00292 m_eff_Nj += recon_jets[5].momentum().perp(); 00293 if( min_dPhi_3 > 0.4 && min_dPhi_All > 0.2 && eTmiss/m_eff_Nj > 0.15 ) { 00294 _hist_meff_E_tight->fill(m_eff_inc,weight); 00295 if( eTmiss/m_eff_Nj > 0.25 ) 00296 _hist_meff_E_medium->fill(m_eff_inc,weight); 00297 if( eTmiss/m_eff_Nj > 0.3 ) 00298 _hist_meff_E_loose->fill(m_eff_inc,weight); 00299 if(m_eff_inc>1400.) _count_E_tight ->fill(0.5,weight); 00300 if(m_eff_inc>1300.&& eTmiss/m_eff_Nj > 0.25 ) 00301 _count_E_medium->fill(0.5,weight); 00302 if(m_eff_inc>1000.&& eTmiss/m_eff_Nj > 0.3 ) 00303 _count_E_loose ->fill(0.5,weight); 00304 } 00305 } 00306 00307 00308 void finalize() { 00309 00310 double norm = crossSection()/femtobarn*5.8/sumOfWeights(); 00311 // these are number of events at 5.8fb^-1 per 100 GeV 00312 scale( _hist_meff_A_medium , 100. * norm ); 00313 scale( _hist_meff_A_tight , 100. * norm ); 00314 scale( _hist_meff_B_medium , 100. * norm ); 00315 scale( _hist_meff_B_tight , 100. * norm ); 00316 scale( _hist_meff_C_medium , 100. * norm ); 00317 scale( _hist_meff_C_tight , 100. * norm ); 00318 scale( _hist_meff_D , 100. * norm ); 00319 scale( _hist_meff_E_loose , 100. * norm ); 00320 scale( _hist_meff_E_medium , 100. * norm ); 00321 scale( _hist_meff_E_tight , 100. * norm ); 00322 // these are number of events at 5.8fb^-1 00323 scale(_count_A_tight ,norm); 00324 scale(_count_A_medium ,norm); 00325 scale(_count_A_loose ,norm); 00326 scale(_count_B_tight ,norm); 00327 scale(_count_B_medium ,norm); 00328 scale(_count_C_tight ,norm); 00329 scale(_count_C_medium ,norm); 00330 scale(_count_C_loose ,norm); 00331 scale(_count_D_tight ,norm); 00332 scale(_count_E_tight ,norm); 00333 scale(_count_E_medium ,norm); 00334 scale(_count_E_loose ,norm); 00335 } 00336 00337 //@} 00338 00339 private: 00340 00341 Histo1DPtr _count_A_tight; 00342 Histo1DPtr _count_A_medium; 00343 Histo1DPtr _count_A_loose; 00344 Histo1DPtr _count_B_tight; 00345 Histo1DPtr _count_B_medium; 00346 Histo1DPtr _count_C_tight; 00347 Histo1DPtr _count_C_medium; 00348 Histo1DPtr _count_C_loose; 00349 Histo1DPtr _count_D_tight; 00350 Histo1DPtr _count_E_tight; 00351 Histo1DPtr _count_E_medium; 00352 Histo1DPtr _count_E_loose; 00353 00354 Histo1DPtr _hist_meff_A_medium; 00355 Histo1DPtr _hist_meff_A_tight; 00356 Histo1DPtr _hist_meff_B_medium; 00357 Histo1DPtr _hist_meff_B_tight; 00358 Histo1DPtr _hist_meff_C_medium; 00359 Histo1DPtr _hist_meff_C_tight; 00360 Histo1DPtr _hist_meff_D; 00361 Histo1DPtr _hist_meff_E_loose; 00362 Histo1DPtr _hist_meff_E_medium; 00363 Histo1DPtr _hist_meff_E_tight; 00364 00365 }; 00366 00367 00368 // This global object acts as a hook for the plugin system 00369 DECLARE_RIVET_PLUGIN(ATLAS_2012_CONF_2012_109); 00370 00371 } Generated on Thu Feb 6 2014 17:38:41 for The Rivet MC analysis system by 1.7.6.1 |