Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:09:40 

0001 // -*- C++ -*-
0002 // MuonIsolationDQM.cc
0003 // Package:    Muon Isolation DQM
0004 // Class:      MuonIsolationDQM
0005 //
0006 /*
0007 
0008  Description: Muon Isolation DQM class
0009 
0010  Implementation: This code will accept a data set and generate plots of
0011     various quantities relevent to the Muon Isolation module. We will 
0012     be using the IsoDeposit class, *not* the MuonIsolation struct.
0013      
0014     The sequence of events is... 
0015         * initalize statics (which variables to plot, axis limtis, etc.)
0016         * run over events
0017             * for each event, run over the muons in that event
0018                 *record IsoDeposit data
0019         * transfer data to histograms, profile plots
0020         * save histograms to a root file
0021 */
0022 //#define DEBUG
0023 
0024 //Class header file
0025 #include "DQMOffline/Muon/interface/MuonIsolationDQM.h"
0026 
0027 //System included files
0028 #include <memory>
0029 #include <string>
0030 #include <typeinfo>
0031 #include <utility>
0032 
0033 //Root included files
0034 #include "TH1.h"
0035 #include "TH2.h"
0036 #include "TProfile.h"
0037 
0038 //Event framework included files
0039 #include "FWCore/Framework/interface/MakerMacros.h"
0040 #include "FWCore/Framework/interface/Event.h"
0041 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0042 
0043 //Other included files
0044 
0045 //Using declarations
0046 using std::pair;
0047 using std::string;
0048 using std::vector;
0049 
0050 using namespace std;
0051 //
0052 //-----------------Constructors---------------------
0053 //
0054 
0055 MuonIsolationDQM::MuonIsolationDQM(const edm::ParameterSet& iConfig) {
0056 #ifdef DEBUG
0057   cout << " Initialise Constructor " << endl;
0058 #endif
0059   requireSTAMuon = iConfig.getUntrackedParameter<bool>("requireSTAMuon");
0060   requireTRKMuon = iConfig.getUntrackedParameter<bool>("requireTRKMuon");
0061   requireGLBMuon = iConfig.getUntrackedParameter<bool>("requireGLBMuon");
0062   dirName = iConfig.getParameter<std::string>("directory");
0063 
0064   vtxBin_ = iConfig.getParameter<int>("vtxBin");
0065   vtxMin_ = iConfig.getParameter<double>("vtxMin");
0066   vtxMax_ = iConfig.getParameter<double>("vtxMax");
0067 
0068   //--------Initialize tags-------
0069   theMuonCollectionLabel_ =
0070       consumes<edm::View<reco::Muon> >(iConfig.getUntrackedParameter<edm::InputTag>("Global_Muon_Label"));
0071   theVertexCollectionLabel_ =
0072       consumes<reco::VertexCollection>(iConfig.getUntrackedParameter<edm::InputTag>("vertexLabel"));
0073 
0074   //-------Initialize Counterse----------------
0075   nEvents = 0;
0076   nSTAMuons = 0;
0077   nTRKMuons = 0;
0078   nGLBMuons = 0;
0079 
0080   InitStatics();
0081 
0082   //------"allocate" space for the data vectors-------
0083   h_1D.resize(NUM_VARS);
0084   h_2D.resize(NUM_VARS_2D);
0085   h_1D_NVTX.resize(NUM_VARS_NVTX);
0086 }
0087 
0088 //
0089 //----------Destructor-----------------
0090 //
0091 MuonIsolationDQM::~MuonIsolationDQM() {
0092 #ifdef DEBUG
0093   cout << "Calling destructor" << endl;
0094 #endif
0095   //Deallocate memory
0096 }
0097 
0098 //
0099 //------------Methods-----------------------
0100 //
0101 void MuonIsolationDQM::InitStatics() {
0102 #ifdef DEBUG
0103   cout << " InitStatistics() " << endl;
0104 #endif
0105   //-----------Initialize primitives-----------
0106   S_BIN_WIDTH = 1.0;  //in GeV
0107   L_BIN_WIDTH = 2.0;  //in GeV
0108   LOG_BINNING_ENABLED = 1;
0109   NUM_LOG_BINS = 15;
0110   LOG_BINNING_RATIO = 1.1;
0111   //ratio by which each bin is wider than the last for log binning
0112   //i.e.  bin widths are (x), (r*x), (r^2*x), ..., (r^(nbins)*x)
0113 
0114   //-------Initialize Titles---------
0115   title_sam = "";   //"[Sample b-jet events] ";
0116   title_cone = "";  //" [in R=0.3 IsoDeposit Cone]";
0117   //The above two pieces of info will be printed on the title of the whole page,
0118   //not for each individual histogram
0119   //  title_cd = "C.D. of ";
0120 
0121   //-------"Allocate" memory for vectors
0122   main_titles.resize(NUM_VARS);
0123   axis_titles.resize(NUM_VARS);
0124   names.resize(NUM_VARS);
0125   param.resize(NUM_VARS, vector<double>(3));
0126   isContinuous.resize(NUM_VARS);
0127 
0128   titles_2D.resize(NUM_VARS_2D);
0129   names_2D.resize(NUM_VARS_2D);
0130 
0131   main_titles_NVtxs.resize(NUM_VARS_NVTX);
0132   axis_titles_NVtxs.resize(NUM_VARS_NVTX);
0133   names_NVtxs.resize(NUM_VARS_NVTX);
0134 
0135 #ifdef DEBUG
0136   cout << "InitStatistics(): vectors resized " << endl;
0137 #endif
0138   //-----Titles of the plots-----------
0139   main_titles[0] = "Total Tracker Momentum, #Delta R = 0.3";
0140   main_titles[1] = "Total EM Cal Energy, #Delta R = 0.3";
0141   main_titles[2] = "Total Had Cal Energy, #Delta R = 0.3";
0142   main_titles[3] = "Total HO Cal Energy, #Delta R = 0.3";
0143   main_titles[4] = "Number of Tracker Tracks, #Delta R = 0.3";
0144   main_titles[5] = "Number of Jets around Muon, #Delta R = 0.3";
0145   main_titles[6] = "Tracker p_{T} within veto cone, #Delta R = 0.3";
0146   main_titles[7] = "EM E_{T} within veto cone, #Delta R = 0.3";
0147   main_titles[8] = "Had E_{T} within veto cone, #Delta R = 0.3";
0148   main_titles[9] = "HO E_{T} within veto cone, #Delta R = 0.3";
0149   main_titles[10] = "Average Momentum per Track, #Delta R = 0.3";
0150   main_titles[11] = "Weighted Energy, #Delta R = 0.3";
0151 
0152   main_titles[12] = "Total Tracker Momentum, #Delta R = 0.5";
0153   main_titles[13] = "Total EM Cal Energy, #Delta R = 0.5";
0154   main_titles[14] = "Total Had Cal Energy, #Delta R = 0.5";
0155   main_titles[15] = "Total HO Cal Energy, #Delta R = 0.5";
0156   main_titles[16] = "Number of Tracker Tracks, #Delta R = 0.5";
0157   main_titles[17] = "Number of Jets around Muon, #Delta R = 0.5";
0158   main_titles[18] = "Tracker p_{T} within veto cone, #Delta R = 0.5";
0159   main_titles[19] = "EM E_{T} within veto cone, #Delta R = 0.5";
0160   main_titles[20] = "Had E_{T} within veto cone, #Delta R = 0.5";
0161   main_titles[21] = "HO E_{T} within veto cone, #Delta R = 0.5";
0162   main_titles[22] = "Average Momentum per Track, #Delta R = 0.5";
0163   main_titles[23] = "Weighted Energy, #Delta R = 0.5";
0164 
0165   main_titles[24] = "Relative Detector-Based Isolation, #Delta R = 0.3";
0166   main_titles[25] = "Relative Detector-Based Isolation, #Delta R = 0.5";
0167 
0168   //-----Titles of the plots-----------
0169   main_titles[26] = "Sum PF Charged Hadron Pt, #Delta R = 0.3";
0170   main_titles[27] = "Sum PF Neutral Hadron Pt, #Delta R = 0.3";
0171   main_titles[28] = "Sum PF Photon Et, #Delta R = 0.3";
0172   main_titles[29] = "Sum PF Neutral Hadron Pt (Higher Pt threshold), #Delta R = 0.3";
0173   main_titles[30] = "Sum PF Photon Et (Higher Pt threshold), #Delta R = 0.3";
0174   main_titles[31] = "Sum PF Charged Particles Pt not from PV  (for Pu corrections), #Delta R = 0.3";
0175 
0176   //-----Titles of the plots-----------
0177   main_titles[32] = "Sum PF Charged Hadron Pt, #Delta R = 0.4";
0178   main_titles[33] = "Sum PF Neutral Hadron Pt, #Delta R = 0.4";
0179   main_titles[34] = "Sum PF Photon Et, #Delta R = 0.4";
0180   main_titles[35] = "Sum PF Neutral Hadron Pt (Higher Pt threshold), #Delta R = 0.4";
0181   main_titles[36] = "Sum PF Photon Et (Higher Pt threshold), #Delta R = 0.4";
0182   main_titles[37] = "Sum PF Charged Particles Pt not from PV  (for Pu corrections), #Delta R = 0.4";
0183 
0184   main_titles[38] = "Relative PF Isolation, #Delta R = 0.3";
0185   main_titles[39] = "Relative PF Isolation, #Delta R = 0.4";
0186 
0187   main_titles[40] = "Relative PF Isolation (Higher Pt threshold), #Delta R = 0.3";
0188   main_titles[41] = "Relative PF Isolation (Higher Pt threshold), #Delta R = 0.4";
0189 
0190   main_titles[42] = "Sum DR Isolation Profile for Charged Hadron,  #Delta R = 0.4";
0191 
0192   main_titles[43] = "Sum DR Isolation Profile for Neutral Hadron,  #Delta R = 0.4";
0193 
0194   main_titles[44] = "Sum DR Isolation Profile for Photon,  #Delta R = 0.4";
0195 
0196   main_titles[45] = "Mean DR Isolation Profile for Charged Hadron,  #Delta R = 0.4";
0197 
0198   main_titles[46] = "Mean DR Isolation Profile for Neutral Hadron,  #Delta R = 0.4";
0199 
0200   main_titles[47] = "Mean DR Isolation Profile for Photon,  #Delta R = 0.4";
0201 
0202 #ifdef DEBUG
0203   cout << "InitStatistics(): main titles 1D DONE " << endl;
0204 #endif
0205   titles_2D[0] = "Total Tracker Momentum, #Delta R = 0.3";
0206   titles_2D[1] = "Total EM Cal Energy, #Delta R = 0.3";
0207   titles_2D[2] = "Total Had Cal Energy, #Delta R = 0.3";
0208   titles_2D[3] = "Total HO Cal Energy, #Delta R = 0.3";
0209   titles_2D[4] = "Sum PF Charged Hadron Pt, #Delta R = 0.4";
0210   titles_2D[5] = "Sum PF Neutral Hadron Pt, #Delta R = 0.4";
0211   titles_2D[6] = "Sum PF Photon Et, #Delta R = 0.4";
0212   titles_2D[7] = "Sum PF Charged Pt Not from PV, #Delta R = 0.4";
0213   titles_2D[8] = "Relative Detector-Based Isolation, #Delta R = 0.4";
0214   titles_2D[9] = "Relative PF Isolation, #Delta R = 0.4";
0215 
0216   main_titles_NVtxs[0] = "Sum PF Neutral Hadron Pt, #DeltaR = 0.4 ( 20 < N_{Vtx} < 50)";
0217   main_titles_NVtxs[1] = "Sum PF Neutral Hadron Pt, #DeltaR = 0.4 (50 < N_{Vtx} < 80)";
0218   main_titles_NVtxs[2] = "Sum PF Neutral Hadron Pt, #DeltaR = 0.4 (80 < N_{Vtx})";
0219   main_titles_NVtxs[3] = "Sum PF Photon Et, #DeltaR = 0.4 ( 20 < N_{Vtx} < 50)";
0220   main_titles_NVtxs[4] = "Sum PF Photon Et, #DeltaR = 0.4 (50 < N_{Vtx} < 80)";
0221   main_titles_NVtxs[5] = "Sum PF Photon Et, #DeltaR = 0.4 (80 < N_{Vtx})";
0222 
0223 #ifdef DEBUG
0224   cout << "InitStatistics(): main titles 2D DONE " << endl;
0225 #endif
0226 
0227   //------Titles on the X or Y axis------------
0228   axis_titles[0] = "#Sigma p_{T}   (GeV)";
0229   axis_titles[1] = "#Sigma E_{T}^{EM}   (GeV)";
0230   axis_titles[2] = "#Sigma E_{T}^{Had}   (GeV)";
0231   axis_titles[3] = "#Sigma E_{T}^{HO}   (GeV)";
0232   axis_titles[4] = "N_{Tracks}";
0233   axis_titles[5] = "N_{Jets}";
0234   axis_titles[6] = "#Sigma p_{T,veto} (GeV)";
0235   axis_titles[7] = "#Sigma E_{T,veto}^{EM}   (GeV)";
0236   axis_titles[8] = "#Sigma E_{T,veto}^{Had}   (GeV)";
0237   axis_titles[9] = "#Sigma E_{T,veto}^{HO}   (GeV)";
0238   axis_titles[10] = "#Sigma p_{T} / N_{Tracks} (GeV)";
0239   axis_titles[11] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
0240 
0241   axis_titles[12] = "#Sigma p_{T}   (GeV)";
0242   axis_titles[13] = "#Sigma E_{T}^{EM}   (GeV)";
0243   axis_titles[14] = "#Sigma E_{T}^{Had}   (GeV)";
0244   axis_titles[15] = "#Sigma E_{T}^{HO}   (GeV)";
0245   axis_titles[16] = "N_{Tracks}";
0246   axis_titles[17] = "N_{Jets}";
0247   axis_titles[18] = "#Sigma p_{T,veto} (GeV)";
0248   axis_titles[19] = "#Sigma E_{T,veto}^{EM}   (GeV)";
0249   axis_titles[20] = "#Sigma E_{T,veto}^{Had}   (GeV)";
0250   axis_titles[21] = "#Sigma E_{T,veto}^{HO}   (GeV)";
0251   axis_titles[22] = "#Sigma p_{T} / N_{Tracks} (GeV)";
0252   axis_titles[23] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
0253 
0254   axis_titles[24] = "(#Sigma Tk p_{T} + #Sigma ECAL p_{T} + #Sigma HCAL p_{T})/ Mu p_{T}  (GeV)";
0255   axis_titles[25] = "(#Sigma Tk p_{T} + #Sigma ECAL p_{T} + #Sigma HCAL p_{T})/ Mu p_{T}  (GeV)";
0256 
0257   axis_titles[26] = "#Sigma PFCharged p_{T}";
0258   axis_titles[27] = "#Sigma PFNeutral p_{T}";
0259   axis_titles[28] = "#Sigma PFPhoton p_{T}";
0260   axis_titles[29] = "#Sigma PFNeutral p_{T}";
0261   axis_titles[30] = "#Sigma PFPhoton p_{T}";
0262   axis_titles[31] = "#Sigma PFCharged p_{T}";
0263 
0264   axis_titles[32] = "#Sigma PFCharged p_{T}";
0265   axis_titles[33] = "#Sigma PFNeutral p_{T}";
0266   axis_titles[34] = "#Sigma PFPhoton p_{T}";
0267   axis_titles[35] = "#Sigma PFNeutral p_{T}";
0268   axis_titles[36] = "#Sigma PFPhoton p_{T}";
0269   axis_titles[37] = "#Sigma PFCharged p_{T}";
0270 
0271   axis_titles[38] = "(#Sigma PFCharged p_{T} + #Sigma PFNeutral p_{T} + #Sigma PFPhoton p_{T}) Mu p_{T}  (GeV)";
0272   axis_titles[39] = "(#Sigma PFCharged p_{T} + #Sigma PFNeutral p_{T} + #Sigma PFPhoton p_{T}) Mu p_{T}  (GeV)";
0273   axis_titles[40] = "(#Sigma PFCharged p_{T} + #Sigma PFNeutral p_{T} + #Sigma PFPhoton p_{T}) Mu p_{T}  (GeV)";
0274   axis_titles[41] = "(#Sigma PFCharged p_{T} + #Sigma PFNeutral p_{T} + #Sigma PFPhoton p_{T}) Mu p_{T}  (GeV)";
0275 
0276   axis_titles[42] = "#Sigma DR PFCharged";
0277   axis_titles[43] = "#Sigma DR PFNeutral";
0278   axis_titles[44] = "#Sigma DR PFPhoton";
0279 
0280   axis_titles[45] = "Mean DR PFCharged";
0281   axis_titles[46] = "Mean DR PFNeutral";
0282   axis_titles[47] = "Mean DR PFPhoton";
0283 
0284   axis_titles_NVtxs[0] = "#Sigma PFNeutral p_{T}";
0285   axis_titles_NVtxs[1] = "#Sigma PFNeutral p_{T}";
0286   axis_titles_NVtxs[2] = "#Sigma PFNeutral p_{T}";
0287   axis_titles_NVtxs[3] = "#Sigma PFPhoton p_{T}";
0288   axis_titles_NVtxs[4] = "#Sigma PFPhoton p_{T}";
0289   axis_titles_NVtxs[5] = "#Sigma PFPhoton p_{T}";
0290 
0291 #ifdef DEBUG
0292   cout << "InitStatistics(): main titles 1D DONE " << endl;
0293 #endif
0294 
0295   //-----------Names given for the root file----------
0296   names[0] = "sumPt_R03";
0297   names[1] = "emEt_R03";
0298   names[2] = "hadEt_R03";
0299   names[3] = "hoEt_R03";
0300   names[4] = "nTracks_R03";
0301   names[5] = "nJets_R03";
0302   names[6] = "trackerVetoPt_R03";
0303   names[7] = "emVetoEt_R03";
0304   names[8] = "hadVetoEt_R03";
0305   names[9] = "hoVetoEt_R03";
0306   names[10] = "avgPt_R03";
0307   names[11] = "weightedEt_R03";
0308 
0309   names[12] = "sumPt_R05";
0310   names[13] = "emEt_R05";
0311   names[14] = "hadEt_R05";
0312   names[15] = "hoEt_R05";
0313   names[16] = "nTracks_R05";
0314   names[17] = "nJets_R05";
0315   names[18] = "trackerVetoPt_R05";
0316   names[19] = "emVetoEt_R05";
0317   names[20] = "hadVetoEt_R05";
0318   names[21] = "hoVetoEt_R05";
0319   names[22] = "avgPt_R05";
0320   names[23] = "weightedEt_R05";
0321 
0322   names[24] = "relDetIso_R03";
0323   names[25] = "relDetIso_R05";
0324 
0325   names[26] = "pfChargedPt_R03";
0326   names[27] = "pfNeutralPt_R03";
0327   names[28] = "pfPhotonPt_R03";
0328   names[29] = "pfNeutralPt_HT_R03";
0329   names[30] = "pfPhotonPt_HT_R03";
0330   names[31] = "pfChargedPt_PU_R03";
0331 
0332   names[32] = "pfChargedPt_R04";
0333   names[33] = "pfNeutralPt_R04";
0334   names[34] = "pfPhotonPt_R04";
0335   names[35] = "pfNeutralPt_HT_R04";
0336   names[36] = "pfPhotonPt_HT_R04";
0337   names[37] = "pfChargedPt_PU_R04";
0338 
0339   names[38] = "relPFIso_R03";
0340   names[39] = "relPFIso_R04";
0341 
0342   names[40] = "relPFIso_HT_R03";
0343   names[41] = "relPFIso_HT_R04";
0344 
0345   names[42] = "SumDR_PFCharged_R04";
0346   names[43] = "SumDR_PFNeutral_R04";
0347   names[44] = "SumDR_PFPhoton_R04";
0348 
0349   names[45] = "MeanDR_PFCharged_R04";
0350   names[46] = "MeanDR_PFNeutral_R04";
0351   names[47] = "MeanDR_PFPhoton_R04";
0352 
0353 #ifdef DEBUG
0354   cout << "InitStatistics(): names 1D DONE " << endl;
0355 #endif
0356 
0357   names_2D[0] = "SumPt_R03";
0358   names_2D[1] = "emEt_R03";
0359   names_2D[2] = "hadEt_R03";
0360   names_2D[3] = "hoEt_R03";
0361   names_2D[4] = "pfChargedPt_R04";
0362   names_2D[5] = "pfNeutralPt_R04";
0363   names_2D[6] = "pfPhotonPt_R04";
0364   names_2D[7] = "pfChargedPUPt_R04";
0365   names_2D[8] = "relDetIso_R03";
0366   names_2D[9] = "relPFIso_R04";
0367 
0368 #ifdef DEBUG
0369   cout << "InitStatistics(): names 2D DONE " << endl;
0370 #endif
0371 
0372   names_NVtxs[0] = "pfNeutralPt_R04_PV20to50";
0373   names_NVtxs[1] = "pfNeutralPt_R04_PV50to80";
0374   names_NVtxs[2] = "pfNeutralPt_R04_PV80toInf";
0375   names_NVtxs[3] = "pfPhotonPt_R04_PV20to50";
0376   names_NVtxs[4] = "pfPhotonPt_R04_PV50to80";
0377   names_NVtxs[5] = "pfPhotonPt_R04_PV80toInf";
0378 
0379   //----------Parameters for binning of histograms---------
0380   //param[var][0] is the number of bins
0381   //param[var][1] is the low edge of the low bin
0382   //param[var][2] is the high edge of the high bin
0383   //
0384   // maximum value------,
0385   //                    |
0386   //                    V
0387   param[0][0] = (int)(20.0 / S_BIN_WIDTH);
0388   param[0][1] = 0.0;
0389   param[0][2] = param[0][0] * S_BIN_WIDTH;
0390   param[1][0] = (int)(20.0 / S_BIN_WIDTH);
0391   param[1][1] = 0.0;
0392   param[1][2] = param[1][0] * S_BIN_WIDTH;
0393   param[2][0] = (int)(20.0 / S_BIN_WIDTH);
0394   param[2][1] = 0.0;
0395   param[2][2] = param[2][0] * S_BIN_WIDTH;
0396   param[3][0] = 20;
0397   param[3][1] = 0.0;
0398   param[3][2] = 2.0;
0399   param[4][0] = 16;
0400   param[4][1] = -0.5;
0401   param[4][2] = param[4][0] - 0.5;
0402   param[5][0] = 4;
0403   param[5][1] = -0.5;
0404   param[5][2] = param[5][0] - 0.5;
0405   param[6][0] = (int)(40.0 / S_BIN_WIDTH);
0406   param[6][1] = 0.0;
0407   param[6][2] = param[6][0] * S_BIN_WIDTH;
0408   param[7][0] = 20;
0409   param[7][1] = 0.0;
0410   param[7][2] = 10.0;
0411   param[8][0] = (int)(20.0 / S_BIN_WIDTH);
0412   param[8][1] = 0.0;
0413   param[8][2] = param[8][0] * S_BIN_WIDTH;
0414   param[9][0] = 20;
0415   param[9][1] = 0.0;
0416   param[9][2] = 5.0;
0417   param[10][0] = (int)(15.0 / S_BIN_WIDTH);
0418   param[10][1] = 0.0;
0419   param[10][2] = param[10][0] * S_BIN_WIDTH;
0420   param[11][0] = (int)(20.0 / S_BIN_WIDTH);
0421   param[11][1] = 0.0;
0422   param[11][2] = param[11][0] * S_BIN_WIDTH;
0423 
0424   param[12][0] = (int)(20.0 / S_BIN_WIDTH);
0425   param[12][1] = 0.0;
0426   param[12][2] = param[12][0] * S_BIN_WIDTH;
0427   param[13][0] = (int)(20.0 / S_BIN_WIDTH);
0428   param[13][1] = 0.0;
0429   param[13][2] = param[13][0] * S_BIN_WIDTH;
0430   param[14][0] = (int)(20.0 / S_BIN_WIDTH);
0431   param[14][1] = 0.0;
0432   param[14][2] = param[14][0] * S_BIN_WIDTH;
0433   param[15][0] = 20;
0434   param[15][1] = 0.0;
0435   param[15][2] = 2.0;
0436   param[16][0] = 16;
0437   param[16][1] = -0.5;
0438   param[16][2] = param[16][0] - 0.5;
0439   param[17][0] = 4;
0440   param[17][1] = -0.5;
0441   param[17][2] = param[17][0] - 0.5;
0442   param[18][0] = (int)(40.0 / S_BIN_WIDTH);
0443   param[18][1] = 0.0;
0444   param[18][2] = param[18][0] * S_BIN_WIDTH;
0445   param[19][0] = 20;
0446   param[19][1] = 0.0;
0447   param[19][2] = 10.0;
0448   param[20][0] = (int)(20.0 / S_BIN_WIDTH);
0449   param[20][1] = 0.0;
0450   param[20][2] = param[20][0] * S_BIN_WIDTH;
0451   param[21][0] = 20;
0452   param[21][1] = 0.0;
0453   param[21][2] = 5.0;
0454   param[22][0] = (int)(15.0 / S_BIN_WIDTH);
0455   param[22][1] = 0.0;
0456   param[22][2] = param[22][0] * S_BIN_WIDTH;
0457   param[23][0] = (int)(20.0 / S_BIN_WIDTH);
0458   param[23][1] = 0.0;
0459   param[23][2] = param[23][0] * S_BIN_WIDTH;
0460 
0461   param[24][0] = 50;
0462   param[24][1] = 0.0;
0463   param[24][2] = 1.0;
0464   param[25][0] = 50;
0465   param[25][1] = 0.0;
0466   param[25][2] = 1.0;
0467 
0468   param[26][0] = (int)(20.0 / S_BIN_WIDTH);
0469   param[26][1] = 0.0;
0470   param[26][2] = param[26][0] * S_BIN_WIDTH;
0471   param[27][0] = (int)(20.0 / S_BIN_WIDTH);
0472   param[27][1] = 0.0;
0473   param[27][2] = param[27][0] * S_BIN_WIDTH;
0474   param[28][0] = (int)(20.0 / S_BIN_WIDTH);
0475   param[28][1] = 0.0;
0476   param[28][2] = param[28][0] * S_BIN_WIDTH;
0477   param[29][0] = (int)(20.0 / S_BIN_WIDTH);
0478   param[29][1] = 0.0;
0479   param[29][2] = param[29][0] * S_BIN_WIDTH;
0480   param[30][0] = (int)(20.0 / S_BIN_WIDTH);
0481   param[30][1] = 0.0;
0482   param[30][2] = param[30][0] * S_BIN_WIDTH;
0483   param[31][0] = (int)(20.0 / S_BIN_WIDTH);
0484   param[31][1] = 0.0;
0485   param[31][2] = param[31][0] * S_BIN_WIDTH;
0486 
0487   param[32][0] = (int)(20.0 / S_BIN_WIDTH);
0488   param[32][1] = 0.0;
0489   param[32][2] = param[32][0] * S_BIN_WIDTH;
0490   param[33][0] = (int)(20.0 / S_BIN_WIDTH);
0491   param[33][1] = 0.0;
0492   param[33][2] = param[33][0] * S_BIN_WIDTH;
0493   param[34][0] = (int)(20.0 / S_BIN_WIDTH);
0494   param[34][1] = 0.0;
0495   param[34][2] = param[34][0] * S_BIN_WIDTH;
0496   param[35][0] = (int)(20.0 / S_BIN_WIDTH);
0497   param[35][1] = 0.0;
0498   param[35][2] = param[35][0] * S_BIN_WIDTH;
0499   param[36][0] = (int)(20.0 / S_BIN_WIDTH);
0500   param[36][1] = 0.0;
0501   param[36][2] = param[36][0] * S_BIN_WIDTH;
0502   param[37][0] = (int)(20.0 / S_BIN_WIDTH);
0503   param[37][1] = 0.0;
0504   param[37][2] = param[37][0] * S_BIN_WIDTH;
0505 
0506   param[38][0] = 50;
0507   param[38][1] = 0.0;
0508   param[38][2] = 1.0;
0509   param[39][0] = 50;
0510   param[39][1] = 0.0;
0511   param[39][2] = 1.0;
0512 
0513   param[40][0] = 50;
0514   param[40][1] = 0.0;
0515   param[40][2] = 1.0;
0516   param[41][0] = 50;
0517   param[41][1] = 0.0;
0518   param[41][2] = 1.0;
0519 
0520   param[42][0] = 50;
0521   param[42][1] = 0.0;
0522   param[42][2] = 5;
0523   param[43][0] = 50;
0524   param[43][1] = 0.0;
0525   param[43][2] = 5;
0526   param[44][0] = 50;
0527   param[44][1] = 0.0;
0528   param[44][2] = 5;
0529 
0530   param[45][0] = 50;
0531   param[45][1] = 0.0;
0532   param[45][2] = 0.4;
0533   param[46][0] = 50;
0534   param[46][1] = 0.0;
0535   param[46][2] = 0.4;
0536   param[47][0] = 50;
0537   param[47][1] = 0.0;
0538   param[47][2] = 0.4;
0539 
0540   //--------------Is the variable continuous (i.e. non-integer)?-------------
0541   //---------(Log binning will only be used for continuous variables)--------
0542   isContinuous[0] = 1;
0543   isContinuous[1] = 1;
0544   isContinuous[2] = 1;
0545   isContinuous[3] = 1;
0546   isContinuous[4] = 0;
0547   isContinuous[5] = 0;
0548   isContinuous[6] = 1;
0549   isContinuous[7] = 1;
0550   isContinuous[8] = 1;
0551   isContinuous[9] = 1;
0552   isContinuous[10] = 1;
0553   isContinuous[11] = 1;
0554 
0555   isContinuous[12] = 1;
0556   isContinuous[13] = 1;
0557   isContinuous[14] = 1;
0558   isContinuous[15] = 1;
0559   isContinuous[16] = 0;
0560   isContinuous[17] = 0;
0561   isContinuous[18] = 1;
0562   isContinuous[19] = 1;
0563   isContinuous[20] = 1;
0564   isContinuous[21] = 1;
0565   isContinuous[22] = 1;
0566   isContinuous[23] = 1;
0567 
0568   isContinuous[24] = 1;
0569   isContinuous[25] = 1;
0570   isContinuous[26] = 1;
0571   isContinuous[27] = 1;
0572   isContinuous[28] = 1;
0573   isContinuous[29] = 1;
0574   isContinuous[30] = 1;
0575   isContinuous[31] = 1;
0576   isContinuous[32] = 1;
0577   isContinuous[33] = 1;
0578   isContinuous[34] = 1;
0579   isContinuous[35] = 1;
0580   isContinuous[36] = 1;
0581   isContinuous[37] = 1;
0582   isContinuous[38] = 1;
0583   isContinuous[39] = 1;
0584   isContinuous[40] = 1;
0585   isContinuous[41] = 1;
0586   isContinuous[42] = 1;
0587   isContinuous[43] = 1;
0588   isContinuous[44] = 1;
0589   isContinuous[45] = 1;
0590   isContinuous[46] = 1;
0591   isContinuous[47] = 1;
0592 
0593 #ifdef DEBUG
0594   cout << "InitStatistics(): DONE " << endl;
0595 #endif
0596 }
0597 
0598 // ------------ method called for each event  ------------
0599 void MuonIsolationDQM::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0600   ++nEvents;
0601   edm::LogInfo("Tutorial") << "\nInvestigating event #" << nEvents << "\n";
0602 #ifdef DEBUG
0603   cout << "[MuonIsolationDQM]: analyze()" << endl;
0604 #endif
0605 
0606   // Get Muon Collection
0607   edm::Handle<edm::View<reco::Muon> > muons;
0608   iEvent.getByToken(theMuonCollectionLabel_, muons);
0609 
0610 #ifdef DEBUG
0611   cout << "[MuonIsolationDQM]: Number of muons -> " << muons->size() << endl;
0612 #endif
0613 
0614   int theMuonData = muons->size();
0615   h_nMuons->Fill(theMuonData);
0616 #ifdef DEBUG
0617   cout << "[MuonIsolationDQM]: Vertex is Valid" << endl;
0618 #endif
0619 
0620   //Get Vertex Information
0621   int _numPV = 0;
0622   edm::Handle<reco::VertexCollection> vertexHandle;
0623   iEvent.getByToken(theVertexCollectionLabel_, vertexHandle);
0624 
0625   if (vertexHandle.isValid()) {
0626     reco::VertexCollection vertex = *(vertexHandle.product());
0627     for (reco::VertexCollection::const_iterator v = vertex.begin(); v != vertex.end(); ++v) {
0628       if (v->isFake())
0629         continue;
0630       if (v->ndof() < 4)
0631         continue;
0632       if (fabs(v->z()) > 24.0)
0633         continue;
0634       ++_numPV;
0635     }
0636   }
0637 
0638 #ifdef DEBUG
0639   cout << "[MuonIsolationDQM]: Vertex is Valid" << endl;
0640 #endif
0641   // Get Muon Collection
0642 
0643   //Fill historgams concerning muon isolation
0644   for (edm::View<reco::Muon>::const_iterator muon = muons->begin(); muon != muons->end(); ++muon) {
0645     if (requireSTAMuon && muon->isStandAloneMuon()) {
0646       ++nSTAMuons;
0647       RecordData(*muon);
0648       FillHistos(_numPV);
0649     } else if (requireTRKMuon && muon->isTrackerMuon()) {
0650       ++nTRKMuons;
0651       RecordData(*muon);
0652       FillHistos(_numPV);
0653     } else if (requireGLBMuon && muon->isGlobalMuon()) {
0654       ++nGLBMuons;
0655       RecordData(*muon);
0656       FillHistos(_numPV);
0657       FillNVtxHistos(_numPV);
0658     }
0659   }
0660 }
0661 
0662 //---------------Record data for a signle muon's data---------------------
0663 void MuonIsolationDQM::RecordData(const reco::Muon& muon) {
0664 #ifdef DEBUG
0665   std::cout << "RecordData()" << endl;
0666 #endif
0667   float MuPt = muon.pt();
0668 
0669   theData[0] = muon.isolationR03().sumPt;
0670   theData[1] = muon.isolationR03().emEt;
0671   theData[2] = muon.isolationR03().hadEt;
0672   theData[3] = muon.isolationR03().hoEt;
0673 
0674   theData[4] = muon.isolationR03().nTracks;
0675   theData[5] = muon.isolationR03().nJets;
0676   theData[6] = muon.isolationR03().trackerVetoPt;
0677   theData[7] = muon.isolationR03().emVetoEt;
0678   theData[8] = muon.isolationR03().hadVetoEt;
0679   theData[9] = muon.isolationR03().hoVetoEt;
0680 
0681   // make sure nTracks != 0 before filling this one
0682   if (theData[4] != 0)
0683     theData[10] = (double)theData[0] / (double)theData[4];
0684   else
0685     theData[10] = -99;
0686 
0687   theData[11] = 1.5 * theData[1] + theData[2];
0688 
0689   theData[12] = muon.isolationR05().sumPt;
0690   theData[13] = muon.isolationR05().emEt;
0691   theData[14] = muon.isolationR05().hadEt;
0692   theData[15] = muon.isolationR05().hoEt;
0693 
0694   theData[16] = muon.isolationR05().nTracks;
0695   theData[17] = muon.isolationR05().nJets;
0696   theData[18] = muon.isolationR05().trackerVetoPt;
0697   theData[19] = muon.isolationR05().emVetoEt;
0698   theData[20] = muon.isolationR05().hadVetoEt;
0699   theData[21] = muon.isolationR05().hoVetoEt;
0700 
0701   // make sure nTracks != 0 before filling this one
0702   if (theData[16] != 0)
0703     theData[22] = (double)theData[12] / (double)theData[16];
0704   else
0705     theData[22] = -99;
0706 
0707   theData[23] = 1.5 * theData[13] + theData[14];
0708 
0709   theData[24] = (theData[0] + theData[1] + theData[2]) / MuPt;
0710   theData[25] = (theData[12] + theData[13] + theData[14]) / MuPt;
0711 
0712   theData[26] = muon.pfIsolationR03().sumChargedHadronPt;
0713   theData[27] = muon.pfIsolationR03().sumNeutralHadronEt;
0714   theData[28] = muon.pfIsolationR03().sumPhotonEt;
0715   theData[29] = muon.pfIsolationR03().sumNeutralHadronEtHighThreshold;
0716   theData[30] = muon.pfIsolationR03().sumPhotonEtHighThreshold;
0717   theData[31] = muon.pfIsolationR03().sumPUPt;
0718 
0719   theData[32] = muon.pfIsolationR04().sumChargedHadronPt;
0720   theData[33] = muon.pfIsolationR04().sumNeutralHadronEt;
0721   theData[34] = muon.pfIsolationR04().sumPhotonEt;
0722   theData[35] = muon.pfIsolationR04().sumNeutralHadronEtHighThreshold;
0723   theData[36] = muon.pfIsolationR04().sumPhotonEtHighThreshold;
0724   theData[37] = muon.pfIsolationR04().sumPUPt;
0725 
0726   theData[38] = (theData[26] + theData[27] + theData[28]) / MuPt;
0727   theData[39] = (theData[32] + theData[33] + theData[34]) / MuPt;
0728 
0729   theData[40] = (theData[26] + theData[29] + theData[30]) / MuPt;
0730   theData[41] = (theData[32] + theData[35] + theData[36]) / MuPt;
0731 
0732   theData[42] = muon.pfSumDRIsoProfileR04().sumChargedHadronPt;
0733   theData[43] = muon.pfSumDRIsoProfileR04().sumNeutralHadronEt;
0734   theData[44] = muon.pfSumDRIsoProfileR04().sumPhotonEt;
0735   theData[45] = muon.pfMeanDRIsoProfileR04().sumChargedHadronPt;
0736   theData[46] = muon.pfMeanDRIsoProfileR04().sumNeutralHadronEt;
0737   theData[47] = muon.pfMeanDRIsoProfileR04().sumPhotonEt;
0738 
0739   //--------------Filling the 2D Histos Data -------- //
0740   theData2D[0] = muon.isolationR03().sumPt;
0741   theData2D[1] = muon.isolationR03().emEt;
0742   theData2D[2] = muon.isolationR03().hadEt;
0743   theData2D[3] = muon.isolationR03().hoEt;
0744 
0745   theData2D[4] = muon.pfIsolationR04().sumChargedHadronPt;
0746   theData2D[5] = muon.pfIsolationR04().sumNeutralHadronEt;
0747   theData2D[6] = muon.pfIsolationR04().sumPhotonEt;
0748   theData2D[7] = muon.pfIsolationR04().sumPUPt;
0749 
0750   theData2D[8] = theData2D[0] + theData2D[1] + theData2D[2] + theData2D[3] / MuPt;  //Det RelIso;
0751   theData2D[9] = theData2D[4] + theData2D[5] + theData2D[6] / MuPt;                 //PF  RelIso;
0752 
0753   //-----------Filling the NVTX 1D HISTOS DATA ------------- //
0754   theDataNVtx[0] = muon.pfIsolationR04().sumNeutralHadronEt;
0755   theDataNVtx[1] = theDataNVtx[0];
0756   theDataNVtx[2] = theDataNVtx[0];
0757 
0758   theDataNVtx[3] = muon.pfIsolationR04().sumPhotonEt;
0759   theDataNVtx[4] = theDataNVtx[3];
0760   theDataNVtx[5] = theDataNVtx[3];
0761 }
0762 void MuonIsolationDQM::bookHistograms(DQMStore::IBooker& ibooker,
0763                                       edm::Run const& /*iRun*/,
0764                                       edm::EventSetup const& /* iSetup */) {
0765   ibooker.cd();
0766   ibooker.setCurrentFolder(dirName);
0767 
0768   ibooker.cd();
0769   ibooker.setCurrentFolder(dirName);
0770 
0771   //---initialize number of muons histogram---
0772   h_nMuons = ibooker.book1D("nMuons", title_sam + "Number of Muons", 20, 0., 20.);
0773   h_nMuons->setAxisTitle("Number of Muons", XAXIS);
0774   h_nMuons->setAxisTitle("Fraction of Events", YAXIS);
0775 
0776   //---Initialize 1D Histograms---
0777   for (int var = 0; var < NUM_VARS; var++) {
0778     h_1D[var] = ibooker.book1D(
0779         names[var], title_sam + main_titles[var] + title_cone, (int)param[var][0], param[var][1], param[var][2]);
0780     h_1D[var]->setAxisTitle(axis_titles[var], XAXIS);
0781     //    GetTH1FromMonitorElement(h_1D[var])->Sumw2();
0782   }  //Finish 1D
0783 
0784   //----Initialize 2D Histograms
0785   for (int var = 0; var < NUM_VARS_2D; var++) {
0786     h_2D[var] = ibooker.bookProfile(
0787         names_2D[var] + "_VsPV", titles_2D[var] + " Vs PV", vtxBin_, vtxMin_, vtxMax_, 20, 0.0, 20.0);
0788     h_2D[var]->setAxisTitle("Number of PV", XAXIS);
0789     h_2D[var]->setAxisTitle(titles_2D[var] + " (GeV)", YAXIS);
0790     //    h_2D[var]->getTH1()->Sumw2();
0791   }
0792 
0793   //-----Initialise PU-Binned histograms
0794   for (int var = 0; var < NUM_VARS_NVTX; var++) {
0795     h_1D_NVTX[var] = ibooker.book1D(names_NVtxs[var], main_titles_NVtxs[var], 50, 0.0, 10.0);
0796     h_1D_NVTX[var]->setAxisTitle(axis_titles_NVtxs[var], XAXIS);
0797     ///    GetTH1FromMonitorElement(h_1D_NVTX[var])->Sumw2();
0798   }
0799 }
0800 
0801 void MuonIsolationDQM::NormalizeHistos() {
0802   for (int var = 0; var < NUM_VARS; var++) {
0803     double entries = GetTH1FromMonitorElement(h_1D[var])->GetEntries();
0804     GetTH1FromMonitorElement(h_1D[var])->Scale(1. / entries);
0805   }
0806 }
0807 
0808 void MuonIsolationDQM::FillHistos(int numPV) {
0809 #ifdef DEBUG
0810   cout << "FillHistos( " << numPV << " )" << endl;
0811 #endif
0812 
0813   //----------Fill 1D histograms---------------
0814   for (int var = 0; var < NUM_VARS; var++) {
0815     h_1D[var]->Fill(theData[var]);
0816     //    cd_plots[var]->Fill(theData[var]);//right now, this is a regular PDF (just like h_1D)
0817     //OFBin   if (theData[var] > param[var][2]) {
0818     //OFBin     // fill the overflow bin
0819     //OFBin     overFlowBin = (int) param[var][0] + 1;
0820     //OFBin     overFlow = GetTH1FromMonitorElement(h_1D[var])->GetBinContent(overFlowBin);
0821     //OFBin     GetTH1FromMonitorElement(h_1D[var])->SetBinContent(overFlowBin, overFlow + 1);
0822     //OFBin   }
0823   }  //Finish 1D
0824 
0825   for (int var = 0; var < NUM_VARS_2D; var++) {
0826     h_2D[var]->Fill(numPV, theData2D[var]);
0827   }
0828 
0829 #ifdef DEBUG
0830   cout << "FillHistos( " << numPV << " ): DONE" << endl;
0831 #endif
0832 }
0833 void MuonIsolationDQM::FillNVtxHistos(int PV) {
0834   if (PV >= 20 && PV < 50) {
0835     h_1D_NVTX[0]->Fill(theDataNVtx[0]);
0836     h_1D_NVTX[3]->Fill(theDataNVtx[3]);
0837   }
0838   if (PV >= 50 && PV < 80) {
0839     h_1D_NVTX[1]->Fill(theDataNVtx[1]);
0840     h_1D_NVTX[4]->Fill(theDataNVtx[4]);
0841   }
0842   if (PV >= 80) {
0843     h_1D_NVTX[2]->Fill(theDataNVtx[2]);
0844     h_1D_NVTX[5]->Fill(theDataNVtx[5]);
0845   }
0846 }
0847 
0848 TH1* MuonIsolationDQM::GetTH1FromMonitorElement(MonitorElement* me) { return me->getTH1(); }
0849 
0850 //define this as a plug-in
0851 DEFINE_FWK_MODULE(MuonIsolationDQM);
[ 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!