L3MuonIsolationProducer/test/L3MuonIsolationAnalyzer.cc

0001 /** \class L3MuonIsolationAnalyzer
0002  *  Analyzer of HLT L3 muon isolation performance
0003  *
0004  *  \author J. Alcaraz
0005  */
0006
0007 #include "RecoMuon/L3MuonIsolationProducer/test/L3MuonIsolationAnalyzer.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 L3MuonIsolationAnalyzer::L3MuonIsolationAnalyzer(const ParameterSet& pset)
0031     : theIsolationLabel(pset.getUntrackedParameter<edm::InputTag>("IsolationCollectionLabel")),
0032       theConeCases(pset.getParameter<std::vector<double> >("ConeCases")),
0033       thePtMin(pset.getParameter<double>("PtMin")),
0034       thePtMax(pset.getParameter<double>("PtMax")),
0035       thePtBins(pset.getParameter<unsigned int>("PtBins")),
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|L3MuonIsolationTest") << " L3MuonIsolationTest constructor called";
0042 }
0043
0044 /// Destructor
0045 L3MuonIsolationAnalyzer::~L3MuonIsolationAnalyzer() {}
0046
0047 void L3MuonIsolationAnalyzer::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   hPtSum = new TH1F("ptSum", "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   hEffVsPt = new TH1F("effVsPt", "Efficiency(%) vs E_{T}^{weighted} (GeV)", thePtBins, thePtMin, thePtMax);
0060
0061   hEffVsPtArray.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), "effVsPt-%.2d-%.2d", j, k);
0067       snprintf(chtit,
0068                sizeof(chtit),
0069                "Eff(%%) vs P_{T}(GeV), cone size %.4f, eta in [%.4f,%.4f]",
0070                theConeCases[j],
0071                theCuts[k].etaRange.min(),
0072                theCuts[k].etaRange.max());
0073       hEffVsPtArray.push_back(new TH1F(chnam, chtit, thePtBins, thePtMin, thePtMax));
0074     }
0075   }
0076 }
0077
0078 void L3MuonIsolationAnalyzer::endJob() {
0079   Puts("L3MuonIsolationAnalyzer>>> FINAL PRINTOUTS -> BEGIN");
0080   Puts("L3MuonIsolationAnalyzer>>> Number of analyzed events= %d", numberOfEvents);
0081   Puts("L3MuonIsolationAnalyzer>>> Number of analyzed muons= %d", numberOfMuons);
0082
0083   // Write the histos to file
0084   theRootFile->cd();
0085
0086   hPtSum->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("");
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 = hEffVsPt->GetNbinsX() + 1;
0106   norm = hEffVsPt->GetBinContent(overflow_bin);
0107   if (norm <= 0.)
0108     norm = 1.;
0109   for (unsigned int k = 1; k < overflow_bin; k++) {
0110     float aux = hEffVsPt->GetBinContent(k);
0111     float eff = 100 * aux / norm;
0112     hEffVsPt->SetBinContent(k, eff);
0113     float pt = thePtMin + (k - 0.5) / thePtBins * (thePtMax - thePtMin);
0114     Puts("\tEfficiency for Pt threshold cut %f: %f", pt, eff);
0115   }
0116   hEffVsPt->Write();
0117
0118   for (unsigned int j = 0; j < hEffVsPtArray.size(); j++) {
0119     overflow_bin = hEffVsPtArray[j]->GetNbinsX() + 1;
0120     norm = hEffVsPtArray[j]->GetBinContent(overflow_bin);
0121     if (norm <= 0.)
0122       norm = 1.;
0123     fprintf(theTxtFile, "%s\n", hEffVsPtArray[j]->GetTitle());
0124     fprintf(theTxtFile, "%s\n", "PT threshold(GeV), efficiency(%): ");
0125     for (unsigned int k = 1; k < overflow_bin; k++) {
0126       float aux = hEffVsPtArray[j]->GetBinContent(k);
0127       float eff = 100 * aux / norm;
0128       hEffVsPtArray[j]->SetBinContent(k, eff);
0129       float ptthr = hEffVsPtArray[j]->GetXaxis()->GetBinCenter(k);
0130       fprintf(theTxtFile, "%6.2f %8.3f\n", ptthr, eff);
0131     }
0132     hEffVsPtArray[j]->Write();
0133   }
0134
0135   theRootFile->Close();
0136   fclose(theTxtFile);
0137
0138   Puts("L3MuonIsolationAnalyzer>>> FINAL PRINTOUTS -> END");
0139 }
0140
0141 void L3MuonIsolationAnalyzer::analyze(const Event& event, const EventSetup& eventSetup) {
0142   static const string metname = "L3MuonIsolation";
0143
0144   numberOfEvents++;
0145
0146   // Get the isolation information from the event
0147   Handle<reco::IsoDepositMap> depMap;
0148   event.getByLabel(theIsolationLabel, depMap);
0149
0150   //! How do I know it's made for tracks? Just guessing
0151   typedef edm::View<reco::Track> tracks_type;
0152   Handle<tracks_type> tracksH;
0153
0154   typedef reco::IsoDepositMap::const_iterator depmap_iterator;
0155   depmap_iterator depMI = depMap->begin();
0156   depmap_iterator depMEnd = depMap->end();
0157
0158   for (; depMI != depMEnd; ++depMI) {
0159     if (depMI.id() != tracksH.id()) {
0160       LogTrace(metname) << "Getting tracks with id " << depMI.id();
0161       event.get(depMI.id(), tracksH);
0162     }
0163
0164     typedef reco::IsoDepositMap::container::const_iterator dep_iterator;
0165     dep_iterator depI = depMI.begin();
0166     dep_iterator depEnd = depMI.end();
0167
0168     typedef tracks_type::const_iterator tk_iterator;
0169     tk_iterator tkI = tracksH->begin();
0170     tk_iterator tkEnd = tracksH->end();
0171
0172     for (; tkI != tkEnd && depI != depEnd; ++tkI, ++depI) {
0173       numberOfMuons++;
0174       tracks_type::const_pointer tk = &*tkI;
0175       const IsoDeposit& dep = *depI;
0176
0177       muonisolation::Cuts::CutSpec cuts_here = theCuts(tk->eta());
0178       double conesize = cuts_here.conesize;
0179       float ptsum = dep.depositWithin(conesize);
0180       hPtSum->Fill(ptsum);
0181
0182       hEffVsCone->Fill(theConeCases.size() + 999.);
0183       for (unsigned int j = 0; j < theConeCases.size(); j++) {
0184         if (dep.depositWithin(theConeCases[j]) < cuts_here.threshold)
0185           hEffVsCone->Fill(j + 1.0);
0186       }
0187
0188       hEffVsPt->Fill(thePtMax + 999.);
0189       for (unsigned int j = 0; j < thePtBins; j++) {
0190         float ptthr = thePtMin + (j + 0.5) / thePtBins * (thePtMax - thePtMin);
0191         if (ptsum < ptthr)
0192           hEffVsPt->Fill(ptthr);
0193       }
0194
0195       for (unsigned int j = 0; j < theConeCases.size(); j++) {
0196         float ptsum = dep.depositWithin(theConeCases[j]);
0197         for (unsigned int k = 0; k < theCuts.size(); k++) {
0198           unsigned int n = theCuts.size() * j + k;
0199           if (fabs(tk->eta()) < theCuts[k].etaRange.min())
0200             continue;
0201           if (fabs(tk->eta()) > theCuts[k].etaRange.max())
0202             continue;
0203           hEffVsPtArray[n]->Fill(thePtMax + 999.);
0204           for (unsigned int l = 0; l < thePtBins; l++) {
0205             float ptthr = thePtMin + (l + 0.5) / thePtBins * (thePtMax - thePtMin);
0206             if (ptsum < ptthr)
0207               hEffVsPtArray[n]->Fill(ptthr);
0208           }
0209         }
0210       }
0211     }
0212   }
0213   Puts("L3MuonIsolationAnalyzer>>> Run: %llu, Event %llu, number of muons %d",
0214        event.id().run(),
0215        event.id().event(),
0216        depMap->size());
0217 }
0218
0219 void L3MuonIsolationAnalyzer::Puts(const char* va_(fmt), ...) {
0220   // Do not write more than 256 characters
0221   const unsigned int bufsize = 256;
0222   char chout[bufsize] = "";
0223   va_list ap;
0224   va_start(ap, va_(fmt));
0225   vsnprintf(chout, bufsize, va_(fmt), ap);
0226   va_end(ap);
0227   LogVerbatim("") << chout;
0228 }
0229
0230 DEFINE_FWK_MODULE(L3MuonIsolationAnalyzer);