// -*- C++ -*- // ATLAS W+c analysis ////////////////////////////////////////////////////////////////////////// /* Description of rivet analysis ATLAS_2014_I1282447 W+c production This rivet routine implements the ATLAS W+c analysis. Apart from those histograms, described and published on HEP Data, here are some helper histograms defined, these are: [note that the division histograms might not work, if you are running +/- channels separately] d02-x01-y01_plus/minus d02-x01-y02_plus/minus --> to calculate the W+/W- charge ratio (inclusive cross section) d08-x01-y01_plus/minus --> to calculate the W+/W- charge ratio (jet mutliplicity distribution) the following histograms only work, when a second RIVET run is done on inclusive W production. d05-x01-y01 --> to calcuate the relative WD / W inclusive ratio for the inclusive cross sections d06-x01-y01_winc --> to calculate the W+DpT / Winclusive ratio for the pT dependent cross sections The histograms can be divided using the following code, python divideWCharm.py import yoda hists_wc = yoda.read("Rivet_Wc.yoda") hists_winc = yoda.read("Rivet_Winc.yoda") // inclusive cross section h_winc= hists_winc["/ATLAS_2014_I1282447/d05-x01-y01"] h_d = hists_wc["/ATLAS_2014_I1282447/d01-x01-y02"] h_dstar= hists_wc["/ATLAS_2014_I1282447/d01-x01-y03"] ratio_wd = h_d.divide(h_winc) ratio_wd.path = "/ATLAS_2014_I1282447/d05-x01-y01" ratio_wdstar = h_d.divide(h_winc) ratio_wdstar.path = "/ATLAS_2014_I1282447/d05-x01-y02" // pT differential h_winc plus= hists_wc["/ATLAS_2014_I1282447/d06-x01-y01_winc"] h_winc minus= hists_wc["/ATLAS_2014_I1282447/d06-x01-y02_winc"] h_wd_plus = hists_winc["/ATLAS_2014_I1282447/d06-x01-y01"] h_wd_minus = hists_winc["/ATLAS_2014_I1282447/d06-x01-y02"] h_wdstar_plus = hists_winc["/ATLAS_2014_I1282447/d06-x01-y03"] h_wdstar_minus = hists_winc["/ATLAS_2014_I1282447/d06-x01-y04"] ratio_wd_plus = h_wd_plus.divide(h_winc plus) ratio_wd_plus.path = "/ATLAS_2014_I1282447/d06-x01-y01" ratio_wd_minus = h_wd_plus.divide(h_winc minus) ratio_wd_minus.path = "/ATLAS_2014_I1282447/d06-x01-y02" ratio_wdstar_plus = h_wdstar_plus.divide(h_winc plus) ratio_wdstar_plus.path = "/ATLAS_2014_I1282447/d06-x01-y03" ratio_wdstar_minus = h_wdstar_plus.divide(h_winc minus) ratio_wdstar_minus.path = "/ATLAS_2014_I1282447/d06-x01-y04" */ ////////////////////////////////////////////////////////////////////////// #include "Rivet/Analysis.hh" #include "Rivet/Projections/UnstableFinalState.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/WFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/ChargedFinalState.hh" #include "Rivet/Projections/VetoedFinalState.hh" #include "HepMC/GenEvent.h" #include "HepMC/GenParticle.h" #include "HepMC/GenVertex.h" #include "HepMC/SimpleVector.h" /// @todo Include more projections as required, e.g. ChargedFinalState, FastJets, ZFinder... namespace Rivet { namespace { /// @todo Correct this in Rivet inline double calcMT(const FourMomentum& a, const FourMomentum& b) { return sqrt(2.0 * a.pT() * b.pT() * (1.0 - cos(a.phi() - b.phi())) ); } } class ATLAS_2014_I1282447 : public Analysis { public: /// Constructor ATLAS_2014_I1282447() : Analysis("ATLAS_2014_I1282447") { } public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { /// @todo Initialise and register projections here UnstableFinalState fs; WFinder wfinder_born_el(fs, -2.5, 2.5, 20*GeV, PID::ELECTRON, 25*GeV, 8000*GeV, 15*GeV, 0.3, WFinder::CLUSTERALL , WFinder::TRACK ); addProjection(wfinder_born_el, "WFinder_born_el"); WFinder wfinder_born_mu(fs, -2.5, 2.5, 20*GeV, PID::MUON, 25*GeV, 8000*GeV, 15*GeV, 0.1, WFinder::CLUSTERALL,WFinder::TRACK ); addProjection(wfinder_born_mu, "WFinder_born_mu"); // all hadrons that could be coming from a charm decay -- // -- for safety, use region -3.5 - 3.5 addProjection(UnstableFinalState(-3.5, 3.5, 0*GeV), "hadrons"); // Input for the jets: no neutrinos, no muons, and no electron which passed the electron cuts // also: NO electron, muon or tau (needed due to ATLAS jet truth reconstruction feature) VetoedFinalState veto; veto.addVetoOnThisFinalState(wfinder_born_el); veto.addVetoPairId(PID::ELECTRON); veto.addVetoPairId(PID::MUON); veto.addVetoPairId(PID::TAU); FastJets jets(veto, FastJets::ANTIKT, 0.4); addProjection(jets, "jets"); // Book histograms MSG_INFO("Booking histograms"); // charge separated integrated cross sections MSG_INFO("Booking charged cross sections"); _hist_wcjet_charge = bookHisto1D("d01-x01-y01"); _hist_wd_charge = bookHisto1D("d01-x01-y02"); _hist_wdstar_charge = bookHisto1D("d01-x01-y03"); // charge integrated total cross sections _hist_wcjet_ratio= bookScatter2D("d02-x01-y01"); _hist_wd_ratio = bookScatter2D("d02-x01-y02"); _hist_wcjet_minus= bookHisto1D("d02-x01-y01_minus"); _hist_wd_minus = bookHisto1D("d02-x01-y02_minus"); _hist_wcjet_plus= bookHisto1D("d02-x01-y01_plus"); _hist_wd_plus = bookHisto1D("d02-x01-y02_plus"); MSG_INFO("Booking eta lep wcjet"); // eta distributions _hist_wplus_wcjet_eta_lep = bookHisto1D("d03-x01-y01"); _hist_wminus_wcjet_eta_lep = bookHisto1D("d03-x01-y02"); MSG_INFO("Booking eta lep d"); _hist_wplus_wdminus_eta_lep = bookHisto1D("d04-x01-y01"); _hist_wminus_wdplus_eta_lep = bookHisto1D("d04-x01-y02"); MSG_INFO("Booking eta lep dstar"); _hist_wplus_wdstar_eta_lep = bookHisto1D("d04-x01-y03"); _hist_wminus_wdstar_eta_lep = bookHisto1D("d04-x01-y04"); MSG_INFO("Booking ratio of cross section (WD over W inclusive)"); // postprocess!!!!! // ratio of cross section (WD over W inclusive) // postprocess!!!!! _hist_w_inc = bookHisto1D("d05-x01-y01"); MSG_INFO("Booking ratio of cross section (WD over W inclusive -- function of pT of D meson) "); // postprocess!!!!! // ratio of cross section (WD over W inclsive -- function of pT of D meson) // postprocess!!!!! _hist_wplusd_wplusinc_pt = bookHisto1D("d06-x01-y01"); _hist_wminusd_wminusinc_pt = bookHisto1D("d06-x01-y02"); _hist_wplusdstar_wplusinc_pt = bookHisto1D("d06-x01-y03"); _hist_wminusdstar_wminusinc_pt = bookHisto1D("d06-x01-y04"); _hist_wplus_winc = bookHisto1D("d06-x01-y01_winc"); _hist_wminus_winc = bookHisto1D("d06-x01-y02_winc"); MSG_INFO("Booking _hist_wcjet_jets "); // jet multiplicity of charge integrated W+cjet cross section (+0 or +1 jet in addition to the charm jet) _hist_wcjet_jets = bookHisto1D("d07-x01-y01"); MSG_INFO("Booking _hist_wcjet_jets_ratio "); // jet multiplicity of W+cjet cross section ratio (+0 or +1 jet in addition to the charm jet) _hist_wcjet_jets_ratio = bookScatter2D("d08-x01-y01"); _hist_wcjet_jets_plus = bookHisto1D("d08-x01-y01_minus"); _hist_wcjet_jets_minus = bookHisto1D("d08-x01-y01_plus"); } /// Perform the per-event analysis void analyze(const Event& event) { events++; const double weight = event.weight(); double charge_weight = 0; // account for OS/SS events int lepton_charge=0; double lepton_eta=0; double meson_pt=0; /// Find leptons const WFinder& wfinder_born_el = applyProjection(event, "WFinder_born_el"); const WFinder& wfinder_born_mu = applyProjection(event, "WFinder_born_mu"); if (wfinder_born_el.empty() && wfinder_born_mu.empty() ) { MSG_DEBUG("No W bosons found"); vetoEvent; } bool keepevent=false; //check electrons if (!wfinder_born_el.particles().empty()) { const FourMomentum& nu = wfinder_born_el.constituentNeutrinos()[0].momentum(); if (calcMT(nu, wfinder_born_el.constituentLeptons()[0].momentum()) > 40*GeV && nu.pT()>25*GeV) { keepevent=true; lepton_charge=wfinder_born_el.constituentLeptons()[0].charge(); lepton_eta=fabs(wfinder_born_el.constituentLeptons()[0].pseudorapidity()); } } //check muons if (!wfinder_born_mu.particles().empty()) { const FourMomentum& nu = wfinder_born_mu.constituentNeutrinos()[0].momentum(); if ( calcMT(nu, wfinder_born_mu.constituentLeptons()[0].momentum()) > 40*GeV && nu.pT()>25*GeV) { keepevent=true; lepton_charge=wfinder_born_mu.constituentLeptons()[0].charge(); lepton_eta=fabs(wfinder_born_mu.constituentLeptons()[0].pseudorapidity()); } } if (!keepevent) { MSG_DEBUG("event does not pass mT and MET cuts"); vetoEvent; } if (lepton_charge>0) { _hist_wplus_winc->fill(10., weight); _hist_wplus_winc->fill(16., weight); _hist_wplus_winc->fill(30., weight); _hist_wplus_winc->fill(60., weight); _hist_w_inc->fill(+1, weight); } else if (lepton_charge<0) { _hist_wminus_winc->fill(10., weight); _hist_wminus_winc->fill(16., weight); _hist_wminus_winc->fill(30., weight); _hist_wminus_winc->fill(60., weight); _hist_w_inc->fill(-1, weight); } // Find hadrons in the event const UnstableFinalState& fs = applyProjection(event, "hadrons"); /// FIND Different channels // 1: wcjet // get jets const Jets& jets = applyProjection(event, "jets").jetsByPt(25.0*GeV); // loop over jets to select jets used to match to charm Jets js; int matched_charmHadron=0; double charm_charge=0; int njets=0; int nj=0; int njph=-1; bool mat_jet=false; double ptcharm=0; if (matched_charmHadron>-1 ) { foreach (const Jet& j, jets) { mat_jet=false; if (fabs(j.eta()) > 2.5) continue; if (fabs(j.pt()) <= 25*GeV ) continue; njets++; foreach (const Particle& p, fs.particles()) { const GenParticle* part = p.genParticle(); if ( p.hasCharm()) { if (isFromBDecay(p)) continue; if (p.pT()<5*GeV) continue; if (hasCharmedChildren(part))continue; if ( deltaR(p.momentum() , j.momentum()) < 0.3 ) { mat_jet=true; if ( p.pT()>ptcharm) {charm_charge=part->pdg_id();ptcharm=p.pT();} } } } if (mat_jet) nj++; } if (charm_charge*lepton_charge > 0 ) charge_weight=-1; else charge_weight=+1; if (nj==1) { if (lepton_charge>0) _hist_wcjet_charge->fill(+1, weight*charge_weight); else if (lepton_charge<0) _hist_wcjet_charge->fill(-1, weight*charge_weight); if (lepton_charge>0) _hist_wplus_wcjet_eta_lep->fill(lepton_eta, weight*charge_weight); else if (lepton_charge<0) _hist_wminus_wcjet_eta_lep->fill(lepton_eta, weight*charge_weight); _hist_wcjet_jets->fill(njets, weight*charge_weight); // if (lepton_charge>0) _hist_wcjet_jets_plus->fill(njets, weight*charge_weight); // else if (lepton_charge<0) _hist_wcjet_jets_minus->fill(njets, weight*charge_weight); // check_wc_weight=check_wc_weight+charge_weight; } } // // 1/2: w+d(*) meson foreach (const Particle& p, fs.particles()) { const GenParticle* part = p.genParticle(); if (p.pT()<8*GeV) continue; if (fabs(p.eta())>2.2) continue; // W+D if ( fabs(part->pdg_id()) == 411) { if (lepton_charge*part->pdg_id() > 0) charge_weight=-1; else charge_weight=+1; // fill histos if (lepton_charge>0) _hist_wd_charge->fill(+1, weight*charge_weight); else if (lepton_charge<0) _hist_wd_charge->fill(-1, weight*charge_weight); if (lepton_charge>0) _hist_wd_plus->fill(0, weight*charge_weight); else if (lepton_charge<0) _hist_wd_minus->fill(0, weight*charge_weight); if (lepton_charge>0) _hist_wplus_wdminus_eta_lep->fill(lepton_eta, weight*charge_weight); else if (lepton_charge<0) _hist_wminus_wdplus_eta_lep->fill(lepton_eta, weight*charge_weight); if (lepton_charge>0) _hist_wplusd_wplusinc_pt->fill(p.pT(), weight*charge_weight); else if (lepton_charge<0) _hist_wminusd_wminusinc_pt->fill(p.pT(), weight*charge_weight); // fill histos // if (lepton_charge>0) _hist_wd_inte->fill(+1, weight*charge_weight); // else if (lepton_charge<0) _hist_wd_inte->fill(-1, weight*charge_weight); } // W+Dstar if ( fabs(part->pdg_id()) == 413 ) { if (lepton_charge*part->pdg_id() > 0) charge_weight=-1; else charge_weight=+1; // fill histos if (lepton_charge>0) _hist_wdstar_charge->fill(+1, weight*charge_weight); else if (lepton_charge<0) _hist_wdstar_charge->fill(-1, weight*charge_weight); if (lepton_charge>0) _hist_wd_plus->fill(0, weight*charge_weight); else if (lepton_charge<0) _hist_wd_minus->fill(0, weight*charge_weight); if (lepton_charge>0) _hist_wplus_wdstar_eta_lep->fill(lepton_eta, weight*charge_weight); else if (lepton_charge<0) _hist_wminus_wdstar_eta_lep->fill(lepton_eta, weight*charge_weight); if (lepton_charge>0) _hist_wplusd_wplusinc_pt->fill(p.pT(), weight*charge_weight); else if (lepton_charge<0) _hist_wminusd_wminusinc_pt->fill(p.pT(), weight*charge_weight); } } }// end of analyse function /// Normalise histograms etc., after the run void finalize() { /// @todo Normalise, scale and otherwise manipulate histograms here // norm to cross section // d01 scale( _hist_wcjet_charge, crossSection()/sumOfWeights()); scale(_hist_wd_charge, crossSection()/sumOfWeights()); scale(_hist_wdstar_charge, crossSection()/sumOfWeights()); //d02 scale(_hist_wcjet_plus , crossSection()/sumOfWeights()); scale(_hist_wcjet_minus , crossSection()/sumOfWeights()); scale(_hist_wd_plus , crossSection()/sumOfWeights()); scale(_hist_wd_minus , crossSection()/sumOfWeights()); divide( _hist_wcjet_plus , _hist_wcjet_minus , _hist_wcjet_ratio ); divide( _hist_wd_plus , _hist_wd_minus , _hist_wd_ratio ); //d03 scale(_hist_wplus_wcjet_eta_lep, crossSection()/sumOfWeights()); scale(_hist_wminus_wcjet_eta_lep , crossSection()/sumOfWeights()); //d04 scale(_hist_wplus_wdminus_eta_lep, crossSection()/sumOfWeights()); scale(_hist_wminus_wdplus_eta_lep, crossSection()/sumOfWeights()); scale(_hist_wplus_wdstar_eta_lep, crossSection()/sumOfWeights()); scale(_hist_wminus_wdstar_eta_lep, crossSection()/sumOfWeights()); //d05 scale( _hist_w_inc , crossSection()/sumOfWeights()); // _hist_wdstar_winc //d06 scale( _hist_wplusd_wplusinc_pt, crossSection()/sumOfWeights()); scale( _hist_wminusd_wminusinc_pt, crossSection()/sumOfWeights()); scale( _hist_wplusdstar_wplusinc_pt, crossSection()/sumOfWeights()); scale( _hist_wminusdstar_wminusinc_pt , crossSection()/sumOfWeights()); scale( _hist_wplus_winc , crossSection()/sumOfWeights()); scale( _hist_wminus_winc, crossSection()/sumOfWeights()); //d07 scale(_hist_wcjet_jets , crossSection()/sumOfWeights()); //d08 //_hist_wcjet_jets_ratio scale(_hist_wcjet_jets_minus, crossSection()/sumOfWeights()); scale(_hist_wcjet_jets_plus, crossSection()/sumOfWeights()); divide( _hist_wcjet_jets_plus , _hist_wcjet_jets_minus , _hist_wcjet_jets_ratio ); } //@} private: // Data members like post-cuts event weight counters go here // Check whether particle comes from b-decay bool isFromBDecay(const Particle& p) { bool isfromB=false; if (p.genParticle() == NULL) return false; const GenParticle* part = p.genParticle(); GenVertex* ivtx = part->production_vertex(); while(ivtx) { if (ivtx->particles_in_size() < 1) { isfromB = false; break; }; const HepMC::GenVertex::particles_in_const_iterator iPart_invtx = ivtx->particles_in_const_begin(); part = (*iPart_invtx); if ( !(part) ) { isfromB = false; break; }; isfromB = PID::hasBottom(part->pdg_id()); if (isfromB==true) break; ivtx = part->production_vertex(); if ( (part->pdg_id() == 2212) || !(ivtx) ) break; // reached beam }; return (isfromB); } // Check whether particle has charmed children bool hasCharmedChildren(const GenParticle *part) { bool hasCharmedChild=false; if (part == NULL) return false; GenVertex* ivtx = part->end_vertex(); if (ivtx == NULL) return false; // if (ivtx->particles_out_size() < 2) return false; HepMC::GenVertex::particles_out_const_iterator iPart_invtx = ivtx->particles_out_const_begin(); HepMC::GenVertex::particles_out_const_iterator end_invtx = ivtx->particles_out_const_end(); for ( ; iPart_invtx != end_invtx; iPart_invtx++ ) { const GenParticle* p2 = (*iPart_invtx); if (p2==part) continue; hasCharmedChild = PID::hasCharm(p2->pdg_id()); if (hasCharmedChild==true) break; hasCharmedChild = hasCharmedChildren(p2); if (hasCharmedChild==true) break; } return (hasCharmedChild); } private: /// @name Histograms //@{ int check_passW; int check_passD; int check_passDstar; int check_passWc; int check_all; int check_weight; int check_eta; int check_wc_weight; int check_wc; int check_wdstar; int check_wdplus; int events; int check_0jet; int check_1jet; int check_2jet; int check_3jet; int check_1wcjet; int check_2wcjet; //d01-x01- Histo1DPtr _hist_wcjet_charge ; Histo1DPtr _hist_wd_charge; Histo1DPtr _hist_wdstar_charge; //d02-x01- Scatter2DPtr _hist_wcjet_ratio; Scatter2DPtr _hist_wd_ratio; // Histo1DPtr _hist_wdstar_inte; Histo1DPtr _hist_wcjet_plus; Histo1DPtr _hist_wd_plus; Histo1DPtr _hist_wcjet_minus; Histo1DPtr _hist_wd_minus; //d03-x01- Histo1DPtr _hist_wplus_wcjet_eta_lep; Histo1DPtr _hist_wminus_wcjet_eta_lep; //d04-x01- Histo1DPtr _hist_wplus_wdminus_eta_lep; Histo1DPtr _hist_wminus_wdplus_eta_lep; //d05-x01- Histo1DPtr _hist_wplus_wdstar_eta_lep; Histo1DPtr _hist_wminus_wdstar_eta_lep; // postprocessing histos //d05-x01 Histo1DPtr _hist_wplus_winc;// Histo1DPtr _hist_wminus_winc;// Histo1DPtr _hist_w_inc; //d06-x01 Histo1DPtr _hist_wplusd_wplusinc_pt ; Histo1DPtr _hist_wminusd_wminusinc_pt; Histo1DPtr _hist_wplusdstar_wplusinc_pt; Histo1DPtr _hist_wminusdstar_wminusinc_pt; Histo1DPtr _hist_wplusd_wplusinc_pt_winc ; Histo1DPtr _hist_wminusd_wminusinc_pt_winc; Histo1DPtr _hist_wplusdstar_wplusinc_pt_winc; Histo1DPtr _hist_wminusdstar_wminusinc_pt_winc; // d07-x01 Histo1DPtr _hist_wcjet_jets ; //d08-x01 Scatter2DPtr _hist_wcjet_jets_ratio ; Histo1DPtr _hist_wcjet_jets_plus ; Histo1DPtr _hist_wcjet_jets_minus; //d09-x01 // Histo1DPtr _hist_wsoftmuon_jets; //@} }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_2014_I1282447); }