Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /** \class HLTMuonRateAnalyzer
0002  *  Get L1/HLT efficiency/rate plots
0003  *
0004  *  \author J. Alcaraz
0005  */
0006 
0007 #include "HLTrigger/Muon/test/HLTMuonRateAnalyzer.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 
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/L1Trigger/interface/L1MuonParticle.h"
0018 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.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 HLTMuonRateAnalyzer::HLTMuonRateAnalyzer(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   theCrossSection = pset.getUntrackedParameter<double>("CrossSection");
0048   // Convert it already into /nb/s)
0049   theLuminosity = pset.getUntrackedParameter<double>("Luminosity") * 1.e-33;
0050 
0051   thePtMin = pset.getUntrackedParameter<double>("PtMin");
0052   thePtMax = pset.getUntrackedParameter<double>("PtMax");
0053   theNbins = pset.getUntrackedParameter<unsigned int>("Nbins");
0054 
0055   theRootFileName = pset.getUntrackedParameter<string>("RootFileName");
0056 
0057   theNumberOfEvents = 0.;
0058 }
0059 
0060 /// Destructor
0061 HLTMuonRateAnalyzer::~HLTMuonRateAnalyzer() = default;
0062 
0063 void HLTMuonRateAnalyzer::beginJob() {
0064   // Create the root file
0065   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0066   theFile->cd();
0067 
0068   char chname[256];
0069   char chtitle[256];
0070   snprintf(chname, 255, "eff_%s", theL1CollectionLabel.encode().c_str());
0071   snprintf(chtitle, 255, "Efficiency (%%) vs L1 Pt threshold (GeV), label=%s", theL1CollectionLabel.encode().c_str());
0072   hL1eff = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0073   snprintf(chname, 255, "rate_%s", theL1CollectionLabel.encode().c_str());
0074   snprintf(chtitle,
0075            255,
0076            "Rate (Hz) vs L1 Pt threshold (GeV), label=%s, L=%.2E (cm^{-2} s^{-1})",
0077            theL1CollectionLabel.encode().c_str(),
0078            theLuminosity * 1.e33);
0079   hL1rate = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0080 
0081   for (auto& theHLTCollectionLabel : theHLTCollectionLabels) {
0082     snprintf(chname, 255, "eff_%s", theHLTCollectionLabel.encode().c_str());
0083     snprintf(
0084         chtitle, 255, "Efficiency (%%) vs HLT Pt threshold (GeV), label=%s", theHLTCollectionLabel.encode().c_str());
0085     hHLTeff.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0086     snprintf(chname, 255, "rate_%s", theHLTCollectionLabel.encode().c_str());
0087     snprintf(chtitle,
0088              255,
0089              "Rate (Hz) vs HLT Pt threshold (GeV), label=%s, L=%.2E (cm^{-2} s^{-1})",
0090              theHLTCollectionLabel.encode().c_str(),
0091              theLuminosity * 1.e33);
0092     hHLTrate.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0093   }
0094 }
0095 
0096 void HLTMuonRateAnalyzer::endJob() {
0097   LogInfo("HLTMuonRateAnalyzer") << " (Weighted) number of analyzed events= " << theNumberOfEvents;
0098   theFile->cd();
0099 
0100   if (theNumberOfEvents == 0) {
0101     LogInfo("HLTMuonRateAnalyzer") << " No histograms will be written because number of events=0!!!";
0102     theFile->Close();
0103     return;
0104   }
0105 
0106   // L1 operations
0107   for (unsigned int k = 0; k <= theNbins + 1; k++) {
0108     double this_eff = hL1eff->GetBinContent(k) / theNumberOfEvents;
0109     // Hope that this will be essentially OK for weighted samples
0110     // It should be strictly OK in a binomial scheme when weights = 1
0111     double this_eff_error = hL1eff->GetBinError(k) / theNumberOfEvents * sqrt(1 - this_eff);
0112     hL1eff->SetBinContent(k, 100 * this_eff);
0113     hL1eff->SetBinError(k, 100 * this_eff_error);
0114     double this_rate = theLuminosity * theCrossSection * this_eff;
0115     double this_rate_error = theLuminosity * theCrossSection * this_eff_error;
0116     hL1rate->SetBinContent(k, this_rate);
0117     hL1rate->SetBinError(k, this_rate_error);
0118   }
0119 
0120   // HLT operations
0121   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0122     for (unsigned int k = 0; k <= theNbins + 1; k++) {
0123       // Hope that this will be essentially OK for weighted samples
0124       // It should be strictly OK in a binomial scheme when weights = 1
0125       double this_eff = hHLTeff[i]->GetBinContent(k) / theNumberOfEvents;
0126       double this_eff_error = hHLTeff[i]->GetBinError(k) / theNumberOfEvents;
0127       hHLTeff[i]->SetBinContent(k, 100 * this_eff);
0128       hHLTeff[i]->SetBinError(k, 100 * this_eff_error);
0129       double this_rate = theLuminosity * theCrossSection * this_eff;
0130       double this_rate_error = theLuminosity * theCrossSection * this_eff_error;
0131       hHLTrate[i]->SetBinContent(k, this_rate);
0132       hHLTrate[i]->SetBinError(k, this_rate_error);
0133     }
0134   }
0135 
0136   // Write the histos to file
0137   hL1eff->Write();
0138   hL1rate->Write();
0139   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0140     hHLTeff[i]->Write();
0141     hHLTrate[i]->Write();
0142   }
0143 
0144   theFile->Close();
0145 }
0146 
0147 void HLTMuonRateAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0148   theFile->cd();
0149 
0150   // Get the HepMC product
0151   double this_event_weight = 1.;
0152   try {
0153     Handle<HepMCProduct> genProduct;
0154     event.getByToken(theGenToken, genProduct);
0155     const HepMC::GenEvent* evt = genProduct->GetEvent();
0156     HepMC::WeightContainer weights = evt->weights();
0157     if (weights.size() > 0)
0158       this_event_weight = weights[0];
0159   } catch (...) {
0160     LogInfo("HLTMuonRateAnalyzer") << " NO HepMCProduct found!!!!!!!!!!!!!!!";
0161     LogInfo("HLTMuonRateAnalyzer") << " SETTING EVENT WEIGHT TO 1";
0162   }
0163   theNumberOfEvents += this_event_weight;
0164 
0165   // Get the L1 collection
0166   Handle<TriggerFilterObjectWithRefs> l1cands;
0167   event.getByToken(theL1CollectionToken, l1cands);
0168   if (l1cands.failedToGet())
0169     return;
0170 
0171   // Get the HLT collections
0172   std::vector<Handle<TriggerFilterObjectWithRefs> > hltcands(theHLTCollectionLabels.size());
0173   unsigned int modules_in_this_event = 0;
0174   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0175     event.getByToken(theHLTCollectionTokens[i], hltcands[i]);
0176     if (hltcands[i].failedToGet())
0177       break;
0178     modules_in_this_event++;
0179   }
0180 
0181   // Fix L1 thresholds to obtain HLT plots
0182   unsigned int nL1FoundRef = 0;
0183   double epsilon = 0.001;
0184   vector<L1MuonParticleRef> l1mu;
0185   l1cands->getObjects(TriggerL1Mu, l1mu);
0186   for (auto& k : l1mu) {
0187     L1MuonParticleRef candref = L1MuonParticleRef(k);
0188     // L1 PTs are "quantized" due to LUTs.
0189     // Their meaning: true_pt > ptLUT more than 90% pof the times
0190     double ptLUT = candref->pt();
0191     // Add "epsilon" to avoid rounding errors when ptLUT==L1Threshold
0192     if (ptLUT + epsilon > theL1ReferenceThreshold)
0193       nL1FoundRef++;
0194   }
0195 
0196   for (unsigned int j = 0; j < theNbins; j++) {
0197     double ptcut = thePtMin + j * (thePtMax - thePtMin) / theNbins;
0198 
0199     // L1 filling
0200     unsigned int nFound = 0;
0201     for (auto& k : l1mu) {
0202       L1MuonParticleRef candref = L1MuonParticleRef(k);
0203       double pt = candref->pt();
0204       if (pt > ptcut)
0205         nFound++;
0206     }
0207     if (nFound >= theNumberOfObjects)
0208       hL1eff->Fill(ptcut, this_event_weight);
0209 
0210     // Stop here if L1 reference cuts were not satisfied
0211     if (nL1FoundRef < theNumberOfObjects)
0212       continue;
0213 
0214     // HLT filling
0215     for (unsigned int i = 0; i < modules_in_this_event; i++) {
0216       unsigned nFound = 0;
0217       vector<RecoChargedCandidateRef> vref;
0218       hltcands[i]->getObjects(TriggerMuon, vref);
0219       for (auto& k : vref) {
0220         RecoChargedCandidateRef candref = RecoChargedCandidateRef(k);
0221         TrackRef tk = candref->get<TrackRef>();
0222         double pt = tk->pt();
0223         double err0 = tk->error(0);
0224         double abspar0 = fabs(tk->parameter(0));
0225         // convert to 90% efficiency threshold
0226         if (abspar0 > 0)
0227           pt += theNSigmas[i] * err0 / abspar0 * pt;
0228         if (pt > ptcut)
0229           nFound++;
0230       }
0231       if (nFound >= theNumberOfObjects) {
0232         hHLTeff[i]->Fill(ptcut, this_event_weight);
0233       } else {
0234         break;
0235       }
0236     }
0237   }
0238 }
0239 
0240 DEFINE_FWK_MODULE(HLTMuonRateAnalyzer);
[ 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!