ATLAS_2015_I1404878.cc
Go to the documentation of this file.
00001 // -*- C++ -*- 00002 #include "Rivet/Analysis.hh" 00003 #include "Rivet/Projections/FinalState.hh" 00004 #include "Rivet/Projections/VetoedFinalState.hh" 00005 #include "Rivet/Projections/IdentifiedFinalState.hh" 00006 #include "Rivet/Projections/PromptFinalState.hh" 00007 #include "Rivet/Projections/DressedLeptons.hh" 00008 #include "Rivet/Projections/FastJets.hh" 00009 #include "Rivet/Projections/VisibleFinalState.hh" 00010 00011 namespace Rivet { 00012 00013 class ATLAS_2015_I1404878 : public Analysis { 00014 public: 00015 00016 /// Constructor 00017 DEFAULT_RIVET_ANALYSIS_CTOR(ATLAS_2015_I1404878); 00018 00019 void init() { 00020 // Eta ranges 00021 Cut eta_full = (Cuts::abseta < 4.2) & (Cuts::pT >= 1.0*MeV); 00022 Cut lep_cuts = (Cuts::abseta < 2.5) && (Cuts::pT > 25*GeV); 00023 00024 // All final state particles 00025 FinalState fs(eta_full); 00026 00027 // Get photons to dress leptons 00028 IdentifiedFinalState photons(fs); 00029 photons.acceptIdPair(PID::PHOTON); 00030 00031 // Projection to find the electrons 00032 IdentifiedFinalState el_id(fs); 00033 el_id.acceptIdPair(PID::ELECTRON); 00034 00035 PromptFinalState electrons(el_id); 00036 electrons.acceptTauDecays(true); 00037 declare(electrons, "electrons"); 00038 00039 DressedLeptons dressedelectrons(photons, electrons, 0.1, lep_cuts, true, true); 00040 declare(dressedelectrons, "dressedelectrons"); 00041 00042 DressedLeptons ewdressedelectrons(photons, electrons, 0.1, eta_full, true, true); 00043 declare(ewdressedelectrons, "ewdressedelectrons"); 00044 00045 // Projection to find the muons 00046 IdentifiedFinalState mu_id(fs); 00047 mu_id.acceptIdPair(PID::MUON); 00048 00049 PromptFinalState muons(mu_id); 00050 muons.acceptTauDecays(true); 00051 declare(muons, "muons"); 00052 00053 DressedLeptons dressedmuons(photons, muons, 0.1, lep_cuts, true, true); 00054 declare(dressedmuons, "dressedmuons"); 00055 00056 DressedLeptons ewdressedmuons(photons, muons, 0.1, eta_full, true, true); 00057 declare(ewdressedmuons, "ewdressedmuons"); 00058 00059 // Projection to find neutrinos 00060 IdentifiedFinalState nu_id; 00061 nu_id.acceptNeutrinos(); 00062 PromptFinalState neutrinos(nu_id); 00063 neutrinos.acceptTauDecays(true); 00064 00065 // get MET from generic invisibles 00066 VetoedFinalState inv_fs(fs); 00067 inv_fs.addVetoOnThisFinalState(VisibleFinalState(fs)); 00068 declare(inv_fs, "InvisibleFS"); 00069 00070 // Jet clustering. 00071 VetoedFinalState vfs; 00072 vfs.addVetoOnThisFinalState(ewdressedelectrons); 00073 vfs.addVetoOnThisFinalState(ewdressedmuons); 00074 vfs.addVetoOnThisFinalState(neutrinos); 00075 FastJets jets(vfs, FastJets::ANTIKT, 0.4); 00076 jets.useInvisibles(true); 00077 declare(jets, "jets"); 00078 00079 // Histogram booking 00080 _h["massttbar"] = bookHisto1D( 1, 1, 1); 00081 _h["massttbar_norm"] = bookHisto1D( 2, 1, 1); 00082 _h["ptttbar"] = bookHisto1D( 3, 1, 1); 00083 _h["ptttbar_norm"] = bookHisto1D( 4, 1, 1); 00084 _h["absrapttbar"] = bookHisto1D( 5, 1, 1); 00085 _h["absrapttbar_norm"] = bookHisto1D( 6, 1, 1); 00086 _h["ptpseudotophadron"] = bookHisto1D( 7, 1, 1); 00087 _h["ptpseudotophadron_norm"] = bookHisto1D( 8, 1, 1); 00088 _h["absrappseudotophadron"] = bookHisto1D( 9, 1, 1); 00089 _h["absrappseudotophadron_norm"] = bookHisto1D(10, 1, 1); 00090 _h["absPout"] = bookHisto1D(11, 1, 1); 00091 _h["absPout_norm"] = bookHisto1D(12, 1, 1); 00092 _h["dPhittbar"] = bookHisto1D(13, 1, 1); 00093 _h["dPhittbar_norm"] = bookHisto1D(14, 1, 1); 00094 _h["HTttbar"] = bookHisto1D(15, 1, 1); 00095 _h["HTttbar_norm"] = bookHisto1D(16, 1, 1); 00096 _h["Yboost"] = bookHisto1D(17, 1, 1); 00097 _h["Yboost_norm"] = bookHisto1D(18, 1, 1); 00098 _h["chittbar"] = bookHisto1D(19, 1, 1); 00099 _h["chittbar_norm"] = bookHisto1D(20, 1, 1); 00100 _h["RWt"] = bookHisto1D(21, 1, 1); 00101 _h["RWt_norm"] = bookHisto1D(22, 1, 1); 00102 00103 } 00104 00105 void analyze(const Event& event) { 00106 00107 const double weight = event.weight(); 00108 00109 // Get the selected objects, using the projections. 00110 vector<DressedLepton> electrons = applyProjection<DressedLeptons>(event, "dressedelectrons").dressedLeptons(); 00111 vector<DressedLepton> muons = applyProjection<DressedLeptons>(event, "dressedmuons").dressedLeptons(); 00112 const Jets& jets = applyProjection<FastJets>(event, "jets").jetsByPt(Cuts::pT > 25*GeV && Cuts::abseta < 2.5); 00113 const FinalState& ifs = applyProjection<FinalState>(event, "InvisibleFS"); 00114 00115 // Calculate MET 00116 FourMomentum met; 00117 foreach (const Particle& p, ifs.particles()) met += p.momentum(); 00118 00119 // Count the number of b-tags 00120 Jets bjets, lightjets; 00121 foreach (const Jet& jet, jets){ 00122 bool b_tagged = jet.bTags(Cuts::pT > 5*GeV).size(); 00123 if ( b_tagged && bjets.size() < 2 ) bjets += jet; 00124 else lightjets += jet; 00125 } 00126 00127 bool single_electron = (electrons.size() == 1) && (muons.empty()); 00128 bool single_muon = (muons.size() == 1) && (electrons.empty()); 00129 00130 DressedLepton *lepton = NULL; 00131 if (single_electron) lepton = &electrons[0]; 00132 else if (single_muon) lepton = &muons[0]; 00133 00134 if(!single_electron && !single_muon) vetoEvent; 00135 if (jets.size() < 4 || bjets.size() < 2) vetoEvent; 00136 00137 FourMomentum pbjet1; //Momentum of bjet1 00138 FourMomentum pbjet2; //Momentum of bjet2 00139 if ( deltaR(bjets[0], *lepton) <= deltaR(bjets[1], *lepton) ) { 00140 pbjet1 = bjets[0].momentum(); 00141 pbjet2 = bjets[1].momentum(); 00142 } else { 00143 pbjet1 = bjets[1].momentum(); 00144 pbjet2 = bjets[0].momentum(); 00145 } 00146 00147 00148 double bestWmass = 1000.0*TeV; 00149 double mWPDG = 80.399*GeV; 00150 int Wj1index = -1, Wj2index = -1; 00151 for (unsigned int i = 0; i < (lightjets.size() - 1); ++i) { 00152 for (unsigned int j = i + 1; j < lightjets.size(); ++j) { 00153 double wmass = (lightjets[i].momentum() + lightjets[j].momentum()).mass(); 00154 if (fabs(wmass - mWPDG) < fabs(bestWmass - mWPDG)) { 00155 bestWmass = wmass; 00156 Wj1index = i; 00157 Wj2index = j; 00158 } 00159 } 00160 } 00161 00162 // compute hadronic W boson 00163 FourMomentum pWhadron = lightjets[Wj1index].momentum() + lightjets[Wj2index].momentum(); 00164 double pz = computeneutrinoz(lepton->momentum(), met); 00165 FourMomentum ppseudoneutrino( sqrt(sqr(met.px()) + sqr(met.py()) + sqr(pz)), met.px(), met.py(), pz); 00166 00167 00168 //compute leptonic, hadronic, combined pseudo-top 00169 FourMomentum ppseudotoplepton = lepton->momentum() + ppseudoneutrino + pbjet1; 00170 FourMomentum ppseudotophadron = pbjet2 + pWhadron; 00171 FourMomentum pttbar = ppseudotoplepton + ppseudotophadron; 00172 00173 Vector3 z_versor(0,0,1); 00174 Vector3 vpseudotophadron = ppseudotophadron.vector3(); 00175 Vector3 vpseudotoplepton = ppseudotoplepton.vector3(); 00176 // Observables 00177 double ystar = 0.5 * deltaRap(ppseudotophadron, ppseudotoplepton); 00178 double chi_ttbar = exp(2 * fabs(ystar)); 00179 double deltaPhi_ttbar = deltaPhi(ppseudotoplepton,ppseudotophadron); 00180 double HT_ttbar = ppseudotophadron.pt() + ppseudotoplepton.pt(); 00181 double Yboost = 0.5 * fabs(ppseudotophadron.rapidity() + ppseudotoplepton.rapidity()); 00182 double R_Wt = pWhadron.pt() / ppseudotophadron.pt(); 00183 double absPout = fabs(vpseudotophadron.dot((vpseudotoplepton.cross(z_versor))/(vpseudotoplepton.cross(z_versor).mod()))); 00184 00185 // absolute cross sections 00186 _h["ptpseudotophadron"]->fill( ppseudotophadron.pt(), weight); //pT of pseudo top hadron 00187 _h["ptttbar"]->fill( pttbar.pt(), weight); //fill pT of ttbar in combined channel 00188 _h["absrappseudotophadron"]->fill(ppseudotophadron.absrap(), weight); 00189 _h["absrapttbar"]->fill( pttbar.absrap(), weight); 00190 _h["massttbar"]->fill( pttbar.mass(), weight); 00191 _h["absPout"]->fill( absPout, weight); 00192 _h["chittbar"]->fill( chi_ttbar, weight); 00193 _h["dPhittbar"]->fill( deltaPhi_ttbar, weight); 00194 _h["HTttbar"]->fill( HT_ttbar, weight); 00195 _h["Yboost"]->fill( Yboost, weight); 00196 _h["RWt"]->fill( R_Wt, weight); 00197 // normalised cross sections 00198 _h["ptpseudotophadron_norm"]->fill( ppseudotophadron.pt(), weight); //pT of pseudo top hadron 00199 _h["ptttbar_norm"]->fill( pttbar.pt(), weight); //fill pT of ttbar in combined channel 00200 _h["absrappseudotophadron_norm"]->fill(ppseudotophadron.absrap(), weight); 00201 _h["absrapttbar_norm"]->fill( pttbar.absrap(), weight); 00202 _h["massttbar_norm"]->fill( pttbar.mass(), weight); 00203 _h["absPout_norm"]->fill( absPout, weight); 00204 _h["chittbar_norm"]->fill( chi_ttbar, weight); 00205 _h["dPhittbar_norm"]->fill( deltaPhi_ttbar, weight); 00206 _h["HTttbar_norm"]->fill( HT_ttbar, weight); 00207 _h["Yboost_norm"]->fill( Yboost, weight); 00208 _h["RWt_norm"]->fill( R_Wt, weight); 00209 00210 } 00211 00212 void finalize() { 00213 // Normalize to cross-section 00214 const double sf = (crossSection() / sumOfWeights()); 00215 for (map<string, Histo1DPtr>::iterator hit = _h.begin(); hit != _h.end(); ++hit) { 00216 scale(hit->second, sf); 00217 if (hit->first.find("_norm") != string::npos) normalize(hit->second); 00218 } 00219 } 00220 00221 private: 00222 00223 00224 double computeneutrinoz(const FourMomentum& lepton, FourMomentum& met) const { 00225 //computing z component of neutrino momentum given lepton and met 00226 double pzneutrino; 00227 double m_W = 80.399; // in GeV, given in the paper 00228 double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py()); 00229 double a = sqr ( lepton.E() )- sqr ( lepton.pz() ); 00230 double b = -2*k*lepton.pz(); 00231 double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k ); 00232 double discriminant = sqr(b) - 4 * a * c; 00233 double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns 00234 if (discriminant < 0) pzneutrino = - b / (2 * a); //if the discriminant is negative 00235 else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value 00236 double absquad[2]; 00237 for (int n=0; n<2; ++n) absquad[n] = fabs(quad[n]); 00238 if (absquad[0] < absquad[1]) pzneutrino = quad[0]; 00239 else pzneutrino = quad[1]; 00240 } 00241 00242 return pzneutrino; 00243 } 00244 00245 double _mT(const FourMomentum &l, FourMomentum &nu) const { 00246 return sqrt( 2 * l.pT() * nu.pT() * (1 - cos(deltaPhi(l, nu))) ); 00247 } 00248 00249 /// @name Objects that are used by the event selection decisions 00250 map<string, Histo1DPtr> _h; 00251 }; 00252 00253 // The hook for the plugin system 00254 DECLARE_RIVET_PLUGIN(ATLAS_2015_I1404878); 00255 00256 } Generated on Tue Dec 13 2016 16:32:35 for The Rivet MC analysis system by 
1.7.6.1
|
