Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:33:27 

0001 // -*- C++ -*-
0002 // MuIsoValidation.cc
0003 // Package:    Muon Isolation Validation
0004 // Class:      MuIsoValidation
0005 //
0006 /*
0007 
0008 
0009  Description: Muon Isolation Validation class
0010 
0011  Implementation: This code will accept a data set and generate plots of
0012     various quantities relevent to the Muon Isolation module. We will 
0013     be using the IsoDeposit class, *not* the MuonIsolation struct.
0014      
0015     The sequence of events is... 
0016         * initalize statics (which variables to plot, axis limtis, etc.)
0017         * run over events
0018             * for each event, run over the muons in that event
0019                 *record IsoDeposit data
0020         * transfer data to histograms, profile plots
0021         * save histograms to a root file
0022         
0023     Easy-peasy.
0024     
0025 */
0026 //
0027 // Original Author:  "C. Jess Riedel", UC Santa Barbara
0028 //         Created:  Tue Jul 17 15:58:24 CDT 2007
0029 //
0030 
0031 //Class header file
0032 #include "Validation/MuonIsolation/interface/MuIsoValidation.h"
0033 
0034 //System included files
0035 #include <memory>
0036 #include <string>
0037 #include <typeinfo>
0038 #include <utility>
0039 
0040 //Root included files
0041 #include "TH1.h"
0042 #include "TH2.h"
0043 #include "TProfile.h"
0044 
0045 //Event framework included files
0046 #include "FWCore/Framework/interface/MakerMacros.h"
0047 #include "FWCore/Framework/interface/Event.h"
0048 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0049 
0050 //Other included files
0051 #include "DataFormats/TrackReco/interface/Track.h"
0052 
0053 //Using declarations
0054 using std::pair;
0055 using std::string;
0056 using std::vector;
0057 
0058 //
0059 //-----------------Constructors---------------------
0060 //
0061 MuIsoValidation::MuIsoValidation(const edm::ParameterSet& ps) {
0062   iConfig = ps;
0063 
0064   //  rootfilename = iConfig.getUntrackedParameter<string>("rootfilename"); // comment out for inclusion
0065   requireCombinedMuon = iConfig.getUntrackedParameter<bool>("requireCombinedMuon");
0066   dirName = iConfig.getUntrackedParameter<std::string>("directory");
0067   //subDirName = iConfig.getParameter<std::string>("@module_label");
0068 
0069   //dirName += subDirName;
0070 
0071   //--------Initialize tags-------
0072   Muon_Tag = iConfig.getUntrackedParameter<edm::InputTag>("Global_Muon_Label");
0073   Muon_Token = consumes<edm::View<reco::Muon> >(Muon_Tag);
0074 
0075   //-------Initialize counters----------------
0076   nEvents = 0;
0077   nIncMuons = 0;
0078   //  nCombinedMuons = 0;
0079 
0080   InitStatics();
0081 
0082   //Set up DAQ
0083   subsystemname_ = iConfig.getUntrackedParameter<std::string>("subSystemFolder", "YourSubsystem");
0084 
0085   //------"allocate" space for the data vectors-------
0086   h_1D.resize(NUM_VARS);
0087   cd_plots.resize(NUM_VARS);
0088   p_2D.resize(NUM_VARS, vector<MonitorElement*>(NUM_VARS));
0089 }
0090 
0091 //
0092 //----------Destructor-----------------
0093 //
0094 MuIsoValidation::~MuIsoValidation() {
0095   //Deallocate memory
0096 }
0097 
0098 //
0099 //------------Methods-----------------------
0100 //
0101 void MuIsoValidation::InitStatics() {
0102   //-----------Initialize primatives-----------
0103   S_BIN_WIDTH = 1.0;  //in GeV
0104   L_BIN_WIDTH = 2.0;  //in GeV
0105   LOG_BINNING_ENABLED = 1;
0106   NUM_LOG_BINS = 15;
0107   LOG_BINNING_RATIO = 1.1;  //ratio by which each bin is wider than the last for log binning
0108   //i.e.  bin widths are (x), (r*x), (r^2*x), ..., (r^(nbins)*x)
0109 
0110   //-------Initialize Titles---------
0111   title_sam = "";   //"[Sample b-jet events] ";
0112   title_cone = "";  //" [in R=0.3 IsoDeposit Cone]";
0113   //The above two pieces of info will be printed on the title of the whole page,
0114   //not for each individual histogram
0115   title_cd = "C.D. of ";
0116 
0117   //-------"Allocate" memory for vectors
0118   main_titles.resize(NUM_VARS);
0119   axis_titles.resize(NUM_VARS);
0120   names.resize(NUM_VARS);
0121   param.resize(NUM_VARS, vector<double>(3));
0122   isContinuous.resize(NUM_VARS);
0123   cdCompNeeded.resize(NUM_VARS);
0124 
0125   //-----Titles of the plots-----------
0126   main_titles[0] = "Total Tracker Momentum";
0127   main_titles[1] = "Total EM Cal Energy";
0128   main_titles[2] = "Total Had Cal Energy";
0129   main_titles[3] = "Total HO Cal Energy";
0130   main_titles[4] = "Number of Tracker Tracks";
0131   main_titles[5] = "Number of Jets around Muon";
0132   main_titles[6] = "Tracker p_{T} within veto cone";
0133   main_titles[7] = "EM E_{T} within veto cone";
0134   main_titles[8] = "Had E_{T} within veto cone";
0135   main_titles[9] = "HO E_{T} within veto cone";
0136   main_titles[10] = "Muon p_{T}";
0137   main_titles[11] = "Muon #eta";
0138   main_titles[12] = "Muon #phi";
0139   main_titles[13] = "Average Momentum per Track ";
0140   main_titles[14] = "Weighted Energy";
0141   main_titles[15] = "PF Sum of Charged Hadron Pt";
0142   main_titles[16] = "PF Sum of Total Hadron Pt";
0143   main_titles[17] = "PF Sum of E,Mu Pt";
0144   main_titles[18] = "PF Sum of Neutral Hadron Et";
0145   main_titles[19] = "PF Sum of Photon Et";
0146   main_titles[20] = "PF Sum of Pt from non-PV";
0147 
0148   //------Titles on the X or Y axis------------
0149   axis_titles[0] = "#Sigma p_{T}   (GeV)";
0150   axis_titles[1] = "#Sigma E_{T}^{EM}   (GeV)";
0151   axis_titles[2] = "#Sigma E_{T}^{Had}   (GeV)";
0152   axis_titles[3] = "#Sigma E_{T}^{HO}   (GeV)";
0153   axis_titles[4] = "N_{Tracks}";
0154   axis_titles[5] = "N_{Jets}";
0155   axis_titles[6] = "#Sigma p_{T,veto} (GeV)";
0156   axis_titles[7] = "#Sigma E_{T,veto}^{EM}   (GeV)";
0157   axis_titles[8] = "#Sigma E_{T,veto}^{Had}   (GeV)";
0158   axis_titles[9] = "#Sigma E_{T,veto}^{HO}   (GeV)";
0159   axis_titles[10] = "p_{T,#mu} (GeV)";
0160   axis_titles[11] = "#eta_{#mu}";
0161   axis_titles[12] = "#phi_{#mu}";
0162   axis_titles[13] = "#Sigma p_{T} / N_{Tracks} (GeV)";
0163   axis_titles[14] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
0164   axis_titles[15] = "#Sigma p_{T}^{PFHadCha}   (GeV)";
0165   axis_titles[16] = "#Sigma p_{T}^{PFTotCha}   (GeV)";
0166   axis_titles[17] = "#Sigma p_{T}^{PFEMu}   (GeV)";
0167   axis_titles[18] = "#Sigma E_{T}^{PFHadNeu}   (GeV)";
0168   axis_titles[19] = "#Sigma E_{T}^{PFPhot}   (GeV)";
0169   axis_titles[20] = "#Sigma p_{T}^{PFPU}   (GeV)";
0170 
0171   //-----------Names given for the root file----------
0172   names[0] = "sumPt";
0173   names[1] = "emEt";
0174   names[2] = "hadEt";
0175   names[3] = "hoEt";
0176   names[4] = "nTracks";
0177   names[5] = "nJets";
0178   names[6] = "trackerVetoPt";
0179   names[7] = "emVetoEt";
0180   names[8] = "hadVetoEt";
0181   names[9] = "hoVetoEt";
0182   names[10] = "muonPt";
0183   names[11] = "muonEta";
0184   names[12] = "muonPhi";
0185   names[13] = "avgPt";
0186   names[14] = "weightedEt";
0187   names[15] = "PFsumChargedHadronPt";
0188   names[16] = "PFsumChargedTotalPt";
0189   names[17] = "PFsumEMuPt";
0190   names[18] = "PFsumNeutralHadronEt";
0191   names[19] = "PFsumPhotonEt";
0192   names[20] = "PFsumPUPt";
0193 
0194   //----------Parameters for binning of histograms---------
0195   //param[var][0] is the number of bins
0196   //param[var][1] is the low edge of the low bin
0197   //param[var][2] is the high edge of the high bin
0198   //
0199   // maximum value------,
0200   //                    |
0201   //                    V
0202   param[0][0] = (int)(20.0 / S_BIN_WIDTH);
0203   param[0][1] = 0.0;
0204   param[0][2] = param[0][0] * S_BIN_WIDTH;
0205   param[1][0] = (int)(20.0 / S_BIN_WIDTH);
0206   param[1][1] = 0.0;
0207   param[1][2] = param[1][0] * S_BIN_WIDTH;
0208   param[2][0] = (int)(20.0 / S_BIN_WIDTH);
0209   param[2][1] = 0.0;
0210   param[2][2] = param[2][0] * S_BIN_WIDTH;
0211   param[3][0] = 20;
0212   param[3][1] = 0.0;
0213   param[3][2] = 2.0;
0214   param[4][0] = 16;
0215   param[4][1] = -0.5;
0216   param[4][2] = param[4][0] - 0.5;
0217   param[5][0] = 4;
0218   param[5][1] = -0.5;
0219   param[5][2] = param[5][0] - 0.5;
0220   param[6][0] = (int)(40.0 / S_BIN_WIDTH);
0221   param[6][1] = 0.0;
0222   param[6][2] = param[6][0] * S_BIN_WIDTH;
0223   param[7][0] = 20;
0224   param[7][1] = 0.0;
0225   param[7][2] = 10.0;
0226   param[8][0] = (int)(20.0 / S_BIN_WIDTH);
0227   param[8][1] = 0.0;
0228   param[8][2] = param[8][0] * S_BIN_WIDTH;
0229   param[9][0] = 20;
0230   param[9][1] = 0.0;
0231   param[9][2] = 5.0;
0232   param[10][0] = (int)(40.0 / S_BIN_WIDTH);
0233   param[10][1] = 0.0;
0234   param[10][2] = param[10][0] * S_BIN_WIDTH;
0235   param[11][0] = 24;
0236   param[11][1] = -2.4;
0237   param[11][2] = 2.4;
0238   param[12][0] = 32;
0239   param[12][1] = -3.2;
0240   param[12][2] = 3.2;
0241   param[13][0] = (int)(15.0 / S_BIN_WIDTH);
0242   param[13][1] = 0.0;
0243   param[13][2] = param[13][0] * S_BIN_WIDTH;
0244   param[14][0] = (int)(20.0 / S_BIN_WIDTH);
0245   param[14][1] = 0.0;
0246   param[14][2] = param[14][0] * S_BIN_WIDTH;
0247   param[15][0] = (int)(20.0 / S_BIN_WIDTH);
0248   param[15][1] = 0.0;
0249   param[15][2] = param[15][0] * S_BIN_WIDTH;
0250   param[16][0] = (int)(20.0 / S_BIN_WIDTH);
0251   param[15][1] = 0.0;
0252   param[16][2] = param[16][0] * S_BIN_WIDTH;
0253   param[17][0] = (int)(20.0 / S_BIN_WIDTH) + 1;
0254   param[17][1] = -S_BIN_WIDTH;
0255   param[17][2] = param[17][0] * S_BIN_WIDTH;
0256   param[18][0] = (int)(20.0 / S_BIN_WIDTH);
0257   param[18][1] = 0.0;
0258   param[18][2] = param[18][0] * S_BIN_WIDTH;
0259   param[19][0] = (int)(20.0 / S_BIN_WIDTH);
0260   param[19][1] = 0.0;
0261   param[19][2] = param[19][0] * S_BIN_WIDTH;
0262   param[20][0] = (int)(20.0 / S_BIN_WIDTH);
0263   param[20][1] = 0.0;
0264   param[20][2] = param[20][0] * S_BIN_WIDTH;
0265 
0266   //--------------Is the variable continuous (i.e. non-integer)?-------------
0267   //---------(Log binning will only be used for continuous variables)--------
0268   isContinuous[0] = 1;
0269   isContinuous[1] = 1;
0270   isContinuous[2] = 1;
0271   isContinuous[3] = 1;
0272   isContinuous[4] = 0;
0273   isContinuous[5] = 0;
0274   isContinuous[6] = 1;
0275   isContinuous[7] = 1;
0276   isContinuous[8] = 1;
0277   isContinuous[9] = 1;
0278   isContinuous[10] = 1;
0279   isContinuous[11] = 1;
0280   isContinuous[12] = 1;
0281   isContinuous[13] = 1;
0282   isContinuous[14] = 1;
0283   isContinuous[15] = 1;
0284   isContinuous[16] = 1;
0285   isContinuous[17] = 1;
0286   isContinuous[18] = 1;
0287   isContinuous[19] = 1;
0288   isContinuous[20] = 1;
0289 
0290   //----Should the cumulative distribution be calculated for this variable?-----
0291   cdCompNeeded[0] = 1;
0292   cdCompNeeded[1] = 1;
0293   cdCompNeeded[2] = 1;
0294   cdCompNeeded[3] = 1;
0295   cdCompNeeded[4] = 1;
0296   cdCompNeeded[5] = 1;
0297   cdCompNeeded[6] = 1;
0298   cdCompNeeded[7] = 1;
0299   cdCompNeeded[8] = 1;
0300   cdCompNeeded[9] = 1;
0301   cdCompNeeded[10] = 0;
0302   cdCompNeeded[11] = 0;
0303   cdCompNeeded[12] = 0;
0304   cdCompNeeded[13] = 1;
0305   cdCompNeeded[14] = 1;
0306   cdCompNeeded[15] = 1;
0307   cdCompNeeded[16] = 1;
0308   cdCompNeeded[17] = 1;
0309   cdCompNeeded[18] = 1;
0310   cdCompNeeded[19] = 1;
0311   cdCompNeeded[20] = 1;
0312 }
0313 
0314 // ------------ method called for each event  ------------
0315 void MuIsoValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0316   ++nEvents;
0317   edm::LogInfo("Tutorial") << "\nInvestigating event #" << nEvents << "\n";
0318 
0319   // Get Muon Collection
0320   edm::Handle<edm::View<reco::Muon> > muonsHandle;  //
0321   iEvent.getByToken(Muon_Token, muonsHandle);
0322 
0323   //Fill event entry in histogram of number of muons
0324   edm::LogInfo("Tutorial") << "Number of Muons: " << muonsHandle->size();
0325   theMuonData = muonsHandle->size();
0326   h_nMuons->Fill(theMuonData);
0327 
0328   //Fill historgams concerning muon isolation
0329   uint iMuon = 0;
0330   for (MuonIterator muon = muonsHandle->begin(); muon != muonsHandle->end(); ++muon, ++iMuon) {
0331     ++nIncMuons;
0332     if (requireCombinedMuon) {
0333       if (muon->combinedMuon().isNull())
0334         continue;
0335     }
0336     RecordData(muon);
0337     FillHistos();
0338   }
0339 }
0340 
0341 //---------------Record data for a signle muon's data---------------------
0342 void MuIsoValidation::RecordData(MuonIterator muon) {
0343   theData[0] = muon->isolationR03().sumPt;
0344   theData[1] = muon->isolationR03().emEt;
0345   theData[2] = muon->isolationR03().hadEt;
0346   theData[3] = muon->isolationR03().hoEt;
0347 
0348   theData[4] = muon->isolationR03().nTracks;
0349   theData[5] = muon->isolationR03().nJets;
0350   theData[6] = muon->isolationR03().trackerVetoPt;
0351   theData[7] = muon->isolationR03().emVetoEt;
0352   theData[8] = muon->isolationR03().hadVetoEt;
0353   theData[9] = muon->isolationR03().hoVetoEt;
0354 
0355   theData[10] = muon->pt();
0356   theData[11] = muon->eta();
0357   theData[12] = muon->phi();
0358 
0359   // make sure nTracks != 0 before filling this one
0360   if (theData[4] != 0)
0361     theData[13] = (double)theData[0] / (double)theData[4];
0362   else
0363     theData[13] = -99;
0364 
0365   theData[14] = 1.5 * theData[1] + theData[2];
0366 
0367   // Now PF isolation
0368   theData[15] = -99.;
0369   theData[16] = -99.;
0370   theData[17] = -99.;
0371   theData[18] = -99.;
0372   theData[19] = -99.;
0373   theData[20] = -99.;
0374   if (muon->isPFMuon() && muon->isPFIsolationValid()) {
0375     theData[15] = muon->pfIsolationR03().sumChargedHadronPt;
0376     theData[16] = muon->pfIsolationR03().sumChargedParticlePt;
0377     theData[17] = muon->pfIsolationR03().sumChargedParticlePt - muon->pfIsolationR03().sumChargedHadronPt;
0378     theData[18] = muon->pfIsolationR03().sumNeutralHadronEt;
0379     theData[19] = muon->pfIsolationR03().sumPhotonEt;
0380     theData[20] = muon->pfIsolationR03().sumPUPt;
0381   }
0382 }
0383 
0384 void MuIsoValidation::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
0385   ibooker.setCurrentFolder(dirName);
0386   //---initialize number of muons histogram---
0387   h_nMuons = ibooker.book1D("nMuons", title_sam + "Number of Muons", 20, 0., 20.);
0388   h_nMuons->setAxisTitle("Number of Muons", XAXIS);
0389   h_nMuons->setAxisTitle("Fraction of Events", YAXIS);
0390 
0391   //---Initialize 1D Histograms---
0392   for (int var = 0; var < NUM_VARS; var++) {
0393     h_1D[var] = ibooker.book1D(names[var],
0394                                title_sam + main_titles[var] + title_cone,
0395                                (int)param[var][0],
0396                                param[var][1],
0397                                param[var][2],
0398                                [](TH1* th1) { th1->Sumw2(); });
0399     h_1D[var]->setAxisTitle(axis_titles[var], XAXIS);
0400     h_1D[var]->setAxisTitle("Fraction of Muons", YAXIS);
0401 
0402     if (cdCompNeeded[var]) {
0403       cd_plots[var] = ibooker.book1D(names[var] + "_cd",
0404                                      title_sam + title_cd + main_titles[var] + title_cone,
0405                                      (int)param[var][0],
0406                                      param[var][1],
0407                                      param[var][2],
0408                                      [](TH1* th1) { th1->Sumw2(); });
0409       cd_plots[var]->setAxisTitle(axis_titles[var], XAXIS);
0410       cd_plots[var]->setAxisTitle("Fraction of Muons", YAXIS);
0411     }
0412   }  //Finish 1D
0413 
0414   //---Initialize 2D Histograms---
0415   for (int var1 = 0; var1 < NUM_VARS; var1++) {
0416     for (int var2 = 0; var2 < NUM_VARS; var2++) {
0417       if (var1 == var2)
0418         continue;
0419 
0420       //Monitor elements is weird and takes y axis parameters as well
0421       //as x axis parameters for a 1D profile plot
0422       p_2D[var1][var2] = ibooker.bookProfile(names[var1] + "_" + names[var2],
0423                                              title_sam + main_titles[var2] + " <vs> " + main_titles[var1] + title_cone,
0424                                              (int)param[var1][0],
0425                                              param[var1][1],
0426                                              param[var1][2],
0427                                              (int)param[var2][0],  //documentation says this is disregarded
0428                                              param[var2][1],       //does this do anything?
0429                                              param[var2][2],       //does this do anything?
0430                                              " ",                  //profile errors = spread/sqrt(num_datums)
0431                                              [&](TProfile* tprof) {
0432                                                if (LOG_BINNING_ENABLED && isContinuous[var1]) {
0433                                                  Double_t* bin_edges = new Double_t[NUM_LOG_BINS + 1];
0434                                                  // nbins+1 because there is one more edge than there are bins
0435                                                  MakeLogBinsForProfile(bin_edges, param[var1][1], param[var1][2]);
0436                                                  tprof->SetBins(NUM_LOG_BINS, bin_edges);
0437                                                  delete[] bin_edges;
0438                                                }
0439                                              });
0440 
0441       p_2D[var1][var2]->setAxisTitle(axis_titles[var1], XAXIS);
0442       p_2D[var1][var2]->setAxisTitle(axis_titles[var2], YAXIS);
0443     }
0444   }  //Finish 2D
0445 
0446   //avg pT not defined for zero tracks.
0447   //MonitorElement is inflxible and won't let me change the
0448   //number of bins!  I guess all I'm doing here is changing
0449   //range of the x axis when it is printed, not the actual
0450   //bins that are filled
0451   p_2D[4][9]->setAxisRange(0.5, 15.5, XAXIS);
0452 }
0453 
0454 void MuIsoValidation::MakeLogBinsForProfile(Double_t* bin_edges, const double min, const double max) {
0455   const double& r = LOG_BINNING_RATIO;
0456   const int& nbins = NUM_LOG_BINS;
0457 
0458   const double first_bin_width = (r > 1.0) ?  //so we don't divide by zero
0459                                      (max - min) * (1 - r) / (1 - pow(r, nbins))
0460                                            : (max - min) / nbins;
0461 
0462   bin_edges[0] = min;
0463   bin_edges[1] = min + first_bin_width;
0464   for (int n = 2; n < nbins; ++n) {
0465     bin_edges[n] = bin_edges[n - 1] + (bin_edges[n - 1] - bin_edges[n - 2]) * r;
0466   }
0467   bin_edges[nbins] = max;
0468 }
0469 
0470 void MuIsoValidation::FillHistos() {
0471   //----------Fill 1D histograms---------------
0472   for (int var = 0; var < NUM_VARS; var++) {
0473     h_1D[var]->Fill(theData[var]);
0474     if (cdCompNeeded[var])
0475       cd_plots[var]->Fill(theData[var]);  //right now, this is a regular PDF (just like h_1D)
0476   }  //Finish 1D
0477 
0478   //----------Fill 2D histograms---------------
0479   for (int var1 = 0; var1 < NUM_VARS; ++var1) {
0480     for (int var2 = 0; var2 < NUM_VARS; ++var2) {
0481       if (var1 == var2)
0482         continue;
0483       //change below to regular int interating!
0484       //      h_2D[var1][var2]->Fill(theData[var1], theData[var2]);
0485       p_2D[var1][var2]->Fill(theData[var1], theData[var2]);
0486     }
0487   }  //Finish 2D
0488 }
0489 
0490 //define this as a plug-in
0491 DEFINE_FWK_MODULE(MuIsoValidation);
[ 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!