// -*- C++ -*- #include "Rivet/Analysis.hh" #include "Rivet/Projections/FinalState.hh" #include "Rivet/Projections/ZFinder.hh" #include "Rivet/Projections/FastJets.hh" #include "Rivet/Projections/HeavyHadrons.hh" #include "Rivet/Projections/VetoedFinalState.hh" namespace Rivet { class ATLAS_Zbb : public Analysis { public: /// @name Constructors etc. //@{ /// Constructors ATLAS_Zbb(std::string name="ATLAS_Zbb") : Analysis(name) { _mode = 1; setNeedsCrossSection(true); } //@} public: /// @name Analysis methods //@{ /// Book histograms and initialise projections before the run void init() { FinalState fs; Cut cuts = Cuts::etaIn(-2.5,2.5) & (Cuts::pT > 20.0*GeV); ZFinder zfinder(fs, cuts, _mode==1? PID::ELECTRON : PID::MUON, 76.0*GeV, 106.0*GeV, 0.1, ZFinder::CLUSTERNODECAY, ZFinder::NOTRACK); addProjection(zfinder, "ZFinder"); //FastJets jetpro1( getProjection("ZFinder").remainingFinalState(), FastJets::ANTIKT, 0.4); VetoedFinalState jet_fs(fs); jet_fs.addVetoOnThisFinalState(getProjection("ZFinder")); FastJets jetpro1(jet_fs, FastJets::ANTIKT, 0.4); jetpro1.useInvisibles(); addProjection(jetpro1, "AntiKtJets04"); // addProjection(HeavyHadrons(), "BHadrons"); //Histograms with data binning _h_bjet_Pt = bookHisto1D( 3, 1, 1); _h_bjet_Y = bookHisto1D( 5, 1, 1); _h_bjet_Yboost = bookHisto1D( 7, 1, 1); _h_bjet_DY20 = bookHisto1D( 9, 1, 1); _h_bjet_ZdPhi20 = bookHisto1D(11, 1, 1); _h_bjet_ZdR20 = bookHisto1D(13, 1, 1); _h_bjet_ZPt = bookHisto1D(15, 1, 1); _h_bjet_ZY = bookHisto1D(17, 1, 1); _h_2bjet_dR = bookHisto1D(21, 1, 1); _h_2bjet_Mbb = bookHisto1D(23, 1, 1); _h_2bjet_ZPt = bookHisto1D(25, 1, 1); _h_2bjet_ZY = bookHisto1D(27, 1, 1); _h_hadron_pT = bookHisto1D("hadlev_B_hadron_pT", 50,0,100); _h_matching_hadron_pT = bookHisto1D("hadlev_B_hadron_pT_matching_jet", 50,0,100); _h_bjet_pT_single = bookHisto1D("hadlev_bjet_pT_single_hadron_0GeV", 20,0,200); _h_bjet_pT_double = bookHisto1D("hadlev_bjet_pT_double_hadron_0GeV", 20,0,200); _h_bjet_pT_single_5GeV = bookHisto1D("hadlev_bjet_pT_single_hadron_5GeV", 20,0,200); _h_bjet_pT_double_5GeV = bookHisto1D("hadlev_bjet_pT_double_hadron_5GeV", 20,0,200); _h_bjet_num_matching = bookHisto1D("hadlev_bjet_number_matching_B_hadrons", 6,-0.5,5.5); _h_bjet_num_matching_5GeV = bookHisto1D("hadlev_bjet_number_matching_B_hadrons_5GeV", 6,-0.5,5.5); _h_bjet_bb_over_b = bookScatter2D("hadlev_bjet_bb_over_b_hadron"); _h_bjet_bb_over_b_5GeV = bookScatter2D("hadlev_bjet_bb_over_b_hadron_5GeV"); _h_quark_pT = bookHisto1D("quarklev_b_quark_pT", 50,0,100); _h_matching_quark_pT = bookHisto1D("quarklev_b_quark_pT_matching_jet", 50,0,100); _h_bjet_pT_single_quark = bookHisto1D("quarklev_bjet_pT_single_quark_0GeV", 20,0,200); _h_bjet_pT_double_quark = bookHisto1D("quarklev_bjet_pT_double_quark_0GeV", 20,0,200); _h_bjet_pT_single_quark_5GeV = bookHisto1D("quarklev_bjet_pT_single_quark_5GeV", 20,0,200); _h_bjet_pT_double_quark_5GeV = bookHisto1D("quarklev_bjet_pT_double_quark_5GeV", 20,0,200); _h_bjet_num_matching_quark = bookHisto1D("quarklev_bjet_number_matching_b_quarks", 6,-0.5,5.5); _h_bjet_num_matching_quark_5GeV = bookHisto1D("quarklev_bjet_number_matching_b_quarks_5GeV", 6,-0.5,5.5); _h_bjet_bb_over_b_quark = bookScatter2D("quarklev_bjet_bb_over_b_quarks"); _h_bjet_bb_over_b_quark_5GeV = bookScatter2D("quarklev_bjet_bb_over_b_quarks_5GeV"); } //========================================================================================== /// Perform the per-event analysis void analyze(const Event& e) { //--------------------------- const double weight = e.weight(); // -- check we have a Z: const ZFinder& zfinder = applyProjection(e, "ZFinder"); if(zfinder.bosons().size() != 1) vetoEvent; const ParticleVector boson_s = zfinder.bosons(); const Particle boson_f = boson_s[0] ; const ParticleVector zleps = zfinder.constituents(); //--------------------------- //--------------------------- Particles weakBs, bquarks; // Loop over the unstable particles foreach (const Particle& p, particles(e.genEvent()) ) { const PdgId pid = p.pid(); //look for b quarks: if(abs(pid)==5) { bool good_B = true; const HepMC::GenParticle* pgen = p.genParticle(); HepMC::GenVertex* vgen = pgen -> end_vertex(); if(vgen) { //make sure this is the last b in the chain: for (HepMC::GenVertex::particles_out_const_iterator it = vgen->particles_out_const_begin(); it != vgen->particles_out_const_end(); ++it) { if ( abs( (*it)->pdg_id() ) ==5 ) { good_B = false; } } } if(good_B) { bquarks += p; } } // Look for B hadrons: if (PID::hasBottom(pid)) { bool good_B = true; const HepMC::GenParticle* pgen = p.genParticle(); HepMC::GenVertex* vgen = pgen -> end_vertex(); // make sure there are no more B hadrons in the chain: for (HepMC::GenVertex::particles_out_const_iterator it = vgen->particles_out_const_begin(); it != vgen->particles_out_const_end(); ++it) { if ( PID::hasBottom( (*it)->pdg_id() ) ) { good_B = false; } } if ( good_B ) { weakBs += p; } } }//loop over particles if (weakBs.size() <1 && bquarks.size() <1 ) vetoEvent; //--------------------------- // -- get the b-jets: const Jets& jets = applyProjection(e, "AntiKtJets04").jetsByPt(Cuts::pT >20.0*GeV && Cuts::absrap <2.4); Jets good_jets, b_jets, b_quark_jets; vector hadron_used; for(uint i=0; i2.4) continue; //veto overlaps with Z leptons: bool veto = false; foreach(const Particle& zlep, zleps) { if(deltaR(jet, zlep) < 0.5) veto = true; } if(veto) continue; good_jets.push_back(jet); int nmatch=0, nmatch_5GeV=0, nmatch_quark=0, nmatch_quark_5GeV=0; int iB=-1; foreach(const Particle& bhadron, weakBs) { iB++; if(hadron_used[iB]) continue; if( deltaR(jet, bhadron) <= 0.3 ) { hadron_used[iB]=true; nmatch++; if(bhadron.pT()>5*GeV) { if(nmatch_5GeV==0) b_jets.push_back(jet); nmatch_5GeV++; } } } // end loop on b-hadrons foreach(const Particle& bquark, bquarks) { if( deltaR(jet, bquark) <= 0.3 ) { nmatch_quark++; if(bquark.pT()>5*GeV) nmatch_quark_5GeV++; } } // end loop on b-quarks _h_bjet_num_matching->fill(nmatch); _h_bjet_num_matching_5GeV->fill(nmatch_5GeV); if(nmatch == 1) _h_bjet_pT_single->fill(jet.pT()/GeV); if(nmatch > 1) _h_bjet_pT_double->fill(jet.pT()/GeV); if(nmatch_5GeV == 1) _h_bjet_pT_single_5GeV->fill(jet.pT()/GeV); if(nmatch_5GeV > 1) _h_bjet_pT_double_5GeV->fill(jet.pT()/GeV); _h_bjet_num_matching_quark->fill(nmatch_quark); _h_bjet_num_matching_quark_5GeV->fill(nmatch_quark_5GeV); if(nmatch_quark == 1) _h_bjet_pT_single_quark->fill(jet.pT()/GeV); if(nmatch_quark > 1) _h_bjet_pT_double_quark->fill(jet.pT()/GeV); if(nmatch_quark_5GeV == 1) _h_bjet_pT_single_quark_5GeV->fill(jet.pT()/GeV); if(nmatch_quark_5GeV > 1) _h_bjet_pT_double_quark_5GeV->fill(jet.pT()/GeV); } //now loop the other way, to avoid double counting matched hadrons: foreach(const Particle& bhadron, weakBs) { if( fabs(bhadron.rapidity()) > 2.7) continue; _h_hadron_pT->fill(bhadron.pT(), weight); bool match = false; foreach(const Jet& jet, good_jets) if( deltaR(jet, bhadron) <= 0.3 ) match=true; if( match ) _h_matching_hadron_pT->fill(bhadron.pT()*GeV, weight); } foreach(const Particle& bquark, bquarks) { if( fabs(bquark.rapidity()) > 2.7) continue; _h_quark_pT->fill(bquark.pT(), weight); bool match = false; foreach(const Jet& jet, good_jets) if( deltaR(jet, bquark) <= 0.3 ) match=true; if( match ) _h_matching_quark_pT->fill(bquark.pT()*GeV, weight); } //--------------------------- //and make sure we have at least 1: if(b_jets.empty()) vetoEvent; //--------------------------- // fill the plots: const double ZpT = boson_f.pT()/GeV; const double ZY = boson_f.absrap(); _h_bjet_ZPt->fill(ZpT, weight); _h_bjet_ZY ->fill(ZY, weight); foreach(const Jet& jet, b_jets) { _h_bjet_Pt->fill(jet.pT()/GeV, weight ); _h_bjet_Y ->fill(jet.absrap(), weight ); const double Yboost = 0.5 * fabs(boson_f.rapidity() + jet.rapidity()); _h_bjet_Yboost->fill(Yboost, weight ); if(ZpT > 20.) { const double ZBDY = fabs( boson_f.rapidity() - jet.rapidity() ); const double ZBDPHI = fabs( deltaPhi(jet.phi(), boson_f.phi()) ); const double ZBDR = deltaR(jet, boson_f, RAPIDITY); _h_bjet_DY20->fill( ZBDY, weight); _h_bjet_ZdPhi20->fill(ZBDPHI, weight); _h_bjet_ZdR20->fill( ZBDR, weight); } } //loop over b-jets if (b_jets.size() < 2) return; _h_2bjet_ZPt->fill(ZpT, weight); _h_2bjet_ZY ->fill(ZY, weight); const double BBDR = deltaR(b_jets[0], b_jets[1], RAPIDITY); const double Mbb = (b_jets[0].momentum() + b_jets[1].momentum()).mass(); _h_2bjet_dR ->fill(BBDR, weight); _h_2bjet_Mbb->fill(Mbb, weight); } // end of analysis loop /// Normalise histograms etc., after the run void finalize() { const double normfac = crossSection() / sumOfWeights(); scale( _h_bjet_Pt, normfac); scale( _h_bjet_Y, normfac); scale( _h_bjet_Yboost, normfac); scale( _h_bjet_DY20, normfac); scale( _h_bjet_ZdPhi20, normfac); scale( _h_bjet_ZdR20, normfac); scale( _h_bjet_ZPt, normfac); scale( _h_bjet_ZY, normfac); scale( _h_2bjet_dR, normfac); scale( _h_2bjet_Mbb, normfac); scale( _h_2bjet_ZPt, normfac); scale( _h_2bjet_ZY, normfac); scale( _h_hadron_pT, normfac); scale( _h_matching_hadron_pT, normfac); scale( _h_bjet_pT_single, normfac); scale( _h_bjet_pT_double, normfac); scale( _h_bjet_pT_single_5GeV, normfac); scale( _h_bjet_pT_double_5GeV, normfac); scale( _h_bjet_num_matching, normfac); scale( _h_bjet_num_matching_5GeV, normfac); divide( _h_bjet_pT_double, _h_bjet_pT_single, _h_bjet_bb_over_b); divide( _h_bjet_pT_double_5GeV, _h_bjet_pT_single_5GeV, _h_bjet_bb_over_b_5GeV); scale( _h_quark_pT, normfac); scale( _h_matching_quark_pT, normfac); scale( _h_bjet_pT_single_quark, normfac); scale( _h_bjet_pT_double_quark, normfac); scale( _h_bjet_pT_single_quark_5GeV, normfac); scale( _h_bjet_pT_double_quark_5GeV, normfac); scale( _h_bjet_num_matching_quark, normfac); scale( _h_bjet_num_matching_quark_5GeV, normfac); divide( _h_bjet_pT_double_quark, _h_bjet_pT_single_quark, _h_bjet_bb_over_b_quark); divide( _h_bjet_pT_double_quark_5GeV, _h_bjet_pT_single_quark_5GeV, _h_bjet_bb_over_b_quark_5GeV); } //@} protected: // Data members like post-cuts event weight counters go here size_t _mode; private: Histo1DPtr _h_bjet_Pt; Histo1DPtr _h_bjet_Y; Histo1DPtr _h_bjet_Yboost; Histo1DPtr _h_bjet_DY20; Histo1DPtr _h_bjet_ZdPhi20; Histo1DPtr _h_bjet_ZdR20; Histo1DPtr _h_bjet_ZPt; Histo1DPtr _h_bjet_ZY; Histo1DPtr _h_2bjet_dR; Histo1DPtr _h_2bjet_Mbb; Histo1DPtr _h_2bjet_ZPt; Histo1DPtr _h_2bjet_ZY; Histo1DPtr _h_hadron_pT; Histo1DPtr _h_matching_hadron_pT; Histo1DPtr _h_bjet_pT_single; Histo1DPtr _h_bjet_pT_double; Histo1DPtr _h_bjet_pT_single_5GeV; Histo1DPtr _h_bjet_pT_double_5GeV; Histo1DPtr _h_bjet_num_matching; Histo1DPtr _h_bjet_num_matching_5GeV; Scatter2DPtr _h_bjet_bb_over_b; Scatter2DPtr _h_bjet_bb_over_b_5GeV; Histo1DPtr _h_quark_pT; Histo1DPtr _h_matching_quark_pT; Histo1DPtr _h_bjet_pT_single_quark; Histo1DPtr _h_bjet_pT_double_quark; Histo1DPtr _h_bjet_pT_single_quark_5GeV; Histo1DPtr _h_bjet_pT_double_quark_5GeV; Histo1DPtr _h_bjet_num_matching_quark; Histo1DPtr _h_bjet_num_matching_quark_5GeV; Scatter2DPtr _h_bjet_bb_over_b_quark; Scatter2DPtr _h_bjet_bb_over_b_quark_5GeV; }; class ATLAS_Zbb_EL : public ATLAS_Zbb { public: ATLAS_Zbb_EL() : ATLAS_Zbb("ATLAS_Zbb_EL") { _mode = 1; } }; class ATLAS_Zbb_MU : public ATLAS_Zbb { public: ATLAS_Zbb_MU() : ATLAS_Zbb("ATLAS_Zbb_MU") { _mode = 2; } }; // The hook for the plugin system DECLARE_RIVET_PLUGIN(ATLAS_Zbb); DECLARE_RIVET_PLUGIN(ATLAS_Zbb_MU); DECLARE_RIVET_PLUGIN(ATLAS_Zbb_EL); }