Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:18:39 

0001 /** \class HLTMuonRateAnalyzerWithWeight
0002  *  Get L1/HLT efficiency/rate plots
0003  *
0004  *  \author J. Alcaraz
0005  */
0006 
0007 #include "HLTrigger/Muon/test/HLTMuonRateAnalyzerWithWeight.h"
0008 
0009 // Collaborating Class Header
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "DataFormats/Common/interface/Handle.h"
0015 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0016 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0017 
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0020 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0021 
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 
0024 #include "TFile.h"
0025 #include "TH1F.h"
0026 
0027 using namespace std;
0028 using namespace edm;
0029 using namespace reco;
0030 using namespace trigger;
0031 using namespace l1extra;
0032 
0033 /// Constructor
0034 HLTMuonRateAnalyzerWithWeight::HLTMuonRateAnalyzerWithWeight(const ParameterSet& pset) {
0035   theGenLabel = pset.getUntrackedParameter<InputTag>("GenLabel");
0036   theL1CollectionLabel = pset.getUntrackedParameter<InputTag>("L1CollectionLabel");
0037   theHLTCollectionLabels = pset.getUntrackedParameter<std::vector<InputTag> >("HLTCollectionLabels");
0038   theGenToken = consumes<edm::HepMCProduct>(theGenLabel);
0039   theL1CollectionToken = consumes<trigger::TriggerFilterObjectWithRefs>(theL1CollectionLabel);
0040   for (auto& theHLTCollectionLabel : theHLTCollectionLabels) {
0041     theHLTCollectionTokens.push_back(consumes<trigger::TriggerFilterObjectWithRefs>(theHLTCollectionLabel));
0042   }
0043   theL1ReferenceThreshold = pset.getUntrackedParameter<double>("L1ReferenceThreshold");
0044   theNSigmas = pset.getUntrackedParameter<std::vector<double> >("NSigmas90");
0045 
0046   theNumberOfObjects = pset.getUntrackedParameter<unsigned int>("NumberOfObjects");
0047 
0048   // Convert it already into /nb/s)
0049   theLuminosity = pset.getUntrackedParameter<double>("Luminosity") * 1.e-33;
0050   theIntegratedLumi = pset.getParameter<double>("IntLumi");
0051   type = pset.getParameter<unsigned int>("Type");
0052 
0053   thePtMin = pset.getUntrackedParameter<double>("PtMin");
0054   thePtMax = pset.getUntrackedParameter<double>("PtMax");
0055   theNbins = pset.getUntrackedParameter<unsigned int>("Nbins");
0056 
0057   theRootFileName = pset.getUntrackedParameter<string>("RootFileName");
0058   theNumberOfBCEvents = 0.;
0059   theNumberOfLightEvents = 0.;
0060 }
0061 
0062 /// Destructor
0063 HLTMuonRateAnalyzerWithWeight::~HLTMuonRateAnalyzerWithWeight() = default;
0064 
0065 void HLTMuonRateAnalyzerWithWeight::beginJob() {
0066   // Create the root file
0067   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0068   theFile->cd();
0069   hNumEvents = new TH1F("NumEvents", "Number of Events analyzed", 2, -0.5, 1.5);
0070 
0071   char chname[256];
0072   char chtitle[256];
0073   snprintf(chname, 255, "Lighteff_%s", theL1CollectionLabel.encode().c_str());
0074   snprintf(chtitle,
0075            255,
0076            "Light Quark events Efficiency (%%) vs L1 Pt threshold (GeV), label=%s",
0077            theL1CollectionLabel.encode().c_str());
0078   hLightL1eff = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0079   hLightL1eff->Sumw2();
0080   snprintf(chname, 255, "Lightrate_%s", theL1CollectionLabel.encode().c_str());
0081   snprintf(chtitle,
0082            255,
0083            "Light Quark events Rate (Hz) vs L1 Pt threshold (GeV), label=%s, L=%.2E (cm^{-2} s^{-1})",
0084            theL1CollectionLabel.encode().c_str(),
0085            theLuminosity * 1.e33);
0086   hLightL1rate = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0087   hLightL1rate->Sumw2();
0088   snprintf(chname, 255, "BCeff_%s", theL1CollectionLabel.encode().c_str());
0089   snprintf(chtitle,
0090            255,
0091            "BC Quark events Efficiency (%%) vs L1 Pt threshold (GeV), label=%s",
0092            theL1CollectionLabel.encode().c_str());
0093   hBCL1eff = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0094   hBCL1eff->Sumw2();
0095   snprintf(chname, 255, "BCrate_%s", theL1CollectionLabel.encode().c_str());
0096   snprintf(chtitle,
0097            255,
0098            "BC Quark events Rate (Hz) vs L1 Pt threshold (GeV), label=%s, L=%.2E (cm^{-2} s^{-1})",
0099            theL1CollectionLabel.encode().c_str(),
0100            theLuminosity * 1.e33);
0101   hBCL1rate = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0102   hBCL1rate->Sumw2();
0103 
0104   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0105     snprintf(chname, 255, "Lighteff_%s", theHLTCollectionLabels[i].encode().c_str());
0106     snprintf(chtitle,
0107              255,
0108              "Light Quark events Efficiency (%%) vs HLT Pt threshold (GeV), label=%s",
0109              theHLTCollectionLabels[i].encode().c_str());
0110     hLightHLTeff.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0111     hLightHLTeff[i]->Sumw2();
0112     snprintf(chname, 255, "Light rate_%s", theHLTCollectionLabels[i].encode().c_str());
0113     snprintf(chtitle,
0114              255,
0115              "Light Quark events Rate (Hz) vs HLT Pt threshold (GeV), label=%s, L=%.2E (cm^{-2} s^{-1})",
0116              theHLTCollectionLabels[i].encode().c_str(),
0117              theLuminosity * 1.e33);
0118     hLightHLTrate.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0119     hLightHLTrate[i]->Sumw2();
0120     snprintf(chname, 255, "BCeff_%s", theHLTCollectionLabels[i].encode().c_str());
0121     snprintf(chtitle,
0122              255,
0123              "BC Quark events Efficiency (%%) vs HLT Pt threshold (GeV), label=%s",
0124              theHLTCollectionLabels[i].encode().c_str());
0125     hBCHLTeff.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0126     hBCHLTeff[i]->Sumw2();
0127     snprintf(chname, 255, "BC rate_%s", theHLTCollectionLabels[i].encode().c_str());
0128     snprintf(chtitle,
0129              255,
0130              "BC Quark events Rate (Hz) vs HLT Pt threshold (GeV), label=%s, L=%.2E (cm^{-2} s^{-1})",
0131              theHLTCollectionLabels[i].encode().c_str(),
0132              theLuminosity * 1.e33);
0133     hBCHLTrate.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0134     hBCHLTrate[i]->Sumw2();
0135   }
0136 }
0137 
0138 void HLTMuonRateAnalyzerWithWeight::endJob() {
0139   //std::cout << "in endjob"<<endl;
0140   theFile->cd();
0141   //std::cout << "in file"<<endl;
0142   if (theNumberOfBCEvents == 0 && theNumberOfLightEvents == 0) {
0143     theFile->Close();
0144     return;
0145   }
0146   //std::cout << "we have events"<<endl;
0147 
0148   // L1 operations
0149   for (unsigned int k = 0; k <= theNbins + 1; k++) {
0150     if (theNumberOfLightEvents != 0) {
0151       double this_eff = hLightL1eff->GetBinContent(k) / theNumberOfLightEvents;
0152       double this_eff_error = hLightL1eff->GetBinError(k) / theNumberOfLightEvents * sqrt(1 - this_eff);
0153       hLightL1eff->SetBinContent(k, 100 * this_eff);
0154       hLightL1eff->SetBinError(k, 100 * this_eff_error);
0155       double this_rate = theLuminosity * this_eff * theNumberOfLightEvents / theIntegratedLumi;
0156       double this_rate_error = theLuminosity * this_eff_error * theNumberOfLightEvents / theIntegratedLumi;
0157       hLightL1rate->SetBinContent(k, this_rate);
0158       hLightL1rate->SetBinError(k, this_rate_error);
0159     }
0160     if (theNumberOfBCEvents != 0) {
0161       double this_eff = hBCL1eff->GetBinContent(k) / theNumberOfBCEvents;
0162       double this_eff_error = hBCL1eff->GetBinError(k) / theNumberOfBCEvents * sqrt(1 - this_eff);
0163       hBCL1eff->SetBinContent(k, 100 * this_eff);
0164       hBCL1eff->SetBinError(k, 100 * this_eff_error);
0165       double this_rate = theLuminosity * this_eff * theNumberOfBCEvents / theIntegratedLumi;
0166       double this_rate_error = theLuminosity * this_eff_error * theNumberOfBCEvents / theIntegratedLumi;
0167       hBCL1rate->SetBinContent(k, this_rate);
0168       hBCL1rate->SetBinError(k, this_rate_error);
0169     }
0170   }
0171 
0172   // HLT operations
0173   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0174     for (unsigned int k = 0; k <= theNbins + 1; k++) {
0175       // Hope that this will be essentially OK for weighted samples
0176       // It should be strictly OK in a binomial scheme when weights = 1
0177       if (theNumberOfLightEvents != 0) {
0178         double this_eff = hLightHLTeff[i]->GetBinContent(k) / theNumberOfLightEvents;
0179         double this_eff_error = hLightHLTeff[i]->GetBinError(k) / theNumberOfLightEvents;
0180         hLightHLTeff[i]->SetBinContent(k, 100 * this_eff);
0181         hLightHLTeff[i]->SetBinError(k, 100 * this_eff_error);
0182         double this_rate = theLuminosity * this_eff * theNumberOfLightEvents / theIntegratedLumi;
0183         double this_rate_error = theLuminosity * this_eff_error * theNumberOfLightEvents / theIntegratedLumi;
0184         hLightHLTrate[i]->SetBinContent(k, this_rate);
0185         hLightHLTrate[i]->SetBinError(k, this_rate_error);
0186       }
0187       if (theNumberOfBCEvents != 0) {
0188         double this_eff = hBCHLTeff[i]->GetBinContent(k) / theNumberOfBCEvents;
0189         double this_eff_error = hBCHLTeff[i]->GetBinError(k) / theNumberOfBCEvents;
0190         hBCHLTeff[i]->SetBinContent(k, 100 * this_eff);
0191         hBCHLTeff[i]->SetBinError(k, 100 * this_eff_error);
0192         double this_rate = theLuminosity * this_eff * theNumberOfBCEvents / theIntegratedLumi;
0193         double this_rate_error = theLuminosity * this_eff_error * theNumberOfBCEvents / theIntegratedLumi;
0194         hBCHLTrate[i]->SetBinContent(k, this_rate);
0195         hBCHLTrate[i]->SetBinError(k, this_rate_error);
0196       }
0197     }
0198   }
0199 
0200   // Write the histos to file
0201   hNumEvents->Write();
0202   hLightL1eff->Write();
0203   hLightL1rate->Write();
0204   hBCL1eff->Write();
0205   hBCL1rate->Write();
0206   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0207     hLightHLTeff[i]->Write();
0208     hLightHLTrate[i]->Write();
0209     hBCHLTeff[i]->Write();
0210     hBCHLTrate[i]->Write();
0211   }
0212   theFile->Close();
0213 }
0214 
0215 void HLTMuonRateAnalyzerWithWeight::analyze(const Event& event, const EventSetup& eventSetup) {
0216   theFile->cd();
0217 
0218   // Get the HepMC product
0219   double this_event_weight = 1.;
0220   bool bcevent = false;
0221   try {
0222     Handle<HepMCProduct> genProduct;
0223     event.getByToken(theGenToken, genProduct);
0224     const HepMC::GenEvent* evt = genProduct->GetEvent();
0225     HepMC::WeightContainer weights = evt->weights();
0226     bcevent = isbc(*evt);
0227     hNumEvents->Fill(1. * bcevent);
0228     if (weights.size() > 0)
0229       this_event_weight = weights[0];
0230     if (type == 3)
0231       this_event_weight *= parentWeight(*evt);
0232     LogInfo("HLTMuonRateAnalyzerWithWeight") << " This event weight is " << this_event_weight;
0233   } catch (...) {
0234     LogWarning("HLTMuonRateAnalyzerWithWeight") << " NO HepMCProduct found!!!!!!!!!!!!!!!";
0235     LogWarning("HLTMuonRateAnalyzerWithWeight") << " SETTING EVENT WEIGHT TO 1";
0236   }
0237   if (bcevent)
0238     theNumberOfBCEvents += this_event_weight;
0239   else
0240     theNumberOfLightEvents += this_event_weight;
0241   // Get the L1 collection
0242   Handle<TriggerFilterObjectWithRefs> l1cands;
0243   event.getByToken(theL1CollectionToken, l1cands);
0244   if (l1cands.failedToGet()) {
0245     LogInfo("HLTMuonRateAnalyzerWithWeight") << " No L1 collection";
0246     // Do nothing
0247     return;
0248   }
0249 
0250   // Get the HLT collections
0251   std::vector<Handle<TriggerFilterObjectWithRefs> > hltcands(theHLTCollectionLabels.size());
0252 
0253   unsigned int modules_in_this_event = 0;
0254   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0255     event.getByToken(theHLTCollectionTokens[i], hltcands[i]);
0256     if (hltcands[i].failedToGet()) {
0257       LogInfo("HLTMuonRateAnalyzerWithWeight") << " No " << theHLTCollectionLabels[i];
0258       break;
0259     }
0260     modules_in_this_event++;
0261   }
0262 
0263   // Fix L1 thresholds to obtain HLT plots
0264   unsigned int nL1FoundRef = 0;
0265   double epsilon = 0.001;
0266   vector<L1MuonParticleRef> l1mu;
0267   l1cands->getObjects(TriggerL1Mu, l1mu);
0268   for (auto& k : l1mu) {
0269     L1MuonParticleRef candref = L1MuonParticleRef(k);
0270     // L1 PTs are "quantized" due to LUTs.
0271     // Their meaning: true_pt > ptLUT more than 90% pof the times
0272     double ptLUT = candref->pt();
0273     // Add "epsilon" to avoid rounding errors when ptLUT==L1Threshold
0274     if (ptLUT + epsilon > theL1ReferenceThreshold)
0275       nL1FoundRef++;
0276   }
0277 
0278   for (unsigned int j = 0; j < theNbins; j++) {
0279     double ptcut = thePtMin + j * (thePtMax - thePtMin) / theNbins;
0280 
0281     // L1 filling
0282     unsigned int nFound = 0;
0283     for (auto& k : l1mu) {
0284       L1MuonParticleRef candref = L1MuonParticleRef(k);
0285       double pt = candref->pt();
0286       if (pt > ptcut)
0287         nFound++;
0288     }
0289     if (nFound >= theNumberOfObjects) {
0290       if (bcevent)
0291         hBCL1eff->Fill(ptcut, this_event_weight);
0292       else
0293         hLightL1eff->Fill(ptcut, this_event_weight);
0294     }
0295     // Stop here if L1 reference cuts were not satisfied
0296     if (nL1FoundRef < theNumberOfObjects)
0297       continue;
0298 
0299     // HLT filling
0300     for (unsigned int i = 0; i < modules_in_this_event; i++) {
0301       unsigned nFound = 0;
0302       vector<RecoChargedCandidateRef> vref;
0303       hltcands[i]->getObjects(TriggerMuon, vref);
0304       for (auto& k : vref) {
0305         RecoChargedCandidateRef candref = RecoChargedCandidateRef(k);
0306         TrackRef tk = candref->get<TrackRef>();
0307         double pt = tk->pt();
0308         double err0 = tk->error(0);
0309         double abspar0 = fabs(tk->parameter(0));
0310         // convert to 90% efficiency threshold
0311         if (abspar0 > 0)
0312           pt += theNSigmas[i] * err0 / abspar0 * pt;
0313         if (pt > ptcut)
0314           nFound++;
0315       }
0316       if (nFound >= theNumberOfObjects) {
0317         if (bcevent)
0318           hBCHLTeff[i]->Fill(ptcut, this_event_weight);
0319         else
0320           hLightHLTeff[i]->Fill(ptcut, this_event_weight);
0321       } else {
0322         break;
0323       }
0324     }
0325   }
0326 }
0327 
0328 bool HLTMuonRateAnalyzerWithWeight::isbc(HepMC::GenEvent const& Gevt) {
0329   bool mybc = false;
0330   int npart = 0;
0331   int nb = 0;
0332   int nc = 0;
0333   for (HepMC::GenEvent::particle_const_iterator particle = Gevt.particles_begin(); particle != Gevt.particles_end();
0334        ++particle) {
0335     ++npart;
0336     int id = abs((*particle)->pdg_id());
0337     //  int status=(*particle)->status();
0338     if (id == 5 || id == 4) {
0339       if (npart == 6 || npart == 7) {
0340         mybc = true;
0341         break;
0342       } else {
0343         HepMC::GenVertex* parent = (*particle)->production_vertex();
0344         for (auto ic = parent->particles_in_const_begin(); ic != parent->particles_in_const_end(); ic++) {
0345           int pid = (*ic)->pdg_id();
0346           if (pid == 21 && id == 5)
0347             nb++;
0348           else if (pid == 21 && id == 4)
0349             nc++;
0350         }
0351       }
0352     }
0353   }
0354   if (nb > 1 || nc > 1)
0355     mybc = true;
0356   return mybc;
0357 }
0358 
0359 double HLTMuonRateAnalyzerWithWeight::parentWeight(HepMC::GenEvent const& Gevt) {
0360   double AdditionalWeight = 1.;
0361   if (type != 3)
0362     return AdditionalWeight;
0363   for (HepMC::GenEvent::particle_const_iterator particle = Gevt.particles_begin(); particle != Gevt.particles_end();
0364        ++particle) {
0365     int id = abs((*particle)->pdg_id());
0366     double pt = (*particle)->momentum().perp();
0367     if (id == 13 && pt > 10) {
0368       HepMC::GenVertex* parent = (*particle)->production_vertex();
0369       for (auto ic = parent->particles_in_const_begin(); ic != parent->particles_in_const_end(); ic++) {
0370         int apid = abs((*ic)->pdg_id());
0371         LogInfo("HLTMuonRateAnalyzerWithWeight") << " Absolute parent id is " << apid;
0372         if (apid > 10000)
0373           apid = apid - (apid / 10000) * 10000;
0374         if (apid > 1000)
0375           apid /= 1000;
0376         if (apid > 100 && apid != 130)
0377           apid /= 100;
0378         LogInfo("HLTMuonRateAnalyzerWithWeight") << " It will be treated as " << apid;
0379         if (apid == 5)
0380           AdditionalWeight = 1. / 8.4;  //b mesons
0381         else if (apid == 4)
0382           AdditionalWeight = 1. / 6.0;  //c mesons
0383         else if (apid == 15)
0384           AdditionalWeight = 1. / 8.7;  //taus
0385         else if (apid == 3 || apid == 130)
0386           AdditionalWeight = 1. / 7.3;  //s-mesons
0387         else if (apid == 2)
0388           AdditionalWeight = 1. / 0.8;  //pions
0389       }
0390     }
0391   }
0392   return AdditionalWeight;
0393 }
0394 
0395 DEFINE_FWK_MODULE(HLTMuonRateAnalyzerWithWeight);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!