Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:32:11

0001 #include "Validation/EventGenerator/interface/TTbar_Kinematics.h"
0002 #include "DQMServices/Core/interface/DQMStore.h"
0003 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
0004 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0005 #include "Validation/EventGenerator/interface/DQMHelper.h"
0006 
0007 using namespace edm;
0008 TTbar_Kinematics::TTbar_Kinematics(const edm::ParameterSet &iConfig)
0009     : hepmcCollection_(iConfig.getParameter<edm::InputTag>("hepmcCollection")),
0010       genEventInfoProductTag_(iConfig.getParameter<edm::InputTag>("genEventInfoProductTag")) {
0011   hepmcCollectionToken_ = consumes<HepMCProduct>(hepmcCollection_);
0012   genEventInfoProductTagToken_ = consumes<GenEventInfoProduct>(genEventInfoProductTag_);
0013 }
0014 
0015 TTbar_Kinematics::~TTbar_Kinematics() {}
0016 
0017 //
0018 // member functions
0019 //
0020 
0021 // ------------ method called for each event  ------------
0022 void TTbar_Kinematics::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0023   // --- the MC weights ---
0024   Handle<GenEventInfoProduct> evt_info;
0025   iEvent.getByToken(genEventInfoProductTagToken_, evt_info);
0026   if (!evt_info.isValid())
0027     return;
0028   weight = evt_info->weight();
0029 
0030   ///Gathering the HepMCProduct information
0031   edm::Handle<HepMCProduct> evt;
0032   iEvent.getByToken(hepmcCollectionToken_, evt);
0033 
0034   //Get EVENT
0035   const HepMC::GenEvent *myGenEvent = evt->GetEvent();
0036 
0037   TLorentzVector tlv_Top, tlv_TopBar, tlv_Bottom, tlv_BottomBar, tlv_Wplus, tlv_Wmin, tlv_TTbar;
0038   bool top(false), antitop(false), antibottom(false), bottom(false), Wplus(false), Wmin(false);
0039   for (HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin();
0040        iter != myGenEvent->particles_end();
0041        iter++) {
0042     if ((*iter)->pdg_id() == PdtPdgMini::t || (*iter)->pdg_id() == PdtPdgMini::anti_t) {
0043       if ((*iter)->end_vertex()) {
0044         HepMC::GenVertex::particle_iterator des;
0045         for (des = (*iter)->end_vertex()->particles_begin(HepMC::children);
0046              des != (*iter)->end_vertex()->particles_end(HepMC::children);
0047              ++des) {
0048           if ((*des)->pdg_id() == PdtPdgMini::b) {
0049             tlv_Bottom.SetPxPyPzE(
0050                 (*des)->momentum().px(), (*des)->momentum().py(), (*des)->momentum().pz(), (*des)->momentum().e());
0051             bottom = true;
0052           }
0053           if ((*des)->pdg_id() == PdtPdgMini::anti_b) {
0054             antibottom = true;
0055             tlv_BottomBar.SetPxPyPzE(
0056                 (*des)->momentum().px(), (*des)->momentum().py(), (*des)->momentum().pz(), (*des)->momentum().e());
0057           }
0058           if ((*des)->pdg_id() == PdtPdgMini::W_plus) {
0059             tlv_TopBar.SetPxPyPzE(
0060                 (*iter)->momentum().px(), (*iter)->momentum().py(), (*iter)->momentum().pz(), (*iter)->momentum().e());
0061             antitop = true;
0062             tlv_Wplus.SetPxPyPzE(
0063                 (*des)->momentum().px(), (*des)->momentum().py(), (*des)->momentum().pz(), (*des)->momentum().e());
0064             Wplus = true;
0065           }
0066           if ((*des)->pdg_id() == PdtPdgMini::W_minus) {
0067             tlv_Top.SetPxPyPzE(
0068                 (*iter)->momentum().px(), (*iter)->momentum().py(), (*iter)->momentum().pz(), (*iter)->momentum().e());
0069             top = true;
0070             tlv_Wmin.SetPxPyPzE(
0071                 (*des)->momentum().px(), (*des)->momentum().py(), (*des)->momentum().pz(), (*des)->momentum().e());
0072             Wmin = true;
0073           }
0074         }
0075       }
0076     }
0077   }
0078 
0079   tlv_TTbar = tlv_Top + tlv_TopBar;
0080 
0081   //---topquarkquantities---
0082   nEvt->Fill(0.5, weight);
0083   if (top && antitop) {
0084     hTopPt->Fill(tlv_Top.Pt(), weight);
0085     hTopPt->Fill(tlv_TopBar.Pt(), weight);
0086 
0087     hTopY->Fill(tlv_Top.Rapidity(), weight);
0088     hTopY->Fill(tlv_TopBar.Rapidity(), weight);
0089 
0090     hTopMass->Fill(tlv_Top.M(), weight);
0091     hTopMass->Fill(tlv_TopBar.M(), weight);
0092 
0093     //---ttbarpairquantities---
0094     hTTbarPt->Fill(tlv_TTbar.Pt(), weight);
0095     hTTbarY->Fill(tlv_TTbar.Rapidity(), weight);
0096     hTTbarMass->Fill(tlv_TTbar.M(), weight);
0097   }
0098   if (bottom && antibottom) {
0099     hBottomPt->Fill(tlv_Bottom.Pt(), weight);
0100     hBottomPt->Fill(tlv_BottomBar.Pt(), weight);
0101 
0102     hBottomEta->Fill(tlv_Bottom.Eta(), weight);
0103     hBottomEta->Fill(tlv_BottomBar.Eta(), weight);
0104 
0105     //hBottomY->Fill(math::XYZTLorentzVector(bottom->momentum()).Rapidity(),weight);
0106     //hBottomY->Fill(math::XYZTLorentzVector(antibottom->momentum()).Rapidity(),weight);
0107 
0108     hBottomY->Fill(tlv_Bottom.Rapidity(), weight);
0109     hBottomY->Fill(tlv_BottomBar.Rapidity(), weight);
0110 
0111     hBottomPz->Fill(tlv_Bottom.Pz(), weight);
0112     hBottomPz->Fill(tlv_BottomBar.Pz(), weight);
0113 
0114     hBottomE->Fill(tlv_Bottom.E(), weight);
0115     hBottomE->Fill(tlv_BottomBar.E(), weight);
0116 
0117     hBottomMass->Fill(tlv_Bottom.M(), weight);
0118     hBottomMass->Fill(tlv_BottomBar.M(), weight);
0119 
0120     hBottomPtPz->Fill(tlv_Bottom.Pt(), tlv_Bottom.Pz(), weight);
0121     hBottomPtPz->Fill(tlv_BottomBar.Pt(), tlv_BottomBar.Pz(), weight);
0122 
0123     hBottomEtaPz->Fill(tlv_Bottom.Eta(), tlv_Bottom.Pz(), weight);
0124     hBottomEtaPz->Fill(tlv_BottomBar.Eta(), tlv_BottomBar.Pz(), weight);
0125 
0126     hBottomEtaPt->Fill(tlv_Bottom.Eta(), tlv_Bottom.Pt(), weight);
0127     hBottomEtaPt->Fill(tlv_BottomBar.Eta(), tlv_BottomBar.Pt(), weight);
0128 
0129     hBottomYPz->Fill(tlv_Bottom.Rapidity(), tlv_Bottom.Pz(), weight);
0130     hBottomYPz->Fill(tlv_BottomBar.Rapidity(), tlv_BottomBar.Pz(), weight);
0131 
0132     hBottomMassPz->Fill(tlv_Bottom.M(), tlv_Bottom.Pz(), weight);
0133     hBottomMassPz->Fill(tlv_BottomBar.M(), tlv_BottomBar.Pz(), weight);
0134 
0135     hBottomMassEta->Fill(tlv_Bottom.M(), tlv_Bottom.Eta(), weight);
0136     hBottomMassEta->Fill(tlv_BottomBar.M(), tlv_BottomBar.Eta(), weight);
0137 
0138     hBottomMassY->Fill(tlv_Bottom.M(), tlv_Bottom.Rapidity(), weight);
0139     hBottomMassY->Fill(tlv_BottomBar.M(), tlv_BottomBar.Rapidity(), weight);
0140 
0141     hBottomMassDeltaY->Fill(tlv_Bottom.M(), tlv_Bottom.Eta() - tlv_Bottom.Rapidity(), weight);
0142     hBottomMassDeltaY->Fill(tlv_BottomBar.M(), tlv_BottomBar.Eta() - tlv_BottomBar.Rapidity(), weight);
0143   }
0144   if (Wplus && Wmin) {
0145     hWplusPz->Fill(tlv_Wplus.Pz(), weight);
0146     hWminPz->Fill(tlv_Wmin.Pz(), weight);
0147   }
0148 }
0149 
0150 // ------------ method called once each job just before starting event loop  ------------
0151 void TTbar_Kinematics::bookHistograms(DQMStore::IBooker &i, edm::Run const &r, edm::EventSetup const &e) {
0152   DQMHelper dqm(&i);
0153   i.setCurrentFolder("Generator/TTbar");
0154 
0155   nEvt = dqm.book1dHisto("nEvt", "n analyzed Events", 1, 0., 1., "bins", "Number of Events");
0156 
0157   hTopPt = dqm.book1dHisto(
0158       "TTbar_TopPt", "t quark transverse momentum", 1000, 0., 1000., "P_{t}^{t quark} (GeV)", "Number of Events");
0159   hTopY = dqm.book1dHisto("TTbar_TopY", "t quark rapidity", 200, -5., 5., "Y_{t quark}", "Number of Events");
0160   hTopMass = dqm.book1dHisto("TTbar_TopMass", "t quark mass", 500, 0., 500., "M_{t quark} (GeV)", "Number of Events");
0161 
0162   hTTbarPt = dqm.book1dHisto(
0163       "TTbar_TTbarPt", "tt pair transverse momentum", 1000, 0., 1000., "P_{t}^{tt pair} (GeV)", "Number of Events");
0164   hTTbarY = dqm.book1dHisto("TTbar_TTbarY", "tt pair rapidity", 200, -5., 5., "Y_{tt pair}", "Number of Events");
0165   hTTbarMass =
0166       dqm.book1dHisto("TTbar_TTbarMass", "tt pair mass", 1000, 0., 1000., "M_{tt pair} (GeV)", "Number of Events");
0167 
0168   hBottomPt = dqm.book1dHisto(
0169       "TTbar_BottomPt", "b quark transverse momentum", 1000, 0., 1000., "P_{t}^{b quark} (GeV)", "Number of Events");
0170   hBottomEta = dqm.book1dHisto(
0171       "TTbar_BottomEta", "b quark pseudo-rapidity", 200, -5., 5., "#eta_{b quark} (GeV)", "Number of Events");
0172   hBottomY =
0173       dqm.book1dHisto("TTbar_BottomY", "b quark rapidity", 200, -5., 5., "M_{b quark} (GeV)", "Number of Events");
0174   hBottomPz = dqm.book1dHisto(
0175       "TTbar_BottomPz", "b quark longitudinal momentum", 200, -100., 100., "P_{z}^{b quark} (GeV)", "Number of Events");
0176   hBottomE =
0177       dqm.book1dHisto("TTbar_BottomE", "b quark energy", 1000, 0., 1000., "E_{b quark} (GeV)", "Number of Events");
0178   hBottomMass =
0179       dqm.book1dHisto("TTbar_BottomMass", "b quark mass", 50, 0., 5., "M_{b quark} (GeV)", "Number of Events");
0180 
0181   hBottomPtPz = dqm.book2dHisto("TTbar_BottomPtPz",
0182                                 "b quark longitudinal vs transverse momentum",
0183                                 1000,
0184                                 0.,
0185                                 1000.,
0186                                 200,
0187                                 -100.,
0188                                 100.,
0189                                 "P_{z}^{b quark} (GeV)",
0190                                 "P_{t}^{b quark} (GeV)");
0191   hBottomEtaPz = dqm.book2dHisto("TTbar_BottomEtaPz",
0192                                  "b quark longitudinal momentum vs pseudorapidity",
0193                                  200,
0194                                  -5.,
0195                                  5.,
0196                                  200,
0197                                  -100.,
0198                                  100.,
0199                                  "#eta_{b quark}",
0200                                  "P_{z}^{b quark} (GeV)");
0201   hBottomEtaPt = dqm.book2dHisto("TTbar_BottomEtaPt",
0202                                  " quark transveral   momentum vs pseudorapidity",
0203                                  200,
0204                                  -5.,
0205                                  5.,
0206                                  1000,
0207                                  0.,
0208                                  1000.,
0209                                  "#eta_{b quark}",
0210                                  "P_{t}^{b quark} (GeV)");
0211   hBottomYPz = dqm.book2dHisto("TTbar_BottomYPz",
0212                                "b quark longitudinal momentum vs rapidity",
0213                                200,
0214                                -5.,
0215                                5.,
0216                                200,
0217                                -100.,
0218                                100.,
0219                                "Y_{b quark}",
0220                                "P_{z}^{b quark} (GeV)");
0221   hBottomMassPz = dqm.book2dHisto("TTbar_BottomMassPz",
0222                                   "b quark longitudinal momentum vs mass",
0223                                   50,
0224                                   0.,
0225                                   5.,
0226                                   200,
0227                                   -100.,
0228                                   100.,
0229                                   "M_{b quark} (GeV)",
0230                                   "P_{z}^{b quark} (GeV)");
0231   hBottomMassEta = dqm.book2dHisto("TTbar_BottomMassEta",
0232                                    "b quark pseudorapidity vs mass",
0233                                    50,
0234                                    0.,
0235                                    5.,
0236                                    200,
0237                                    -5.,
0238                                    5.,
0239                                    "M_{b quark} (GeV)",
0240                                    "#eta_{b quark}");
0241   hBottomMassY = dqm.book2dHisto(
0242       "TTbar_BottomMassY", "b quark rapidity vs mass", 50, 0., 5., 200, -5., 5., "M_{b quark} (GeV)", "Y_{b quark}");
0243   hBottomMassDeltaY = dqm.book2dHisto("TTbar_BottomMassDeltaY",
0244                                       "b quark pseudorapidity - rapidity vs mass",
0245                                       50,
0246                                       0.,
0247                                       50.,
0248                                       2000,
0249                                       -5.,
0250                                       5.,
0251                                       "M_{b quark} (GeV)",
0252                                       "Y_{b quark}");
0253 
0254   hWplusPz = dqm.book1dHisto(
0255       "TTbar_WplusPz", "W+ boson longitudinal momentum", 200, -100., 100., "P_{z}^{W+} (GeV)", "Number of Events");
0256   hWminPz = dqm.book1dHisto(
0257       "TTbar_WminPz", "W- boson longitudinal momentum", 200, -100., 100., "P_{z}^{W-} (GeV)", "Number of Events");
0258 }