Muon/test/MuEnrichRenormalizer.cc

0001 #include <iostream>
0002
0003 #include "HLTrigger/Muon/test/MuEnrichRenormalizer.h"
0004
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "FWCore/Framework/interface/MakerMacros.h"
0009
0010 using namespace edm;
0011 using namespace std;
0012
0013 MuEnrichRenormalizer::MuEnrichRenormalizer(const ParameterSet& pset) : type(pset.getParameter<int>("type")) {
0014   theGenToken = consumes<edm::HepMCProduct>(edm::InputTag("VtxSmeared"));
0015   anaIntlumi = 0;
0016   anaLight = 0;
0017   anaBC = 0;
0018   rwbc = 1.;
0019   rwlight = 1.;
0020   if (type == 1) {
0021     genIntlumi = 0.02465;  // nb-1 generated events
0022     genLight = 426000;     // total generated light quark events
0023     genBC = 192085;        // total generated bc quark events
0024   } else if (type == 2) {
0025     genIntlumi = 0.2465;  // nb-1 generated events
0026     genLight = 168300;    // total generated light quark events
0027     genBC = 257561;       // total generated bc quark events
0028   } else if (type == 3) {
0029     genIntlumi = 0.242;
0030     genLight = 10800;
0031     genBC = 141596;
0032   }
0033 }
0034
0035 void MuEnrichRenormalizer::beginJob() {}
0036
0037 void MuEnrichRenormalizer::analyze(const Event& e, const EventSetup&) {
0038   Handle<HepMCProduct> EvtHandle;
0039   e.getByToken(theGenToken, EvtHandle);
0040   const HepMC::GenEvent* Gevt = EvtHandle->GetEvent();
0041   bool mybc = false;
0042   int npart = 0;
0043   int nb = 0;
0044   int nc = 0;
0045   if (Gevt != nullptr) {
0046     for (HepMC::GenEvent::particle_const_iterator particle = Gevt->particles_begin(); particle != Gevt->particles_end();
0047          ++particle) {
0048       ++npart;
0049       int id = abs((*particle)->pdg_id());
0050       //int status=(*particle)->status();
0051       if (id == 5 || id == 4) {
0052         if (npart == 6 || npart == 7) {
0053           mybc = true;
0054           break;
0055         } else {
0056           HepMC::GenVertex* parent = (*particle)->production_vertex();
0057           for (auto ic = parent->particles_in_const_begin(); ic != parent->particles_in_const_end(); ic++) {
0058             int pid = (*ic)->pdg_id();
0059             if (pid == 21 && id == 5)
0060               nb++;
0061             else if (pid == 21 && id == 4)
0062               nc++;
0063           }
0064         }
0065       }
0066     }
0067   }
0068   if (nb > 1 || nc > 1)
0069     mybc = true;
0070   if (mybc)
0071     ++anaBC;
0072   else
0073     ++anaLight;
0074   return;
0075 }
0076
0077 void MuEnrichRenormalizer::endJob() {
0078   double rLight = double(anaLight) / genLight;
0079   double rBC = double(anaBC) / genBC;
0080   edm::LogVerbatim("MuEnrichRenormalizer") << "Generated uds events=" << genLight;
0081   edm::LogVerbatim("MuEnrichRenormalizer") << "Generated bc events=" << genBC;
0082   edm::LogVerbatim("MuEnrichRenormalizer") << "Analyzed uds events=" << anaLight;
0083   edm::LogVerbatim("MuEnrichRenormalizer") << "Analyzed bc events=" << anaBC;
0084   edm::LogVerbatim("MuEnrichRenormalizer") << "Ratio of analyzed to simulated uds=" << rLight;
0085   edm::LogVerbatim("MuEnrichRenormalizer") << "Ratio of analyzed to simulated bc =" << rBC;
0086   if (rBC < rLight) {
0087     rwlight = rBC / rLight;
0088     anaIntlumi = genIntlumi * rBC;
0089   } else {
0090     rwbc = rLight / rBC;
0091     anaIntlumi = genIntlumi * rLight;
0092   }
0093   edm::LogVerbatim("MuEnrichRenormalizer") << "rwbc=" << rwbc << " rwlight=" << rwlight;
0094   edm::LogVerbatim("MuEnrichRenormalizer") << "Corresponding Lumi =" << anaIntlumi << " nb-1";
0095   return;
0096 }
0097
0098 DEFINE_FWK_MODULE(MuEnrichRenormalizer);