Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 // MuIsoValidation.cc
00003 // Package:    Muon Isolation Validation
00004 // Class:      MuIsoValidation
00005 // 
00006 /*
00007 
00008 
00009  Description: Muon Isolation Validation class
00010 
00011  Implementation: This code will accept a data set and generate plots of
00012         various quantities relevent to the Muon Isolation module. We will 
00013         be using the IsoDeposit class, *not* the MuonIsolation struct.
00014          
00015         The sequence of events is... 
00016                 * initalize statics (which variables to plot, axis limtis, etc.)
00017                 * run over events
00018                         * for each event, run over the muons in that event
00019                                 *record IsoDeposit data
00020                 * transfer data to histograms, profile plots
00021                 * save histograms to a root file
00022                 
00023         Easy-peasy.
00024         
00025 */
00026 //
00027 // Original Author:  "C. Jess Riedel", UC Santa Barbara
00028 //         Created:  Tue Jul 17 15:58:24 CDT 2007
00029 //
00030 
00031 //Class header file
00032 #include "Validation/MuonIsolation/interface/MuIsoValidation.h"
00033 
00034 //System included files
00035 #include <memory>
00036 #include <string>
00037 #include <typeinfo>
00038 #include <utility>
00039 
00040 //Root included files
00041 #include "TH1.h"
00042 #include "TH2.h"
00043 #include "TProfile.h"
00044 
00045 //Event framework included files
00046 #include "FWCore/Framework/interface/MakerMacros.h"
00047 #include "FWCore/Framework/interface/Event.h"
00048 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00049 
00050 //Other included files
00051 #include "DataFormats/TrackReco/interface/Track.h"
00052 
00053 //Using declarations
00054 using std::vector;
00055 using std::pair;
00056 using std::string;
00057 
00058 
00059 
00060 //
00061 //-----------------Constructors---------------------
00062 //
00063 MuIsoValidation::MuIsoValidation(const edm::ParameterSet& iConfig)
00064 {
00065   
00066   //  rootfilename = iConfig.getUntrackedParameter<string>("rootfilename"); // comment out for inclusion
00067   requireCombinedMuon = iConfig.getUntrackedParameter<bool>("requireCombinedMuon");
00068   dirName = iConfig.getParameter<std::string>("directory");
00069   //  subDirName = iConfig.getParameter<std::string>("@module_label");
00070   
00071   //  dirName += subDirName;
00072   
00073   //--------Initialize tags-------
00074   Muon_Tag = iConfig.getUntrackedParameter<edm::InputTag>("Global_Muon_Label");
00075   
00076   //-------Initialize counters----------------
00077   nEvents = 0;
00078   nIncMuons = 0;   
00079   //  nCombinedMuons = 0;
00080   
00081   InitStatics();
00082   
00083   //Set up DAQ
00084   dbe = 0;
00085   dbe = edm::Service<DQMStore>().operator->();
00086   
00087   //------"allocate" space for the data vectors-------
00088   
00089   /*
00090     h_1D        is a 2D vector with indices [var][muon#]
00091     cd_plots    is a 2D vector with indices [var][muon#]  
00092     h_2D        is a 3D vector with indices [var][var][muon#]
00093     p_2D        is a 3D vector with indices [var][var][muon#]
00094   */
00095   //NOTE:the total number of muons and events is initially unknown, 
00096   //       so that dimension is not initialized. Hence, theMuonData
00097   //     needs no resizing.
00098   
00099   h_1D.resize    (NUM_VARS);
00100   cd_plots.resize(NUM_VARS);
00101   //  h_2D.resize(NUM_VARS, vector<MonitorElement*>     (NUM_VARS));
00102   p_2D.resize(NUM_VARS, vector<MonitorElement*>(NUM_VARS));
00103   
00104   dbe->cd();
00105 }
00106 
00107 //
00108 //----------Destructor-----------------
00109 //
00110 MuIsoValidation::~MuIsoValidation(){
00111   
00112   //Deallocate memory
00113   
00114 }
00115 
00116 //
00117 //------------Methods-----------------------
00118 //
00119 void MuIsoValidation::InitStatics(){
00120   
00121   //-----------Initialize primatives-----------
00122   S_BIN_WIDTH = 1.0;//in GeV
00123   L_BIN_WIDTH = 2.0;//in GeV
00124   LOG_BINNING_ENABLED = 1;
00125   NUM_LOG_BINS = 15;
00126   LOG_BINNING_RATIO = 1.1;//ratio by which each bin is wider than the last for log binning
00127   //i.e.  bin widths are (x), (r*x), (r^2*x), ..., (r^(nbins)*x)
00128   
00129   
00130   //-------Initialize Titles---------
00131   title_sam = "";//"[Sample b-jet events] ";
00132   title_cone = "";//" [in R=0.3 IsoDeposit Cone]";
00133   //The above two pieces of info will be printed on the title of the whole page,
00134   //not for each individual histogram
00135   title_cd = "C.D. of ";
00136   
00137   //-------"Allocate" memory for vectors
00138   main_titles.resize(NUM_VARS);
00139   axis_titles.resize(NUM_VARS);
00140   names.resize(NUM_VARS);
00141   param.resize(NUM_VARS, vector<double>(3) );
00142   isContinuous.resize(NUM_VARS);
00143   cdCompNeeded.resize(NUM_VARS);
00144   
00145   //-----Titles of the plots-----------
00146   main_titles[0 ] = "Total Tracker Momentum";
00147   main_titles[1 ] = "Total EM Cal Energy";
00148   main_titles[2 ] = "Total Had Cal Energy";
00149   main_titles[3 ] = "Total HO Cal Energy";
00150   main_titles[4 ] = "Number of Tracker Tracks";
00151   main_titles[5 ] = "Number of Jets around Muon";
00152   main_titles[6 ] = "Tracker p_{T} within veto cone";
00153   main_titles[7 ] = "EM E_{T} within veto cone";
00154   main_titles[8 ] = "Had E_{T} within veto cone";
00155   main_titles[9 ] = "HO E_{T} within veto cone";
00156   main_titles[10] = "Muon p_{T}";
00157   main_titles[11] = "Muon #eta";
00158   main_titles[12] = "Muon #phi";
00159   main_titles[13] = "Average Momentum per Track ";
00160   main_titles[14] = "Weighted Energy";
00161   main_titles[15] = "PF Sum of Charged Hadron Pt";
00162   main_titles[16] = "PF Sum of Total Hadron Pt";
00163   main_titles[17] = "PF Sum of E,Mu Pt";
00164   main_titles[18] = "PF Sum of Neutral Hadron Et";
00165   main_titles[19] = "PF Sum of Photon Et";
00166   main_titles[20] = "PF Sum of Pt from non-PV";
00167 
00168   //------Titles on the X or Y axis------------
00169   axis_titles[0 ] = "#Sigma p_{T}   (GeV)";
00170   axis_titles[1 ] = "#Sigma E_{T}^{EM}   (GeV)";
00171   axis_titles[2 ] = "#Sigma E_{T}^{Had}   (GeV)";
00172   axis_titles[3 ] = "#Sigma E_{T}^{HO}   (GeV)";
00173   axis_titles[4 ] = "N_{Tracks}";
00174   axis_titles[5 ] = "N_{Jets}";
00175   axis_titles[6 ] = "#Sigma p_{T,veto} (GeV)";
00176   axis_titles[7 ] = "#Sigma E_{T,veto}^{EM}   (GeV)";
00177   axis_titles[8 ] = "#Sigma E_{T,veto}^{Had}   (GeV)";
00178   axis_titles[9 ] = "#Sigma E_{T,veto}^{HO}   (GeV)";
00179   axis_titles[10] = "p_{T,#mu} (GeV)";
00180   axis_titles[11] = "#eta_{#mu}";
00181   axis_titles[12] = "#phi_{#mu}";
00182   axis_titles[13] = "#Sigma p_{T} / N_{Tracks} (GeV)";
00183   axis_titles[14] = "(1.5) X #Sigma E_{T}^{EM} + #Sigma E_{T}^{Had}";
00184   axis_titles[15] = "#Sigma p_{T}^{PFHadCha}   (GeV)";
00185   axis_titles[16] = "#Sigma p_{T}^{PFTotCha}   (GeV)";
00186   axis_titles[17] = "#Sigma p_{T}^{PFEMu}   (GeV)";
00187   axis_titles[18] = "#Sigma E_{T}^{PFHadNeu}   (GeV)";
00188   axis_titles[19] = "#Sigma E_{T}^{PFPhot}   (GeV)";
00189   axis_titles[20] = "#Sigma p_{T}^{PFPU}   (GeV)";
00190 
00191   //-----------Names given for the root file----------
00192   names[0 ] = "sumPt";
00193   names[1 ] = "emEt";
00194   names[2 ] = "hadEt";
00195   names[3 ] = "hoEt";
00196   names[4 ] = "nTracks";
00197   names[5 ] = "nJets";
00198   names[6 ] = "trackerVetoPt";
00199   names[7 ] = "emVetoEt";
00200   names[8 ] = "hadVetoEt";
00201   names[9 ] = "hoVetoEt";
00202   names[10] = "muonPt";
00203   names[11] = "muonEta";
00204   names[12] = "muonPhi";
00205   names[13] = "avgPt";
00206   names[14] = "weightedEt";
00207   names[15] = "PFsumChargedHadronPt";
00208   names[16] = "PFsumChargedTotalPt";
00209   names[17] = "PFsumEMuPt";
00210   names[18] = "PFsumNeutralHadronEt";
00211   names[19] = "PFsumPhotonEt";
00212   names[20] = "PFsumPUPt";
00213 
00214   //----------Parameters for binning of histograms---------
00215   //param[var][0] is the number of bins
00216   //param[var][1] is the low edge of the low bin
00217   //param[var][2] is the high edge of the high bin
00218   //
00219   // maximum value------,
00220   //                    |
00221   //                    V                  
00222   param[0 ][0]= (int)( 20.0/S_BIN_WIDTH); param[0 ][1]=  0.0; param[0 ][2]= param[0 ][0]*S_BIN_WIDTH;
00223   param[1 ][0]= (int)( 20.0/S_BIN_WIDTH); param[1 ][1]=  0.0; param[1 ][2]= param[1 ][0]*S_BIN_WIDTH;
00224   param[2 ][0]= (int)( 20.0/S_BIN_WIDTH); param[2 ][1]=  0.0; param[2 ][2]= param[2 ][0]*S_BIN_WIDTH;
00225   param[3 ][0]=                       20; param[3 ][1]=  0.0; param[3 ][2]=                      2.0;
00226   param[4 ][0]=                       16; param[4 ][1]= -0.5; param[4 ][2]=         param[4 ][0]-0.5;
00227   param[5 ][0]=                        4; param[5 ][1]= -0.5; param[5 ][2]=         param[5 ][0]-0.5;
00228   param[6 ][0]= (int)( 40.0/S_BIN_WIDTH); param[6 ][1]=  0.0; param[6 ][2]= param[6 ][0]*S_BIN_WIDTH;
00229   param[7 ][0]=                       20; param[7 ][1]=  0.0; param[7 ][2]=                     10.0;
00230   param[8 ][0]= (int)( 20.0/S_BIN_WIDTH); param[8 ][1]=  0.0; param[8 ][2]= param[8 ][0]*S_BIN_WIDTH;
00231   param[9 ][0]=                       20; param[9 ][1]=  0.0; param[9 ][2]=                      5.0;
00232   param[10][0]= (int)( 40.0/S_BIN_WIDTH); param[10][1]=  0.0; param[10][2]= param[10][0]*S_BIN_WIDTH;
00233   param[11][0]=                       24; param[11][1]= -2.4; param[11][2]=                      2.4;
00234   param[12][0]=                       32; param[12][1]= -3.2; param[12][2]=                      3.2;
00235   param[13][0]= (int)( 15.0/S_BIN_WIDTH); param[13][1]=  0.0; param[13][2]= param[13][0]*S_BIN_WIDTH;
00236   param[14][0]= (int)( 20.0/S_BIN_WIDTH); param[14][1]=  0.0; param[14][2]= param[14][0]*S_BIN_WIDTH;
00237   param[15][0]= (int)( 20.0/S_BIN_WIDTH); param[15][1]=  0.0; param[15][2]= param[15][0]*S_BIN_WIDTH;
00238   param[16][0]= (int)( 20.0/S_BIN_WIDTH); param[15][1]=  0.0; param[16][2]= param[16][0]*S_BIN_WIDTH;
00239   param[17][0]= (int)( 20.0/S_BIN_WIDTH)+1; param[17][1]= -S_BIN_WIDTH; param[17][2]= param[17][0]*S_BIN_WIDTH;
00240   param[18][0]= (int)( 20.0/S_BIN_WIDTH); param[18][1]=  0.0; param[18][2]= param[18][0]*S_BIN_WIDTH;
00241   param[19][0]= (int)( 20.0/S_BIN_WIDTH); param[19][1]=  0.0; param[19][2]= param[19][0]*S_BIN_WIDTH;
00242   param[20][0]= (int)( 20.0/S_BIN_WIDTH); param[20][1]=  0.0; param[20][2]= param[20][0]*S_BIN_WIDTH;
00243 
00244   //--------------Is the variable continuous (i.e. non-integer)?-------------
00245   //---------(Log binning will only be used for continuous variables)--------
00246   isContinuous[0 ] = 1;
00247   isContinuous[1 ] = 1;
00248   isContinuous[2 ] = 1;
00249   isContinuous[3 ] = 1;
00250   isContinuous[4 ] = 0;
00251   isContinuous[5 ] = 0;
00252   isContinuous[6 ] = 1;
00253   isContinuous[7 ] = 1;
00254   isContinuous[8 ] = 1;
00255   isContinuous[9 ] = 1;
00256   isContinuous[10] = 1;
00257   isContinuous[11] = 1;
00258   isContinuous[12] = 1;
00259   isContinuous[13] = 1;
00260   isContinuous[14] = 1;
00261   isContinuous[15] = 1;
00262   isContinuous[16] = 1;
00263   isContinuous[17] = 1;
00264   isContinuous[18] = 1;
00265   isContinuous[19] = 1;
00266   isContinuous[20] = 1;
00267 
00268   //----Should the cumulative distribution be calculated for this variable?-----
00269   cdCompNeeded[0 ] = 1;
00270   cdCompNeeded[1 ] = 1;
00271   cdCompNeeded[2 ] = 1;
00272   cdCompNeeded[3 ] = 1;
00273   cdCompNeeded[4 ] = 1;
00274   cdCompNeeded[5 ] = 1;
00275   cdCompNeeded[6 ] = 1;
00276   cdCompNeeded[7 ] = 1;
00277   cdCompNeeded[8 ] = 1;
00278   cdCompNeeded[9 ] = 1;
00279   cdCompNeeded[10] = 0;
00280   cdCompNeeded[11] = 0;
00281   cdCompNeeded[12] = 0;
00282   cdCompNeeded[13] = 1;
00283   cdCompNeeded[14] = 1;
00284   cdCompNeeded[15] = 1;
00285   cdCompNeeded[16] = 1;
00286   cdCompNeeded[17] = 1;
00287   cdCompNeeded[18] = 1;
00288   cdCompNeeded[19] = 1;
00289   cdCompNeeded[20] = 1;
00290 }
00291 
00292 
00293 // ------------ method called for each event  ------------
00294 void MuIsoValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup){
00295   
00296   ++nEvents;
00297   edm::LogInfo("Tutorial") << "\nInvestigating event #" << nEvents<<"\n";
00298   
00299   // Get Muon Collection 
00300   edm::Handle<edm::View<reco::Muon> > muonsHandle; // 
00301   iEvent.getByLabel(Muon_Tag, muonsHandle);
00302   
00303   //Fill event entry in histogram of number of muons
00304   edm::LogInfo("Tutorial") << "Number of Muons: " << muonsHandle->size();
00305   theMuonData = muonsHandle->size();
00306   h_nMuons->Fill(theMuonData);
00307   
00308   //Fill historgams concerning muon isolation 
00309   uint iMuon=0;
00310   dbe->setCurrentFolder(dirName.c_str());
00311   for (MuonIterator muon = muonsHandle->begin(); muon != muonsHandle->end(); ++muon, ++iMuon ) {
00312     ++nIncMuons;
00313     if (requireCombinedMuon) {
00314       if (muon->combinedMuon().isNull()) continue;
00315     }
00316     //    ++nCombinedMuons;
00317     RecordData(muon);
00318     FillHistos();
00319   }
00320   dbe->cd();
00321   
00322 }
00323 
00324 //---------------Record data for a signle muon's data---------------------
00325 void MuIsoValidation::RecordData(MuonIterator muon){
00326   
00327   
00328   theData[0] = muon->isolationR03().sumPt;
00329   theData[1] = muon->isolationR03().emEt;
00330   theData[2] = muon->isolationR03().hadEt;
00331   theData[3] = muon->isolationR03().hoEt;
00332 
00333   theData[4] = muon->isolationR03().nTracks;
00334   theData[5] = muon->isolationR03().nJets;
00335   theData[6] = muon->isolationR03().trackerVetoPt;
00336   theData[7] = muon->isolationR03().emVetoEt;
00337   theData[8] = muon->isolationR03().hadVetoEt;
00338   theData[9] = muon->isolationR03().hoVetoEt;
00339   
00340   theData[10] = muon->pt();
00341   theData[11] = muon->eta();
00342   theData[12] = muon->phi();
00343 
00344   // make sure nTracks != 0 before filling this one
00345   if (theData[4] != 0) theData[13] = (double)theData[0] / (double)theData[4];
00346   else theData[13] = -99;
00347 
00348   theData[14] = 1.5 * theData[1] + theData[2];
00349 
00350   // Now PF isolation 
00351   theData[15] = -99.;
00352   theData[16] = -99.;
00353   theData[17] = -99.;
00354   theData[18] = -99.;
00355   theData[19] = -99.;
00356   theData[20] = -99.;
00357   if ( muon->isPFMuon() && muon->isPFIsolationValid() ) {
00358     theData[15] = muon->pfIsolationR03().sumChargedHadronPt;
00359     theData[16] = muon->pfIsolationR03().sumChargedParticlePt;
00360     theData[17] = muon->pfIsolationR03().sumChargedParticlePt-muon->pfIsolationR03().sumChargedHadronPt;
00361     theData[18] = muon->pfIsolationR03().sumNeutralHadronEt;
00362     theData[19] = muon->pfIsolationR03().sumPhotonEt;
00363     theData[20] = muon->pfIsolationR03().sumPUPt;
00364   }
00365 
00366 }
00367 
00368 // ------------ method called once each job just before starting event loop  ------------
00369 void 
00370 MuIsoValidation::beginJob()
00371 {
00372   
00373   edm::LogInfo("Tutorial") << "\n#########################################\n\n"
00374                            << "Lets get started! " 
00375                            << "\n\n#########################################\n";
00376   dbe->setCurrentFolder(dirName.c_str());
00377   InitHistos();
00378   dbe->cd();
00379   
00380 }
00381 
00382 // ------------ method called once each job just after ending the event loop  ------------
00383 void 
00384 MuIsoValidation::endJob() {
00385   
00386   // check if ME still there (and not killed by MEtoEDM for memory saving)
00387   if( dbe )
00388     {
00389       // check existence of first histo in the list
00390       if (! dbe->get(dirName+"/nMuons")) return;
00391     }
00392   else
00393     return;
00394 
00395   edm::LogInfo("Tutorial") << "\n#########################################\n\n"
00396                            << "Total Number of Events: " << nEvents
00397                            << "\nTotal Number of Muons: " << nIncMuons
00398                            << "\n\n#########################################\n"
00399                            << "\nInitializing Histograms...\n";
00400   
00401   edm::LogInfo("Tutorial") << "\nIntializing Finished.  Filling...\n";
00402   NormalizeHistos();
00403   edm::LogInfo("Tutorial") << "\nFilled.  Saving...\n";
00404   //  dbe->save(rootfilename); // comment out for incorporation
00405   edm::LogInfo("Tutorial") << "\nSaved.  Peace, homie, I'm out.\n";
00406   
00407 }
00408 
00409 void MuIsoValidation::InitHistos(){
00410   
00411   //---initialize number of muons histogram---
00412   h_nMuons = dbe->book1D("nMuons", title_sam + "Number of Muons", 20, 0., 20.);
00413   h_nMuons->setAxisTitle("Number of Muons",XAXIS);
00414   h_nMuons->setAxisTitle("Fraction of Events",YAXIS);
00415   
00416   
00417   //---Initialize 1D Histograms---
00418   for(int var = 0; var < NUM_VARS; var++){
00419     h_1D[var] = dbe->book1D(
00420                             names[var], 
00421                             title_sam + main_titles[var] + title_cone, 
00422                             (int)param[var][0], 
00423                             param[var][1], 
00424                             param[var][2]
00425                             );
00426     h_1D[var]->setAxisTitle(axis_titles[var],XAXIS);
00427     h_1D[var]->setAxisTitle("Fraction of Muons",YAXIS);
00428     GetTH1FromMonitorElement(h_1D[var])->Sumw2();
00429 
00430     if (cdCompNeeded[var]) {
00431       cd_plots[var] = dbe->book1D(
00432                                   names[var] + "_cd", 
00433                                   title_sam + title_cd + main_titles[var] + title_cone, 
00434                                   (int)param[var][0], 
00435                                   param[var][1], 
00436                                   param[var][2]
00437                                   );
00438       cd_plots[var]->setAxisTitle(axis_titles[var],XAXIS);
00439       cd_plots[var]->setAxisTitle("Fraction of Muons",YAXIS);
00440       GetTH1FromMonitorElement(cd_plots[var])->Sumw2();
00441     }
00442   }//Finish 1D
00443   
00444   //---Initialize 2D Histograms---
00445   for(int var1 = 0; var1 < NUM_VARS; var1++){
00446     for(int var2 = 0; var2 < NUM_VARS; var2++){
00447       if(var1 == var2) continue;
00448       
00449       /*      h_2D[var1][var2] = dbe->book2D(
00450               names[var1] + "_" + names[var2] + "_s",
00451               //title is in "y-var vs. x-var" format
00452               title_sam + main_titles[var2] + " <vs> " + main_titles[var1] + title_cone, 
00453               (int)param[var1][0],
00454               param[var1][1],
00455               param[var1][2],
00456               (int)param[var2][0],
00457               param[var2][1],
00458               param[var2][2]
00459               );
00460       */
00461       //Monitor elements is weird and takes y axis parameters as well
00462       //as x axis parameters for a 1D profile plot
00463       p_2D[var1][var2] = dbe->bookProfile(
00464                                           names[var1] + "_" + names[var2],
00465                                           title_sam + main_titles[var2] + " <vs> " + main_titles[var1] + title_cone,
00466                                           (int)param[var1][0],
00467                                           param[var1][1],
00468                                           param[var1][2],
00469                                           (int)param[var2][0], //documentation says this is disregarded
00470                                           param[var2][1],      //does this do anything?
00471                                           param[var2][2],      //does this do anything?
00472                                           " "                  //profile errors = spread/sqrt(num_datums)
00473                                           );
00474       
00475       if(LOG_BINNING_ENABLED && isContinuous[var1]){
00476         Double_t * bin_edges = new Double_t[NUM_LOG_BINS+1];
00477         // nbins+1 because there is one more edge than there are bins
00478         MakeLogBinsForProfile(bin_edges, param[var1][1], param[var1][2]);
00479         GetTProfileFromMonitorElement(p_2D[var1][var2])->SetBins(NUM_LOG_BINS, bin_edges);
00480         delete[] bin_edges;
00481       }
00482       /*      h_2D[var1][var2]->setAxisTitle(axis_titles[var1],XAXIS);
00483               h_2D[var1][var2]->setAxisTitle(axis_titles[var2],YAXIS);
00484               GetTH2FromMonitorElement(h_2D[var1][var2])->Sumw2();
00485       */
00486       
00487       p_2D[var1][var2]->setAxisTitle(axis_titles[var1],XAXIS);
00488       p_2D[var1][var2]->setAxisTitle(axis_titles[var2],YAXIS);
00489       //                        GetTProfileFromMonitorElement(p_2D[var1][var2])->Sumw2();
00490     }
00491   }//Finish 2D
00492   
00493   
00494   
00495   //avg pT not defined for zero tracks.
00496   //MonitorElement is inflxible and won't let me change the
00497   //number of bins!  I guess all I'm doing here is changing 
00498   //range of the x axis when it is printed, not the actual
00499   //bins that are filled
00500   p_2D[4][9]->setAxisRange(0.5,15.5,XAXIS);
00501   
00502 }
00503 
00504 void MuIsoValidation::MakeLogBinsForProfile(Double_t* bin_edges,        const double min, 
00505                                             const double max){
00506   
00507   const double &r = LOG_BINNING_RATIO;
00508   const int &nbins = NUM_LOG_BINS;
00509   
00510   const double first_bin_width = (r > 1.0) ? //so we don't divide by zero
00511     (max - min)*(1-r)/(1-pow(r,nbins)) :
00512     (max - min)/nbins;
00513   
00514   bin_edges[0] = min;
00515   bin_edges[1] = min + first_bin_width;
00516   for(int n = 2; n<nbins; ++n){
00517     bin_edges[n] = bin_edges[n-1] + (bin_edges[n-1] - bin_edges[n-2])*r;
00518   }
00519   bin_edges[nbins] = max;
00520 }
00521 
00522 void MuIsoValidation::NormalizeHistos() {
00523   for(int var=0; var<NUM_VARS; var++){   
00524     //turn cd_plots into CDF's
00525     //underflow -> bin #0.  overflow -> bin #(nbins+1)
00526     //0th bin doesn't need changed
00527     
00528     //----normalize------
00529     double entries = GetTH1FromMonitorElement(h_1D[var])->GetEntries();
00530     if (entries==0)continue;
00531     GetTH1FromMonitorElement(h_1D[var])->Scale(1./entries);
00532 
00533     if (cdCompNeeded[var]) {
00534       int n_max = int(param[var][0])+1;
00535       for(int n=1; n<=n_max; ++n){
00536         cd_plots[var]->setBinContent(n, cd_plots[var]->getBinContent(n) + cd_plots[var]->getBinContent(n-1)); //Integrate.
00537       }
00538       //----normalize------
00539       GetTH1FromMonitorElement(cd_plots[var])->Scale(1./entries);
00540     }
00541   }
00542 }
00543 
00544 void MuIsoValidation::FillHistos() {
00545   
00546   int overFlowBin;
00547   double overFlow = 0;
00548   
00549   //----------Fill 1D histograms---------------
00550   for(int var=0; var<NUM_VARS; var++){  
00551     h_1D[var]->Fill(theData[var]);
00552     if (cdCompNeeded[var]) cd_plots[var]->Fill(theData[var]);//right now, this is a regular PDF (just like h_1D)
00553     if (theData[var] > param[var][2]) {
00554       // fill the overflow bin
00555       overFlowBin = (int) param[var][0] + 1;
00556       overFlow = GetTH1FromMonitorElement(h_1D[var])->GetBinContent(overFlowBin);
00557       GetTH1FromMonitorElement(h_1D[var])->SetBinContent(overFlowBin, overFlow + 1);
00558     }
00559   }//Finish 1D
00560   
00561   //----------Fill 2D histograms---------------
00562   for(int var1=0; var1<NUM_VARS; ++var1){
00563     for(int var2=0; var2<NUM_VARS; ++var2){
00564       if(var1 == var2) continue;
00565       //change below to regular int interating!
00566       //      h_2D[var1][var2]->Fill(theData[var1], theData[var2]);
00567       p_2D[var1][var2]->Fill(theData[var1], theData[var2]);
00568     }
00569   }//Finish 2D
00570 }
00571 
00572 TH1* MuIsoValidation::GetTH1FromMonitorElement(MonitorElement* me) {
00573   return me->getTH1();
00574 }
00575 
00576 TH2* MuIsoValidation::GetTH2FromMonitorElement(MonitorElement* me) {
00577   return me->getTH2F();
00578 }
00579 
00580 TProfile* MuIsoValidation::GetTProfileFromMonitorElement(MonitorElement* me) {
00581   return me->getTProfile();
00582 }
00583 
00584 
00585 //define this as a plug-in
00586 DEFINE_FWK_MODULE(MuIsoValidation);
/data/refman/pasoursint/CMSSW_5_3_0/src/Validation/MuonIsolation/src/MuIsoValidation.cc
