L2MuonIsolationProducer/test/L2MuonIsolationAnalyzer.cc

0001 /** \class L2MuonIsolationAnalyzer
0002  *  Analyzer of HLT L2 muon isolation performance
0003  *
0004  *  \author J. Alcaraz
0005  */
0006
0007 #include "RecoMuon/L2MuonIsolationProducer/test/L2MuonIsolationAnalyzer.h"
0008
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include "DataFormats/RecoCandidate/interface/IsoDeposit.h"
0017 #include "DataFormats/RecoCandidate/interface/IsoDepositFwd.h"
0018 #include "DataFormats/TrackReco/interface/Track.h"
0019
0020 #include "TFile.h"
0021 #include "TH1F.h"
0022
0023 #include "Varargs.h"
0024
0025 using namespace std;
0026 using namespace edm;
0027 using namespace reco;
0028
0029 /// Constructor
0030 L2MuonIsolationAnalyzer::L2MuonIsolationAnalyzer(const ParameterSet& pset)
0031     : theIsolationLabel(pset.getUntrackedParameter<edm::InputTag>("IsolationCollectionLabel")),
0032       theConeCases(pset.getParameter<std::vector<double> >("ConeCases")),
0033       theEtMin(pset.getParameter<double>("EtMin")),
0034       theEtMax(pset.getParameter<double>("EtMax")),
0035       theEtBins(pset.getParameter<unsigned int>("EtBins")),
0036       theCuts(pset.getParameter<std::vector<double> >("EtaBounds"),
0037               pset.getParameter<std::vector<double> >("ConeSizes"),
0038               pset.getParameter<std::vector<double> >("Thresholds")),
0039       theRootFileName(pset.getUntrackedParameter<string>("rootFileName")),
0040       theTxtFileName(pset.getUntrackedParameter<string>("txtFileName")) {
0041   LogDebug("Muon|RecoMuon|L2MuonIsolationTest") << " L2MuonIsolationTest constructor called";
0042 }
0043
0044 /// Destructor
0045 L2MuonIsolationAnalyzer::~L2MuonIsolationAnalyzer() {}
0046
0047 void L2MuonIsolationAnalyzer::beginJob() {
0048   // Create output files
0049   theTxtFile = fopen(theTxtFileName.data(), "w");
0050
0051   theRootFile = TFile::Open(theRootFileName.c_str(), "RECREATE");
0052   theRootFile->cd();
0053
0054   numberOfEvents = 0;
0055   numberOfMuons = 0;
0056
0057   hEtSum = new TH1F("etSum", "Sum E_{T}^{weighted} (GeV)", 100, 0, 50);
0058   hEffVsCone = new TH1F("effVsCone", "Efficiency(%) vs Cone Set", theConeCases.size(), 0.5, theConeCases.size() + 0.5);
0059   hEffVsEt = new TH1F("effVsEt", "Efficiency(%) vs E_{T}^{weighted} (GeV)", theEtBins, theEtMin, theEtMax);
0060
0061   hEffVsEtArray.clear();
0062   char chnam[256];
0063   char chtit[1000];
0064   for (unsigned int j = 0; j < theConeCases.size(); j++) {
0065     for (unsigned int k = 0; k < theCuts.size(); k++) {
0066       snprintf(chnam, sizeof(chnam), "effVsEt-%.2d-%.2d", j, k);
0067       snprintf(chtit,
0068                sizeof(chtit),
0069                "Eff(%%) vs E_{T}(GeV), cone size %.4f, eta in [%.4f,%.4f]",
0070                theConeCases[j],
0071                theCuts[k].etaRange.min(),
0072                theCuts[k].etaRange.max());
0073       hEffVsEtArray.push_back(new TH1F(chnam, chtit, theEtBins, theEtMin, theEtMax));
0074     }
0075   }
0076 }
0077
0078 void L2MuonIsolationAnalyzer::endJob() {
0079   Puts("L2MuonIsolationAnalyzer>>> FINAL PRINTOUTS -> BEGIN");
0080   Puts("L2MuonIsolationAnalyzer>>> Number of analyzed events= %d", numberOfEvents);
0081   Puts("L2MuonIsolationAnalyzer>>> Number of analyzed muons= %d", numberOfMuons);
0082
0083   // Write the histos to file
0084   theRootFile->cd();
0085
0086   hEtSum->Write();
0087
0088   unsigned int overflow_bin;
0089   float norm;
0090
0091   overflow_bin = hEffVsCone->GetNbinsX() + 1;
0092   norm = hEffVsCone->GetBinContent(overflow_bin);
0093   if (norm <= 0.)
0094     norm = 1.;
0095   Puts("EffVsCone Normalization= %f", norm);
0096   for (unsigned int k = 1; k < overflow_bin; k++) {
0097     float aux = hEffVsCone->GetBinContent(k);
0098     float eff = 100 * aux / norm;
0099     hEffVsCone->SetBinContent(k, eff);
0100     Puts("\tEfficiency in cone index= %d (size=%f): %f", k, theConeCases[k - 1], eff);
0101   }
0102   hEffVsCone->Write();
0103
0104   Puts("");
0105   overflow_bin = hEffVsEt->GetNbinsX() + 1;
0106   norm = hEffVsEt->GetBinContent(overflow_bin);
0107   if (norm <= 0.)
0108     norm = 1.;
0109   Puts("EffVsEt Normalization= %f", norm);
0110   for (unsigned int k = 1; k < overflow_bin; k++) {
0111     float aux = hEffVsEt->GetBinContent(k);
0112     float eff = 100 * aux / norm;
0113     hEffVsEt->SetBinContent(k, eff);
0114     float et = theEtMin + (k - 0.5) / theEtBins * (theEtMax - theEtMin);
0115     Puts("\tEfficiency for Et threshold cut %f: %f", et, eff);
0116   }
0117   hEffVsEt->Write();
0118
0119   for (unsigned int j = 0; j < hEffVsEtArray.size(); j++) {
0120     overflow_bin = hEffVsEtArray[j]->GetNbinsX() + 1;
0121     norm = hEffVsEtArray[j]->GetBinContent(overflow_bin);
0122     if (norm <= 0.)
0123       norm = 1.;
0124     Puts("EffVsEt[%d] Normalization= %f", j, norm);
0125     fprintf(theTxtFile, "%s\n", hEffVsEtArray[j]->GetTitle());
0126     fprintf(theTxtFile, "%s\n", "ET threshold(GeV), efficiency(%): ");
0127     for (unsigned int k = 1; k < overflow_bin; k++) {
0128       float aux = hEffVsEtArray[j]->GetBinContent(k);
0129       float eff = 100 * aux / norm;
0130       hEffVsEtArray[j]->SetBinContent(k, eff);
0131       float etthr = hEffVsEtArray[j]->GetXaxis()->GetBinCenter(k);
0132       fprintf(theTxtFile, "%6.2f %8.3f\n", etthr, eff);
0133     }
0134     hEffVsEtArray[j]->Write();
0135   }
0136
0137   theRootFile->Close();
0138   fclose(theTxtFile);
0139
0140   Puts("L2MuonIsolationAnalyzer>>> FINAL PRINTOUTS -> END");
0141 }
0142
0143 void L2MuonIsolationAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0144   static const string metname = "L3MuonIsolation";
0145
0146   numberOfEvents++;
0147
0148   // Get the isolation information from the event
0149   Handle<reco::IsoDepositMap> depMap;
0150   event.getByLabel(theIsolationLabel, depMap);
0151
0152   //! How do I know it's made for tracks? Just guessing
0153   typedef edm::View<reco::Track> tracks_type;
0154   Handle<tracks_type> tracksH;
0155
0156   typedef reco::IsoDepositMap::const_iterator depmap_iterator;
0157   depmap_iterator depMI = depMap->begin();
0158   depmap_iterator depMEnd = depMap->end();
0159
0160   for (; depMI != depMEnd; ++depMI) {
0161     if (depMI.id() != tracksH.id()) {
0162       LogTrace(metname) << "Getting tracks with id " << depMI.id();
0163       event.get(depMI.id(), tracksH);
0164     }
0165
0166     typedef reco::IsoDepositMap::container::const_iterator dep_iterator;
0167     dep_iterator depI = depMI.begin();
0168     dep_iterator depEnd = depMI.end();
0169
0170     typedef tracks_type::const_iterator tk_iterator;
0171     tk_iterator tkI = tracksH->begin();
0172     tk_iterator tkEnd = tracksH->end();
0173
0174     for (; tkI != tkEnd && depI != depEnd; ++tkI, ++depI) {
0175       numberOfMuons++;
0176       tracks_type::const_pointer tk = &*tkI;
0177       const IsoDeposit& dep = *depI;
0178       //Puts("L2MuonIsolationAnalyzer>>> Track pt: %f, Deposit eta%f", tk->pt(), dep.eta());
0179
0180       muonisolation::Cuts::CutSpec cuts_here = theCuts(tk->eta());
0181       double conesize = cuts_here.conesize;
0182       float etsum = dep.depositWithin(conesize);
0183       hEtSum->Fill(etsum);
0184
0185       hEffVsCone->Fill(theConeCases.size() + 999.);
0186       for (unsigned int j = 0; j < theConeCases.size(); j++) {
0187         if (dep.depositWithin(theConeCases[j]) < cuts_here.threshold)
0188           hEffVsCone->Fill(j + 1.0);
0189       }
0190
0191       hEffVsEt->Fill(theEtMax + 999.);
0192       for (unsigned int j = 0; j < theEtBins; j++) {
0193         float etthr = theEtMin + (j + 0.5) / theEtBins * (theEtMax - theEtMin);
0194         if (etsum < etthr)
0195           hEffVsEt->Fill(etthr);
0196       }
0197
0198       for (unsigned int j = 0; j < theConeCases.size(); j++) {
0199         float etsum = dep.depositWithin(theConeCases[j]);
0200         for (unsigned int k = 0; k < theCuts.size(); k++) {
0201           unsigned int n = theCuts.size() * j + k;
0202           if (fabs(tk->eta()) < theCuts[k].etaRange.min())
0203             continue;
0204           if (fabs(tk->eta()) > theCuts[k].etaRange.max())
0205             continue;
0206           hEffVsEtArray[n]->Fill(theEtMax + 999.);
0207           for (unsigned int l = 0; l < theEtBins; l++) {
0208             float etthr = theEtMin + (l + 0.5) / theEtBins * (theEtMax - theEtMin);
0209             if (etsum < etthr)
0210               hEffVsEtArray[n]->Fill(etthr);
0211           }
0212         }
0213       }
0214     }
0215   }
0216   Puts("L2MuonIsolationAnalyzer>>> Run: %llu, Event %llu, number of muons %d",
0217        event.id().run(),
0218        event.id().event(),
0219        depMap->size());
0220 }
0221
0222 void L2MuonIsolationAnalyzer::Puts(const char* va_(fmt), ...) {
0223   // Do not write more than 256 characters
0224   const unsigned int bufsize = 256;
0225   char chout[bufsize] = "";
0226   va_list ap;
0227   va_start(ap, va_(fmt));
0228   vsnprintf(chout, bufsize, va_(fmt), ap);
0229   va_end(ap);
0230   LogVerbatim("") << chout;
0231   //LogError("") << chout;
0232 }
0233
0234 DEFINE_FWK_MODULE(L2MuonIsolationAnalyzer);