00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "Validation/MuonIsolation/interface/MuIsoValidation.h"
00033
00034
00035 #include <memory>
00036 #include <string>
00037 #include <typeinfo>
00038 #include <utility>
00039
00040
00041 #include "TH1.h"
00042 #include "TH2.h"
00043 #include "TProfile.h"
00044
00045
00046 #include "FWCore/Framework/interface/MakerMacros.h"
00047 #include "FWCore/Framework/interface/Event.h"
00048 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00049
00050
00051 #include "DataFormats/TrackReco/interface/Track.h"
00052
00053
00054 using std::vector;
00055 using std::pair;
00056 using std::string;
00057
00058
00059
00060
00061
00062
00063 MuIsoValidation::MuIsoValidation(const edm::ParameterSet& iConfig)
00064 {
00065
00066
00067 requireCombinedMuon = iConfig.getUntrackedParameter<bool>("requireCombinedMuon");
00068 dirName = iConfig.getParameter<std::string>("directory");
00069
00070
00071
00072
00073
00074 Muon_Tag = iConfig.getUntrackedParameter<edm::InputTag>("Global_Muon_Label");
00075
00076
00077 nEvents = 0;
00078 nIncMuons = 0;
00079
00080
00081 InitStatics();
00082
00083
00084 dbe = 0;
00085 dbe = edm::Service<DQMStore>().operator->();
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099 h_1D.resize (NUM_VARS);
00100 cd_plots.resize(NUM_VARS);
00101
00102 p_2D.resize(NUM_VARS, vector<MonitorElement*>(NUM_VARS));
00103
00104 dbe->cd();
00105 }
00106
00107
00108
00109
00110 MuIsoValidation::~MuIsoValidation(){
00111
00112
00113
00114 }
00115
00116
00117
00118
00119 void MuIsoValidation::InitStatics(){
00120
00121
00122 S_BIN_WIDTH = 1.0;
00123 L_BIN_WIDTH = 2.0;
00124 LOG_BINNING_ENABLED = 1;
00125 NUM_LOG_BINS = 15;
00126 LOG_BINNING_RATIO = 1.1;
00127
00128
00129
00130
00131 title_sam = "";
00132 title_cone = "";
00133
00134
00135 title_cd = "C.D. of ";
00136
00137
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
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
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
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
00215
00216
00217
00218
00219
00220
00221
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
00245
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
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
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
00300 edm::Handle<edm::View<reco::Muon> > muonsHandle;
00301 iEvent.getByLabel(Muon_Tag, muonsHandle);
00302
00303
00304 edm::LogInfo("Tutorial") << "Number of Muons: " << muonsHandle->size();
00305 theMuonData = muonsHandle->size();
00306 h_nMuons->Fill(theMuonData);
00307
00308
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
00317 RecordData(muon);
00318 FillHistos();
00319 }
00320 dbe->cd();
00321
00322 }
00323
00324
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
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
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
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
00383 void
00384 MuIsoValidation::endJob() {
00385
00386
00387 if( dbe )
00388 {
00389
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
00405 edm::LogInfo("Tutorial") << "\nSaved. Peace, homie, I'm out.\n";
00406
00407 }
00408
00409 void MuIsoValidation::InitHistos(){
00410
00411
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
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 }
00443
00444
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
00450
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
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],
00470 param[var2][1],
00471 param[var2][2],
00472 " "
00473 );
00474
00475 if(LOG_BINNING_ENABLED && isContinuous[var1]){
00476 Double_t * bin_edges = new Double_t[NUM_LOG_BINS+1];
00477
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
00483
00484
00485
00486
00487 p_2D[var1][var2]->setAxisTitle(axis_titles[var1],XAXIS);
00488 p_2D[var1][var2]->setAxisTitle(axis_titles[var2],YAXIS);
00489
00490 }
00491 }
00492
00493
00494
00495
00496
00497
00498
00499
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) ?
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
00525
00526
00527
00528
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));
00537 }
00538
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
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]);
00553 if (theData[var] > param[var][2]) {
00554
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 }
00560
00561
00562 for(int var1=0; var1<NUM_VARS; ++var1){
00563 for(int var2=0; var2<NUM_VARS; ++var2){
00564 if(var1 == var2) continue;
00565
00566
00567 p_2D[var1][var2]->Fill(theData[var1], theData[var2]);
00568 }
00569 }
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
00586 DEFINE_FWK_MODULE(MuIsoValidation);