ATLAS_2012_I1203852.cc
Go to the documentation of this file.
00001 // -*- C++ -*- 00002 #include "Rivet/Analysis.hh" 00003 #include "Rivet/Projections/FastJets.hh" 00004 #include "Rivet/Projections/FinalState.hh" 00005 #include "Rivet/Projections/IdentifiedFinalState.hh" 00006 #include "Rivet/Projections/WFinder.hh" 00007 #include "Rivet/Projections/LeadingParticlesFinalState.hh" 00008 #include "Rivet/Projections/UnstableFinalState.hh" 00009 #include "Rivet/Projections/VetoedFinalState.hh" 00010 #include "Rivet/Projections/DressedLeptons.hh" 00011 #include "Rivet/Projections/MergedFinalState.hh" 00012 #include "Rivet/Projections/MissingMomentum.hh" 00013 #include "Rivet/Projections/InvMassFinalState.hh" 00014 00015 #define ZMASS 91.1876 // GeV 00016 00017 namespace Rivet { 00018 00019 using namespace Cuts; 00020 00021 00022 /// Generic Z candidate 00023 struct Zstate : public ParticlePair { 00024 Zstate() { } 00025 Zstate(ParticlePair _particlepair) : ParticlePair(_particlepair) { } 00026 FourMomentum mom() const { return first.momentum() + second.momentum(); } 00027 operator FourMomentum() const { return mom(); } 00028 static bool cmppT(const Zstate& lx, const Zstate& rx) { return lx.mom().pT() < rx.mom().pT(); } 00029 }; 00030 00031 00032 /// 4l to ZZ assignment -- algorithm 00033 void identifyZstates(Zstate & Z1, Zstate &Z2, const Particles & leptons_sel4l) { 00034 00035 ///////////////////////////////////////////////////////////////////////////// 00036 /// ZZ->4l pairing 00037 /// - Exactly two same flavour opposite charged leptons 00038 /// - Ambiguities in pairing are resolved by choosing the combination 00039 /// that results in the smaller value of the sum |mll - mZ| for the two pairs 00040 ///////////////////////////////////////////////////////////////////////////// 00041 00042 Particles part_pos_el, part_neg_el, part_pos_mu, part_neg_mu; 00043 foreach (const Particle& l , leptons_sel4l) { 00044 if (l.abspid() == PID::ELECTRON) { 00045 if (l.pid() < 0) part_neg_el.push_back(l); 00046 if (l.pid() > 0) part_pos_el.push_back(l); 00047 } 00048 else if (l.abspid() == PID::MUON) { 00049 if (l.pid() < 0) part_neg_mu.push_back(l); 00050 if (l.pid() > 0) part_pos_mu.push_back(l); 00051 } 00052 } 00053 00054 // ee/mm channel 00055 if ( part_neg_el.size() == 2 || part_neg_mu.size() == 2) { 00056 00057 Zstate Zcand_1, Zcand_2, Zcand_3, Zcand_4; 00058 if (part_neg_el.size() == 2) { // ee 00059 Zcand_1 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[0] ) ); 00060 Zcand_2 = Zstate( ParticlePair( part_neg_el[0], part_pos_el[1] ) ); 00061 Zcand_3 = Zstate( ParticlePair( part_neg_el[1], part_pos_el[0] ) ); 00062 Zcand_4 = Zstate( ParticlePair( part_neg_el[1], part_pos_el[1] ) ); 00063 } else { // mumu 00064 Zcand_1 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[0] ) ); 00065 Zcand_2 = Zstate( ParticlePair( part_neg_mu[0], part_pos_mu[1] ) ); 00066 Zcand_3 = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[0] ) ); 00067 Zcand_4 = Zstate( ParticlePair( part_neg_mu[1], part_pos_mu[1] ) ); 00068 } 00069 00070 // We can have the following pairs: (Z1 + Z4) || (Z2 + Z3) 00071 double minValue_1, minValue_2; 00072 minValue_1 = fabs( Zcand_1.mom().mass() - ZMASS ) + fabs( Zcand_4.mom().mass() - ZMASS); 00073 minValue_2 = fabs( Zcand_2.mom().mass() - ZMASS ) + fabs( Zcand_3.mom().mass() - ZMASS); 00074 if (minValue_1 < minValue_2 ) { 00075 Z1 = Zcand_1; 00076 Z2 = Zcand_4; 00077 } else { 00078 Z1 = Zcand_2; 00079 Z2 = Zcand_3; 00080 } 00081 00082 // emu channel 00083 } else if (part_neg_mu.size() == 1 && part_neg_el.size() == 1) { 00084 Z1 = Zstate ( ParticlePair (part_neg_mu[0], part_pos_mu[0] ) ); 00085 Z2 = Zstate ( ParticlePair (part_neg_el[0], part_pos_el[0] ) ); 00086 } 00087 00088 } 00089 00090 00091 /// @name ZZ analysis 00092 class ATLAS_2012_I1203852 : public Analysis { 00093 public: 00094 00095 /// Default constructor 00096 ATLAS_2012_I1203852() 00097 : Analysis("ATLAS_2012_I1203852") 00098 { } 00099 00100 void init() { 00101 FinalState fs(-5.0, 5.0, 0.0*GeV); 00102 00103 // Final states to form Z bosons 00104 vids.push_back(make_pair(PID::ELECTRON, PID::POSITRON)); 00105 vids.push_back(make_pair(PID::MUON, PID::ANTIMUON)); 00106 00107 vnuids.push_back(make_pair(PID::NU_E, PID::NU_EBAR) ); 00108 vnuids.push_back(make_pair(PID::NU_MU, PID::NU_MUBAR) ); 00109 vnuids.push_back(make_pair(PID::NU_TAU, PID::NU_TAUBAR) ); 00110 00111 IdentifiedFinalState Photon(fs); 00112 Photon.acceptIdPair(PID::PHOTON); 00113 00114 IdentifiedFinalState bare_EL(fs); 00115 bare_EL.acceptIdPair(PID::ELECTRON); 00116 00117 IdentifiedFinalState bare_MU(fs); 00118 bare_MU.acceptIdPair(PID::MUON); 00119 00120 00121 // Selection 1: ZZ-> llll selection 00122 Cut etaranges_lep = etaIn(-3.16, 3.16) & (pT > 7*GeV); 00123 00124 DressedLeptons electron_sel4l(Photon, bare_EL, 0.1, true/*Do clustering*/, 00125 etaranges_lep, false/*false = ignore photons from hadron or tau*/); 00126 addProjection(electron_sel4l, "ELECTRON_sel4l"); 00127 DressedLeptons muon_sel4l(Photon, bare_MU, 0.1,true, etaranges_lep, false); 00128 addProjection(muon_sel4l, "MUON_sel4l"); 00129 00130 00131 // Selection 2: ZZ-> llnunu selection 00132 Cut etaranges_lep2 = etaIn(-2.5, 2.5) & (pT > 10*GeV); 00133 00134 DressedLeptons electron_sel2l2nu(Photon, bare_EL, 0.1,true, etaranges_lep2, false); 00135 addProjection(electron_sel2l2nu, "ELECTRON_sel2l2nu"); 00136 DressedLeptons muon_sel2l2nu(Photon, bare_MU, 0.1,true, etaranges_lep2, false); 00137 addProjection(muon_sel2l2nu, "MUON_sel2l2nu"); 00138 00139 00140 /// Get all neutrinos. These will not be used to form jets. 00141 /// We'll use the highest 2 pT neutrinos to calculate the MET 00142 IdentifiedFinalState neutrino_fs(-4.5, 4.5, 0.0*GeV); 00143 neutrino_fs.acceptNeutrinos(); 00144 addProjection(neutrino_fs, "NEUTRINO_FS"); 00145 00146 VetoedFinalState jetinput; 00147 jetinput.addVetoOnThisFinalState(bare_MU); 00148 jetinput.addVetoOnThisFinalState(neutrino_fs); 00149 00150 FastJets jetpro(fs, FastJets::ANTIKT, 0.4); 00151 addProjection(jetpro, "jet"); 00152 00153 // Both ZZ on-shell histos 00154 _h_ZZ_xsect = bookHisto1D(1, 1, 1); 00155 _h_ZZ_ZpT = bookHisto1D(3, 1, 1); 00156 _h_ZZ_phill = bookHisto1D(5, 1, 1); 00157 _h_ZZ_mZZ = bookHisto1D(7, 1, 1); 00158 00159 // One Z off-shell (ZZstar) histos 00160 _h_ZZs_xsect = bookHisto1D(1, 1, 2); 00161 00162 // ZZ -> llnunu histos 00163 _h_ZZnunu_xsect = bookHisto1D(1, 1, 3); 00164 _h_ZZnunu_ZpT = bookHisto1D(4, 1, 1); 00165 _h_ZZnunu_phill = bookHisto1D(6, 1, 1); 00166 _h_ZZnunu_mZZ = bookHisto1D(8, 1, 1); 00167 } 00168 00169 00170 /// Do the analysis 00171 void analyze(const Event& e) { 00172 const double weight = e.weight(); 00173 00174 //////////////////////////////////////////////////////////////////// 00175 // preselection of leptons for ZZ-> llll final state 00176 //////////////////////////////////////////////////////////////////// 00177 00178 Particles leptons_sel4l; 00179 00180 const vector<DressedLepton>& mu_sel4l = applyProjection<DressedLeptons>(e, "MUON_sel4l").clusteredLeptons(); 00181 const vector<DressedLepton>& el_sel4l = applyProjection<DressedLeptons>(e, "ELECTRON_sel4l").clusteredLeptons(); 00182 00183 vector<DressedLepton> leptonsFS_sel4l; 00184 leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), mu_sel4l.begin(), mu_sel4l.end() ); 00185 leptonsFS_sel4l.insert( leptonsFS_sel4l.end(), el_sel4l.begin(), el_sel4l.end() ); 00186 00187 //////////////////////////////////////////////////////////////////// 00188 // OVERLAP removal dR(l,l)>0.2 00189 //////////////////////////////////////////////////////////////////// 00190 foreach ( const DressedLepton& l1, leptonsFS_sel4l) { 00191 bool isolated = true; 00192 foreach (DressedLepton& l2, leptonsFS_sel4l) { 00193 const double dR = deltaR(l1, l2); 00194 if (dR < 0.2 && l1 != l2) { isolated = false; break; } 00195 } 00196 if (isolated) leptons_sel4l.push_back(l1); 00197 } 00198 00199 ////////////////////////////////////////////////////////////////// 00200 // Exactly two opposite charged leptons 00201 ////////////////////////////////////////////////////////////////// 00202 00203 // calculate total 'flavour' charge 00204 double totalcharge = 0; 00205 foreach (Particle& l, leptons_sel4l) totalcharge += l.pid(); 00206 00207 // Analyze 4 lepton events 00208 if (leptons_sel4l.size() == 4 && totalcharge == 0 ) { 00209 Zstate Z1, Z2; 00210 00211 // Identifies Z states from 4 lepton pairs 00212 identifyZstates(Z1, Z2,leptons_sel4l); 00213 00214 //////////////////////////////////////////////////////////////////////////// 00215 // Z MASS WINDOW 00216 // -ZZ: for both Z: 66<mZ<116 GeV 00217 // -ZZ*: one Z on-shell: 66<mZ<116 GeV, one Z off-shell: mZ>20 GeV 00218 /////////////////////////////////////////////////////////////////////////// 00219 00220 Zstate leadPtZ = std::max(Z1, Z2, Zstate::cmppT); 00221 00222 double mZ1 = Z1.mom().mass(); 00223 double mZ2 = Z2.mom().mass(); 00224 double ZpT = leadPtZ.mom().pT(); 00225 double phill = fabs(deltaPhi(leadPtZ.first, leadPtZ.second)); 00226 if (phill > M_PI) phill = 2*M_PI-phill; 00227 double mZZ = (Z1.mom() + Z2.mom()).mass(); 00228 00229 if (mZ1 > 20*GeV && mZ2 > 20*GeV) { 00230 // ZZ* selection 00231 if (inRange(mZ1, 66*GeV, 116*GeV) || inRange(mZ2, 66*GeV, 116*GeV)) { 00232 _h_ZZs_xsect -> fill(sqrtS()*GeV, weight); 00233 } 00234 00235 // ZZ selection 00236 if (inRange(mZ1, 66*GeV, 116*GeV) && inRange(mZ2, 66*GeV, 116*GeV)) { 00237 _h_ZZ_xsect -> fill(sqrtS()*GeV, weight); 00238 _h_ZZ_ZpT -> fill(ZpT , weight); 00239 _h_ZZ_phill -> fill(phill , weight); 00240 _h_ZZ_mZZ -> fill(mZZ , weight); 00241 } 00242 } 00243 } 00244 00245 //////////////////////////////////////////////////////////////////// 00246 /// preselection of leptons for ZZ-> llnunu final state 00247 //////////////////////////////////////////////////////////////////// 00248 00249 Particles leptons_sel2l2nu; // output 00250 const vector<DressedLepton>& mu_sel2l2nu = applyProjection<DressedLeptons>(e, "MUON_sel2l2nu").clusteredLeptons(); 00251 const vector<DressedLepton>& el_sel2l2nu = applyProjection<DressedLeptons>(e, "ELECTRON_sel2l2nu").clusteredLeptons(); 00252 00253 vector<DressedLepton> leptonsFS_sel2l2nu; 00254 leptonsFS_sel2l2nu.insert( leptonsFS_sel2l2nu.end(), mu_sel2l2nu.begin(), mu_sel2l2nu.end() ); 00255 leptonsFS_sel2l2nu.insert( leptonsFS_sel2l2nu.end(), el_sel2l2nu.begin(), el_sel2l2nu.end() ); 00256 00257 // Lepton preselection for ZZ-> llnunu 00258 if ((mu_sel2l2nu.empty() || el_sel2l2nu.empty()) // cannot have opposite flavour 00259 && (leptonsFS_sel2l2nu.size() == 2) // exactly two leptons 00260 && (leptonsFS_sel2l2nu[0].charge() * leptonsFS_sel2l2nu[1].charge() < 1 ) // opposite charge 00261 && (deltaR(leptonsFS_sel2l2nu[0], leptonsFS_sel2l2nu[1]) > 0.3) // overlap removal 00262 && (leptonsFS_sel2l2nu[0].pT() > 20*GeV && leptonsFS_sel2l2nu[1].pT() > 20*GeV)) { // trigger requirement 00263 leptons_sel2l2nu.insert(leptons_sel2l2nu.end(), leptonsFS_sel2l2nu.begin(), leptonsFS_sel2l2nu.end()); 00264 } 00265 if (leptons_sel2l2nu.empty()) vetoEvent; // no further analysis, fine to veto 00266 00267 Particles leptons_sel2l2nu_jetveto; 00268 foreach (const DressedLepton& l, mu_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.constituentLepton()); 00269 foreach (const DressedLepton& l, el_sel2l2nu) leptons_sel2l2nu_jetveto.push_back(l.constituentLepton()); 00270 double ptll = (leptons_sel2l2nu[0].momentum() + leptons_sel2l2nu[1].momentum()).pT(); 00271 00272 // Find Z1-> ll 00273 FinalState fs2(-3.2, 3.2); 00274 InvMassFinalState imfs(fs2, vids, 20*GeV, sqrtS()); 00275 imfs.calc(leptons_sel2l2nu); 00276 if (imfs.particlePairs().size() != 1) vetoEvent; 00277 const ParticlePair& Z1constituents = imfs.particlePairs()[0]; 00278 FourMomentum Z1 = Z1constituents.first.momentum() + Z1constituents.second.momentum(); 00279 00280 // Find Z2-> nunu 00281 const Particles& neutrinoFS = applyProjection<IdentifiedFinalState>(e, "NEUTRINO_FS").particlesByPt(); 00282 FinalState fs3(-4.5, 4.5); 00283 InvMassFinalState imfs_Znunu(fs3, vnuids, 20*GeV, sqrtS(), ZMASS); 00284 imfs_Znunu.calc(neutrinoFS); 00285 if (imfs_Znunu.particlePairs().size() == 0 ) vetoEvent; 00286 00287 // Z to neutrinos 00288 FourMomentum Z2 = imfs_Znunu.particlePairs()[0].first.momentum() + imfs_Znunu.particlePairs()[0].second.momentum(); 00289 double met_Znunu = Z2.pT(); 00290 00291 // mTZZ 00292 const double mT2_1st_term = add_quad(ZMASS, ptll) + add_quad(ZMASS, met_Znunu); 00293 const double mT2_2nd_term = Z1.px() + Z2.px(); 00294 const double mT2_3rd_term = Z1.py() + Z2.py(); 00295 const double mTZZ = sqrt(sqr(mT2_1st_term) - sqr(mT2_2nd_term) - sqr(mT2_3rd_term)); 00296 00297 if (!inRange(Z2.mass(), 66*GeV, 116*GeV)) vetoEvent; 00298 if (!inRange(Z1.mass(), 76*GeV, 106*GeV)) vetoEvent; 00299 00300 ///////////////////////////////////////////////////////////// 00301 // AXIAL MET < 75 GeV 00302 //////////////////////////////////////////////////////////// 00303 00304 double dPhiZ1Z2 = fabs(deltaPhi(Z1, Z2)); 00305 if (dPhiZ1Z2 > M_PI) dPhiZ1Z2 = 2*M_PI - dPhiZ1Z2; 00306 const double axialEtmiss = -Z2.pT()*cos(dPhiZ1Z2); 00307 if (axialEtmiss < 75*GeV) vetoEvent; 00308 00309 const double ZpT = Z1.pT(); 00310 double phill = fabs(deltaPhi(Z1constituents.first, Z1constituents.second)); 00311 if (phill > M_PI) phill = 2*M_PI - phill; 00312 00313 00314 //////////////////////////////////////////////////////////////////////////// 00315 // JETS 00316 // -"j": found by "jetpro" projection && pT() > 25 GeV && |eta| < 4.5 00317 // -"goodjets": "j" && dR(electron/muon,jet) > 0.3 00318 // 00319 // JETVETO: veto all events with at least one good jet 00320 /////////////////////////////////////////////////////////////////////////// 00321 vector<Jet> good_jets; 00322 foreach (const Jet& j, applyProjection<FastJets>(e, "jet").jetsByPt(25)) { 00323 if (j.abseta() > 4.5) continue; 00324 bool isLepton = 0; 00325 foreach (const Particle& l, leptons_sel2l2nu_jetveto) { 00326 const double dR = deltaR(l.momentum(), j.momentum()); 00327 if (dR < 0.3) { isLepton = true; break; } 00328 } 00329 if (!isLepton) good_jets.push_back(j); 00330 } 00331 size_t n_sel_jets = good_jets.size(); 00332 if (n_sel_jets != 0) vetoEvent; 00333 00334 00335 ///////////////////////////////////////////////////////////// 00336 // Fractional MET and lepton pair difference: "RatioMet"< 0.4 00337 //////////////////////////////////////////////////////////// 00338 double ratioMet = fabs(Z2.pT() - Z1.pT()) / Z1.pT(); 00339 if (ratioMet > 0.4 ) vetoEvent; 00340 00341 00342 // End of ZZllnunu selection: now fill histograms 00343 _h_ZZnunu_xsect->fill(sqrtS()/GeV, weight); 00344 _h_ZZnunu_ZpT ->fill(ZpT, weight); 00345 _h_ZZnunu_phill->fill(phill, weight); 00346 _h_ZZnunu_mZZ ->fill(mTZZ, weight); 00347 00348 } 00349 00350 00351 /// Finalize 00352 void finalize() { 00353 const double norm = crossSection()/sumOfWeights()/femtobarn; 00354 00355 scale(_h_ZZ_xsect, norm); 00356 normalize(_h_ZZ_ZpT); 00357 normalize(_h_ZZ_phill); 00358 normalize(_h_ZZ_mZZ); 00359 00360 scale(_h_ZZs_xsect, norm); 00361 scale(_h_ZZnunu_xsect, norm); 00362 normalize(_h_ZZnunu_ZpT); 00363 normalize(_h_ZZnunu_phill); 00364 normalize(_h_ZZnunu_mZZ); 00365 } 00366 00367 00368 private: 00369 00370 Histo1DPtr _h_ZZ_xsect, _h_ZZ_ZpT, _h_ZZ_phill, _h_ZZ_mZZ; 00371 Histo1DPtr _h_ZZs_xsect; 00372 Histo1DPtr _h_ZZnunu_xsect, _h_ZZnunu_ZpT, _h_ZZnunu_phill, _h_ZZnunu_mZZ; 00373 vector< pair<PdgId,PdgId> > vids, vnuids; 00374 00375 }; 00376 00377 00378 // The hook for the plugin system 00379 DECLARE_RIVET_PLUGIN(ATLAS_2012_I1203852); 00380 00381 } Generated on Tue Sep 30 2014 19:45:42 for The Rivet MC analysis system by 
1.7.6.1
|
