Pythia6Interface/test/HZZ4muAnalyzer.cc

0001 #include <iostream>
0002
0003 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0004 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0005 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0006
0007 // essentials !!!
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "DataFormats/Common/interface/Handle.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0014 #include "TH1.h"
0015
0016 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0017
0018 class HZZ4muAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0019 public:
0020   //
0021   explicit HZZ4muAnalyzer(const edm::ParameterSet&);
0022   ~HZZ4muAnalyzer() = default;  // no need to delete ROOT stuff
0023                                 // as it'll be deleted upon closing TFile
0024
0025   void analyze(const edm::Event&, const edm::EventSetup&) override;
0026   void beginJob() override;
0027
0028 private:
0029   const edm::EDGetTokenT<GenEventInfoProduct> tokenGenInfo_;
0030   const edm::EDGetTokenT<edm::HepMCProduct> tokenHepMC_;
0031   TH1D* fHist2muMass;
0032   TH1D* fHist4muMass;
0033   TH1D* fHistZZMass;
0034 };
0035
0036 HZZ4muAnalyzer::HZZ4muAnalyzer(const edm::ParameterSet& pset)
0037     : tokenGenInfo_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0038       tokenHepMC_(consumes<edm::HepMCProduct>(edm::InputTag("VtxSmeared"))),
0039       fHist2muMass(nullptr),
0040       fHist4muMass(nullptr),
0041       fHistZZMass(nullptr) {
0042   // actually, pset is NOT in use - we keep it here just for illustratory putposes
0043   usesResource(TFileService::kSharedResource);
0044 }
0045
0046 void HZZ4muAnalyzer::beginJob() {
0047   edm::Service<TFileService> fs;
0048   fHist2muMass = fs->make<TH1D>("Hist2muMass", "2-mu inv. mass", 100, 60., 120.);
0049   fHist4muMass = fs->make<TH1D>("Hist4muMass", "4-mu inv. mass", 100, 170., 210.);
0050   fHistZZMass = fs->make<TH1D>("HistZZMass", "ZZ inv. mass", 100, 170., 210.);
0051 }
0052
0053 void HZZ4muAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
0054   // here's an example of accessing GenEventInfoProduct
0055   const edm::Handle<GenEventInfoProduct>& GenInfoHandle = e.getHandle(tokenGenInfo_);
0056   double qScale = GenInfoHandle->qScale();
0057   double pthat = (GenInfoHandle->hasBinningValues() ? (GenInfoHandle->binningValues())[0] : 0.0);
0058   std::cout << " qScale = " << qScale << " pthat = " << pthat << std::endl;
0059   //
0060   // this (commented out) code below just exemplifies how to access certain info
0061   //
0062   //double evt_weight1 = GenInfoHandle->weights()[0]; // this is "stanrd Py6 evt weight;
0063   // corresponds to PYINT1/VINT(97)
0064   //double evt_weight2 = GenInfoHandle->weights()[1]; // in case you run in CSA mode or otherwise
0065   // use PYEVWT routine, this will be weight
0066   // as returned by PYEVWT, i.e. PYINT1/VINT(99)
0067   //std::cout << " evt_weight1 = " << evt_weight1 << std::endl;
0068   //std::cout << " evt_weight2 = " << evt_weight2 << std::endl;
0069   //double weight = GenInfoHandle->weight();
0070   //std::cout << " as returned by the weight() method, integrated event weight = " << weight << std::endl;
0071
0072   // here's an example of accessing particles in the event record (HepMCProduct)
0073   //
0074   // find initial (unsmeared, unfiltered,...) HepMCProduct
0075   const edm::Handle<edm::HepMCProduct>& EvtHandle = e.getHandle(tokenHepMC_);
0076
0077   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0078
0079   // this a pointer - and not an array/vector/...
0080   // because this example explicitely assumes
0081   // that there one and only Higgs in the record
0082   //
0083   HepMC::GenVertex* HiggsDecVtx = 0;
0084
0085   // find the 1st vertex with outgoing Higgs
0086   // and get Higgs decay vertex from there;
0087   //
0088   // in principal, one can look for the vertex
0089   // with incoming Higgs as well...
0090   //
0091   for (HepMC::GenEvent::vertex_const_iterator vit = Evt->vertices_begin(); vit != Evt->vertices_end(); vit++) {
0092     for (HepMC::GenVertex::particles_out_const_iterator pout = (*vit)->particles_out_const_begin();
0093          pout != (*vit)->particles_out_const_end();
0094          pout++) {
0095       if ((*pout)->pdg_id() == 25) {
0096         if ((*pout)->end_vertex() != 0) {
0097           HiggsDecVtx = (*pout)->end_vertex();
0098           break;
0099         }
0100       }
0101     }
0102     if (HiggsDecVtx != 0) {
0103       break;  // break the initial loop over vertices
0104     }
0105   }
0106
0107   if (HiggsDecVtx == 0) {
0108     std::cout << " There is NO Higgs in this event ! " << std::endl;
0109     return;
0110   }
0111
0112   if (e.id().event() == 1) {
0113     std::cout << " " << std::endl;
0114     std::cout << " We do some example printouts in the event 1 " << std::endl;
0115     std::cout << " Higgs decay found at the vertex " << HiggsDecVtx->barcode() << " (barcode)" << std::endl;
0116
0117     std::vector<HepMC::GenParticle*> HiggsChildren;
0118
0119     for (HepMC::GenVertex::particles_out_const_iterator H0in = HiggsDecVtx->particles_out_const_begin();
0120          H0in != HiggsDecVtx->particles_out_const_end();
0121          H0in++) {
0122       HiggsChildren.push_back(*H0in);
0123     }
0124     std::cout << " Number of Higgs (immediate) children = " << HiggsChildren.size() << std::endl;
0125     for (unsigned int ic = 0; ic < HiggsChildren.size(); ic++) {
0126       HiggsChildren[ic]->print();
0127     }
0128   }
0129
0130   // select and store stable descendants of the Higgs
0131   //
0132   std::vector<HepMC::GenParticle*> StableHiggsDesc;
0133
0134   if (e.id().event() == 1)
0135     std::cout << " Now let us list all descendents of the Higgs" << std::endl;
0136   for (HepMC::GenVertex::particle_iterator des = HiggsDecVtx->particles_begin(HepMC::descendants);
0137        des != HiggsDecVtx->particles_end(HepMC::descendants);
0138        des++) {
0139     if (e.id().event() == 1)
0140       (*des)->print();
0141     if ((*des)->status() == 1)
0142       StableHiggsDesc.push_back(*des);
0143   }
0144
0145   HepMC::FourVector Mom2part;
0146   double XMass2part = 0.;
0147   double XMass4part = 0.;
0148   double XMass2pairs = 0.;
0149   std::vector<HepMC::FourVector> Mom2partCont;
0150
0151   // browse the array of stable descendants
0152   // and do 2-mu inv.mass
0153   //
0154   for (unsigned int i = 0; i < StableHiggsDesc.size(); i++) {
0155     // skip other than mu
0156     //
0157     if (std::abs(StableHiggsDesc[i]->pdg_id()) != 13)
0158       continue;
0159
0160     for (unsigned int j = i + 1; j < StableHiggsDesc.size(); j++) {
0161       // skip other than mu
0162       //
0163       if (std::abs(StableHiggsDesc[j]->pdg_id()) != 13)
0164         continue;
0165       //
0166       // skip same charge combo's
0167       //
0168       if ((StableHiggsDesc[i]->pdg_id() * StableHiggsDesc[j]->pdg_id()) > 0)
0169         continue;
0170       //
0171       // OK, opposite charges, do the job
0172       //
0173       Mom2part = HepMC::FourVector((StableHiggsDesc[i]->momentum().px() + StableHiggsDesc[j]->momentum().px()),
0174                                    (StableHiggsDesc[i]->momentum().py() + StableHiggsDesc[j]->momentum().py()),
0175                                    (StableHiggsDesc[i]->momentum().pz() + StableHiggsDesc[j]->momentum().pz()),
0176                                    (StableHiggsDesc[i]->momentum().e() + StableHiggsDesc[j]->momentum().e()));
0177
0178       XMass2part = Mom2part.m();
0179       fHist2muMass->Fill(XMass2part);
0180       //cout << " counters : " << StableHiggsDesc[i]->barcode() << " "
0181       //                       << StableHiggsDesc[j]->barcode()
0182       //            << " -> 2-part mass = " << XMass2part << endl ;
0183       //
0184       // store if 2-part. inv. mass fits into (roughly) Z-mass interval
0185       //
0186       if (XMass2part > 80. && XMass2part < 100.) {
0187         Mom2partCont.push_back(Mom2part);
0188       }
0189     }
0190   }
0191
0192   // make 4-part inv.mass
0193   //
0194
0195   double px4, py4, pz4, e4;
0196   px4 = py4 = pz4 = e4 = 0.;
0197   if (StableHiggsDesc.size() == 4) {
0198     for (unsigned int i = 0; i < StableHiggsDesc.size(); i++) {
0199       px4 += StableHiggsDesc[i]->momentum().px();
0200       py4 += StableHiggsDesc[i]->momentum().py();
0201       pz4 += StableHiggsDesc[i]->momentum().pz();
0202       e4 += StableHiggsDesc[i]->momentum().e();
0203     }
0204     XMass4part = HepMC::FourVector(px4, py4, pz4, e4).m();
0205     fHist4muMass->Fill(XMass4part);
0206   }
0207   //std::cout << " 4-part inv. mass = " << XMass4part << std::endl ;
0208
0209   // make 2-pairs (ZZ) inv.mass
0210   //
0211   //std::cout << " selected Z-candidates in this event : " << Mom2partCont.size() << std::endl ;
0212   for (unsigned int i = 0; i < Mom2partCont.size(); i++) {
0213     for (unsigned int j = i + 1; j < Mom2partCont.size(); j++) {
0214       // Mom2pairs = Mom2partCont[i] + Mom2partCont[j] ;
0215       XMass2pairs = HepMC::FourVector((Mom2partCont[i].px() + Mom2partCont[j].px()),
0216                                       (Mom2partCont[i].py() + Mom2partCont[j].py()),
0217                                       (Mom2partCont[i].pz() + Mom2partCont[j].pz()),
0218                                       (Mom2partCont[i].e() + Mom2partCont[j].e()))
0219                         .m();
0220       fHistZZMass->Fill(XMass2pairs);
0221       //cout << " 2-pairs (ZZ) inv. mass = " << XMass2pairs << endl ;
0222     }
0223   }
0224
0225   return;
0226 }
0227
0228 typedef HZZ4muAnalyzer HZZ4muExampleAnalyzer;
0229 DEFINE_FWK_MODULE(HZZ4muExampleAnalyzer);