HepMCCandAlgos/test/HZZ4muAnalyzer.cc

0001 #include <iostream>
0002
0003 #include "PhysicsTools/HepMCCandAlgos/test/HZZ4muAnalyzer.h"
0004
0005 // essentials !!!
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009
0010 #include "TFile.h"
0011 #include "TH1.h"
0012
0013 using namespace edm;
0014 using namespace std;
0015
0016 HZZ4muAnalyzer::HZZ4muAnalyzer(const ParameterSet& pset)
0017     : fToken(consumes<HepMCProduct>(InputTag("source"))),
0018       fOutputFileName(pset.getUntrackedParameter<string>("HistOutFile", std::string("TestHiggsMass.root"))),
0019       fOutputFile(0),
0020       fHist2muMass(0),
0021       fHist4muMass(0),
0022       fHistZZMass(0) {}
0023
0024 void HZZ4muAnalyzer::beginJob() {
0025   fOutputFile = new TFile(fOutputFileName.c_str(), "RECREATE");
0026   fHist2muMass = new TH1D("Hist2muMass", "2-mu inv. mass", 100, 60., 120.);
0027   fHist4muMass = new TH1D("Hist4muMass", "4-mu inv. mass", 100, 170., 210.);
0028   fHistZZMass = new TH1D("HistZZMass", "ZZ inv. mass", 100, 170., 210.);
0029
0030   return;
0031 }
0032
0033 void HZZ4muAnalyzer::analyze(const Event& e, const EventSetup&) {
0034   Handle<HepMCProduct> EvtHandle;
0035
0036   // find initial (unsmeared, unfiltered,...) HepMCProduct
0037   //
0038   e.getByToken(fToken, EvtHandle);
0039
0040   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0041
0042   // this a pointer - and not an array/vector/...
0043   // because this example explicitely assumes
0044   // that there one and only Higgs in the record
0045   //
0046   HepMC::GenVertex* HiggsDecVtx = 0;
0047
0048   // find the 1st vertex with outgoing Higgs
0049   // and get Higgs decay vertex from there;
0050   //
0051   // in principal, one can look for the vertex
0052   // with incoming Higgs as well...
0053   //
0054   for (HepMC::GenEvent::vertex_const_iterator vit = Evt->vertices_begin(); vit != Evt->vertices_end(); vit++) {
0055     for (HepMC::GenVertex::particles_out_const_iterator pout = (*vit)->particles_out_const_begin();
0056          pout != (*vit)->particles_out_const_end();
0057          pout++) {
0058       if ((*pout)->pdg_id() == 25 && (*pout)->status() == 2) {
0059         if ((*pout)->end_vertex() != 0) {
0060           HiggsDecVtx = (*pout)->end_vertex();
0061           break;
0062         }
0063       }
0064     }
0065     if (HiggsDecVtx != 0) {
0066       break;  // break the initial loop over vertices
0067     }
0068   }
0069
0070   if (HiggsDecVtx == 0) {
0071     cout << " There is NO Higgs in this event ! " << endl;
0072     return;
0073   }
0074
0075   if (e.id().event() == 1) {
0076     cout << " " << endl;
0077     cout << " We do some example printouts in the event 1 " << endl;
0078     cout << " Higgs decay found at the vertex " << HiggsDecVtx->barcode() << " (barcode)" << endl;
0079
0080     vector<HepMC::GenParticle*> HiggsChildren;
0081
0082     for (HepMC::GenVertex::particles_out_const_iterator H0in = HiggsDecVtx->particles_out_const_begin();
0083          H0in != HiggsDecVtx->particles_out_const_end();
0084          H0in++) {
0085       HiggsChildren.push_back(*H0in);
0086     }
0087     cout << " Number of Higgs (immediate) children = " << HiggsChildren.size() << endl;
0088     for (unsigned int ic = 0; ic < HiggsChildren.size(); ic++) {
0089       HiggsChildren[ic]->print();
0090     }
0091   }
0092
0093   // select and store stable descendants of the Higgs
0094   //
0095   vector<HepMC::GenParticle*> StableHiggsDesc;
0096
0097   if (e.id().event() == 1)
0098     cout << " Now let us list all descendents of the Higgs" << endl;
0099   for (HepMC::GenVertex::particle_iterator des = HiggsDecVtx->particles_begin(HepMC::descendants);
0100        des != HiggsDecVtx->particles_end(HepMC::descendants);
0101        des++) {
0102     if (e.id().event() == 1)
0103       (*des)->print();
0104     if ((*des)->status() == 1)
0105       StableHiggsDesc.push_back(*des);
0106   }
0107
0108   HepMC::FourVector Mom2part;
0109   double XMass2part = 0.;
0110   double XMass4part = 0.;
0111   double XMass2pairs = 0.;
0112   vector<HepMC::FourVector> Mom2partCont;
0113
0114   // browse the array of stable descendants
0115   // and do 2-mu inv.mass
0116   //
0117   for (unsigned int i = 0; i < StableHiggsDesc.size(); i++) {
0118     // skip other than mu
0119     //
0120     if (std::abs(StableHiggsDesc[i]->pdg_id()) != 13)
0121       continue;
0122
0123     for (unsigned int j = i + 1; j < StableHiggsDesc.size(); j++) {
0124       // skip other than mu
0125       //
0126       if (std::abs(StableHiggsDesc[j]->pdg_id()) != 13)
0127         continue;
0128       //
0129       // skip same charge combo's
0130       //
0131       if ((StableHiggsDesc[i]->pdg_id() * StableHiggsDesc[j]->pdg_id()) > 0)
0132         continue;
0133       //
0134       // OK, opposite charges, do the job
0135       //
0136       Mom2part = HepMC::FourVector((StableHiggsDesc[i]->momentum().px() + StableHiggsDesc[j]->momentum().px()),
0137                                    (StableHiggsDesc[i]->momentum().py() + StableHiggsDesc[j]->momentum().py()),
0138                                    (StableHiggsDesc[i]->momentum().pz() + StableHiggsDesc[j]->momentum().pz()),
0139                                    (StableHiggsDesc[i]->momentum().e() + StableHiggsDesc[j]->momentum().e()));
0140
0141       XMass2part = Mom2part.m();
0142       fHist2muMass->Fill(XMass2part);
0143       //cout << " counters : " << StableHiggsDesc[i]->barcode() << " "
0144       //                       << StableHiggsDesc[j]->barcode()
0145       //            << " -> 2-part mass = " << XMass2part << endl ;
0146       //
0147       // store if 2-part. inv. mass fits into (roughly) Z-mass interval
0148       //
0149       if (XMass2part > 80. && XMass2part < 100.) {
0150         Mom2partCont.push_back(Mom2part);
0151       }
0152     }
0153   }
0154
0155   // make 4-part inv.mass
0156   //
0157   double px4, py4, pz4, e4;
0158   px4 = py4 = pz4 = e4 = 0.;
0159   if (StableHiggsDesc.size() == 4) {
0160     for (unsigned int i = 0; i < StableHiggsDesc.size(); i++) {
0161       px4 += StableHiggsDesc[i]->momentum().px();
0162       py4 += StableHiggsDesc[i]->momentum().py();
0163       pz4 += StableHiggsDesc[i]->momentum().pz();
0164       e4 += StableHiggsDesc[i]->momentum().e();
0165     }
0166     XMass4part = HepMC::FourVector(px4, py4, pz4, e4).m();
0167     fHist4muMass->Fill(XMass4part);
0168   }
0169   //cout << " 4-part inv. mass = " << XMass4part << endl ;
0170
0171   // make 2-pairs (ZZ) inv.mass
0172   //
0173   //cout << " selected Z-candidates in this event : " << Mom2partCont.size() << endl ;
0174   for (unsigned int i = 0; i < Mom2partCont.size(); i++) {
0175     for (unsigned int j = i + 1; j < Mom2partCont.size(); j++) {
0176       // Mom2pairs = Mom2partCont[i] + Mom2partCont[j] ;
0177       XMass2pairs = HepMC::FourVector((Mom2partCont[i].px() + Mom2partCont[j].px()),
0178                                       (Mom2partCont[i].py() + Mom2partCont[j].py()),
0179                                       (Mom2partCont[i].pz() + Mom2partCont[j].pz()),
0180                                       (Mom2partCont[i].e() + Mom2partCont[j].e()))
0181                         .m();
0182       fHistZZMass->Fill(XMass2pairs);
0183       //cout << " 2-pairs (ZZ) inv. mass = " << XMass2pairs << endl ;
0184     }
0185   }
0186
0187   return;
0188 }
0189
0190 void HZZ4muAnalyzer::endJob() {
0191   fOutputFile->Write();
0192   fOutputFile->Close();
0193
0194   return;
0195 }
0196
0197 typedef HZZ4muAnalyzer HZZ4muTestAnalyzer;
0198 DEFINE_FWK_MODULE(HZZ4muTestAnalyzer);