Muon/test/HLTMuonTurnOnAnalyzer.cc

0001 /** \class HLTMuonTurnOnAnalyzer
0002  *  Get L1/HLT turn on curves
0003  *
0004  *  \author J. Alcaraz
0005  */
0006
0007 #include "HLTrigger/Muon/test/HLTMuonTurnOnAnalyzer.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/L1Trigger/interface/L1MuonParticle.h"
0017 #include "DataFormats/L1Trigger/interface/L1MuonParticleFwd.h"
0018
0019 #include "DataFormats/TrackReco/interface/Track.h"
0020 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
0021 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
0022
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024
0025 #include "TFile.h"
0026 #include "TH1F.h"
0027
0028 using namespace std;
0029 using namespace edm;
0030 using namespace reco;
0031 using namespace trigger;
0032 using namespace l1extra;
0033
0034 /// Constructor
0035 HLTMuonTurnOnAnalyzer::HLTMuonTurnOnAnalyzer(const ParameterSet& pset) {
0036   theGenLabel = pset.getUntrackedParameter<InputTag>("GenLabel");
0037   useMuonFromGenerator = pset.getUntrackedParameter<bool>("UseMuonFromGenerator");
0038   theL1CollectionLabel = pset.getUntrackedParameter<InputTag>("L1CollectionLabel");
0039   theHLTCollectionLabels = pset.getUntrackedParameter<std::vector<InputTag> >("HLTCollectionLabels");
0040   theGenToken = consumes<edm::HepMCProduct>(theGenLabel);
0041   theL1CollectionToken = consumes<trigger::TriggerFilterObjectWithRefs>(theL1CollectionLabel);
0042   for (auto& theHLTCollectionLabel : theHLTCollectionLabels) {
0043     theHLTCollectionTokens.push_back(consumes<trigger::TriggerFilterObjectWithRefs>(theHLTCollectionLabel));
0044   }
0045   theReferenceThreshold = pset.getUntrackedParameter<double>("ReferenceThreshold");
0046
0047   thePtMin = pset.getUntrackedParameter<double>("PtMin");
0048   thePtMax = pset.getUntrackedParameter<double>("PtMax");
0049   theNbins = pset.getUntrackedParameter<unsigned int>("Nbins");
0050
0051   theRootFileName = pset.getUntrackedParameter<string>("RootFileName");
0052
0053   theNumberOfEvents = 0.;
0054 }
0055
0056 /// Destructor
0057 HLTMuonTurnOnAnalyzer::~HLTMuonTurnOnAnalyzer() = default;
0058
0059 void HLTMuonTurnOnAnalyzer::beginJob() {
0060   // Create the root file
0061   theFile = new TFile(theRootFileName.c_str(), "RECREATE");
0062   theFile->cd();
0063
0064   char chname[256];
0065   char chtitle[256];
0066   snprintf(chname, 255, "eff_%s", theL1CollectionLabel.encode().c_str());
0067   snprintf(chtitle, 255, "Efficiency (%%) vs Generated Pt (GeV), label=%s", theL1CollectionLabel.encode().c_str());
0068   hL1eff = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0069   hL1nor = new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax);
0070
0071   for (auto& theHLTCollectionLabel : theHLTCollectionLabels) {
0072     snprintf(chname, 255, "eff_%s", theHLTCollectionLabel.encode().c_str());
0073     snprintf(chtitle, 255, "Efficiency (%%) vs Generated Pt (GeV), label=%s", theHLTCollectionLabel.encode().c_str());
0074     hHLTeff.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0075     hHLTnor.push_back(new TH1F(chname, chtitle, theNbins, thePtMin, thePtMax));
0076   }
0077 }
0078
0079 void HLTMuonTurnOnAnalyzer::endJob() {
0080   LogInfo("HLTMuonTurnOnAnalyzer") << " (Weighted) number of analyzed events= " << theNumberOfEvents;
0081   theFile->cd();
0082
0083   if (theNumberOfEvents == 0) {
0084     LogInfo("HLTMuonTurnOnAnalyzer") << " No histograms will be written because number of events=0!!!";
0085     theFile->Close();
0086     return;
0087   }
0088
0089   // L1 operations
0090   hL1eff->Divide(hL1nor);
0091   hL1eff->Scale(100.);
0092
0093   // HLT operations
0094   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0095     hHLTeff[i]->Divide(hHLTnor[i]);
0096     hHLTeff[i]->Scale(100.);
0097   }
0098
0099   // Write the histos to file
0100   hL1eff->Write();
0101   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0102     hHLTeff[i]->Write();
0103   }
0104
0105   theFile->Close();
0106 }
0107
0108 void HLTMuonTurnOnAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0109   theFile->cd();
0110
0111   // Get the HepMC product
0112   double this_event_weight = 1.;
0113   Handle<HepMCProduct> genProduct;
0114   event.getByToken(theGenToken, genProduct);
0115
0116   const HepMC::GenEvent* evt = genProduct->GetEvent();
0117   HepMC::WeightContainer weights = evt->weights();
0118   if (weights.size() > 0)
0119     this_event_weight = weights[0];
0120   theNumberOfEvents += this_event_weight;
0121
0122   // Get the L1 collection
0123   Handle<TriggerFilterObjectWithRefs> l1cands;
0124   event.getByToken(theL1CollectionToken, l1cands);
0125
0126   // Get the HLT collections
0127   std::vector<Handle<TriggerFilterObjectWithRefs> > hltcands(theHLTCollectionLabels.size());
0128
0129   unsigned int modules_in_this_event = 0;
0130   for (unsigned int i = 0; i < theHLTCollectionLabels.size(); i++) {
0131     event.getByToken(theHLTCollectionTokens[i], hltcands[i]);
0132     if (hltcands[i].failedToGet())
0133       break;
0134     modules_in_this_event++;
0135   }
0136
0137   // Get the muon with maximum pt at generator level or reconstruction, depending on the choice
0138   bool refmuon_found = false;
0139   double ptuse = -1;
0140
0141   if (useMuonFromGenerator) {
0142     HepMC::GenEvent::particle_const_iterator part;
0143     for (part = evt->particles_begin(); part != evt->particles_end(); ++part) {
0144       int id1 = (*part)->pdg_id();
0145       if (id1 != 13 && id1 != -13)
0146         continue;
0147       float pt1 = (*part)->momentum().perp();
0148       if (pt1 > ptuse) {
0149         refmuon_found = true;
0150         ptuse = pt1;
0151       }
0152     }
0153   } else {
0154     unsigned int i = modules_in_this_event - 1;
0155     vector<RecoChargedCandidateRef> vref;
0156     hltcands[i]->getObjects(TriggerMuon, vref);
0157     for (auto& k : vref) {
0158       RecoChargedCandidateRef candref = RecoChargedCandidateRef(k);
0159       TrackRef tk = candref->get<TrackRef>();
0160       double pt = tk->pt();
0161       if (pt > ptuse) {
0162         refmuon_found = true;
0163         ptuse = pt;
0164       }
0165     }
0166   }
0167
0168   if (!refmuon_found) {
0169     LogInfo("HLTMuonTurnOnAnalyzer") << " NO reference muon found!!!";
0170     LogInfo("HLTMuonTurnOnAnalyzer") << " Skipping event";
0171     return;
0172   }
0173
0174   // Fix L1 thresholds to obtain the efficiecy plot
0175   unsigned int nL1FoundRef = 0;
0176   double epsilon = 0.001;
0177   vector<L1MuonParticleRef> l1mu;
0178   l1cands->getObjects(TriggerL1Mu, l1mu);
0179   for (auto& k : l1mu) {
0180     L1MuonParticleRef candref = L1MuonParticleRef(k);
0181     // L1 PTs are "quantized" due to LUTs.
0182     // Their meaning: true_pt > ptLUT more than 90% pof the times
0183     double ptLUT = candref->pt();
0184     // Add "epsilon" to avoid rounding errors when ptLUT==L1Threshold
0185     if (ptLUT + epsilon > theReferenceThreshold)
0186       nL1FoundRef++;
0187   }
0188   hL1nor->Fill(ptuse, this_event_weight);
0189   if (nL1FoundRef > 0)
0190     hL1eff->Fill(ptuse, this_event_weight);
0191
0192   // HLT filling
0193   unsigned int last_module = modules_in_this_event - 1;
0194   if ((!useMuonFromGenerator) && last_module > 0)
0195     last_module--;
0196   for (unsigned int i = 0; i <= last_module; i++) {
0197     double ptcut = theReferenceThreshold;
0198     unsigned nFound = 0;
0199     vector<RecoChargedCandidateRef> vref;
0200     hltcands[i]->getObjects(TriggerMuon, vref);
0201     for (auto& k : vref) {
0202       RecoChargedCandidateRef candref = RecoChargedCandidateRef(k);
0203       TrackRef tk = candref->get<TrackRef>();
0204       double pt = tk->pt();
0205       if (pt > ptcut)
0206         nFound++;
0207     }
0208     hHLTnor[i]->Fill(ptuse, this_event_weight);
0209     if (nFound > 0)
0210       hHLTeff[i]->Fill(ptuse, this_event_weight);
0211   }
0212 }
0213
0214 DEFINE_FWK_MODULE(HLTMuonTurnOnAnalyzer);