EcalPreshowerMonitorModule/src/ESOccupancyTask.cc

0001 #include <memory>
0002 #include <fstream>
0003 #include <iostream>
0004
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Framework/interface/MakerMacros.h"
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/ParameterSet/interface/FileInPath.h"
0010 #include "DataFormats/DetId/interface/DetId.h"
0011 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0012 #include "DataFormats/EcalDigi/interface/ESDataFrame.h"
0013 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0014 #include "DQMServices/Core/interface/DQMStore.h"
0015 #include "DQM/EcalPreshowerMonitorModule/interface/ESOccupancyTask.h"
0016
0017 #include "TStyle.h"
0018 #include "TH2F.h"
0019
0020 using namespace cms;
0021 using namespace edm;
0022 using namespace std;
0023
0024 ESOccupancyTask::ESOccupancyTask(const edm::ParameterSet& ps) {
0025   rechittoken_ = consumes<ESRecHitCollection>(ps.getParameter<InputTag>("RecHitLabel"));
0026   prefixME_ = ps.getUntrackedParameter<string>("prefixME", "EcalPreshower");
0027
0028   eCount_ = 0;
0029
0030   //Histogram init
0031   for (int i = 0; i < 2; ++i)
0032     for (int j = 0; j < 2; ++j) {
0033       hRecOCC_[i][j] = nullptr;
0034       hRecNHit_[i][j] = nullptr;
0035       hEng_[i][j] = nullptr;
0036       hEvEng_[i][j] = nullptr;
0037       hEnDensity_[i][i] = nullptr;
0038       hGoodRecNHit_[i][j] = nullptr;
0039
0040       hSelEng_[i][j] = nullptr;
0041       hSelOCC_[i][j] = nullptr;
0042       hSelEnDensity_[i][j] = nullptr;
0043     }
0044
0045   for (int i = 0; i < 2; ++i)
0046     hE1E2_[i] = nullptr;
0047 }
0048
0049 void ESOccupancyTask::bookHistograms(DQMStore::IBooker& iBooker, Run const&, EventSetup const&) {
0050   iBooker.setCurrentFolder(prefixME_ + "/ESOccupancyTask");
0051
0052   //Booking Histograms
0053   //Notice: Change ESRenderPlugin under DQM/RenderPlugins/src if you change this histogram name.
0054   char histo[200];
0055   for (int i = 0; i < 2; ++i)
0056     for (int j = 0; j < 2; ++j) {
0057       int iz = (i == 0) ? 1 : -1;
0058       sprintf(histo, "ES RecHit 2D Occupancy Z %d P %d", iz, j + 1);
0059       hRecOCC_[i][j] = iBooker.book2D(histo, histo, 40, 0.5, 40.5, 40, 0.5, 40.5);
0060       hRecOCC_[i][j]->setAxisTitle("Si X", 1);
0061       hRecOCC_[i][j]->setAxisTitle("Si Y", 2);
0062
0063       //Bin 40,40 is used to save eumber of event for scaling.
0064       sprintf(histo, "ES Energy Density Z %d P %d", iz, j + 1);
0065       hEnDensity_[i][j] = iBooker.book2D(histo, histo, 40, 0.5, 40.5, 40, 0.5, 40.5);
0066       hEnDensity_[i][j]->setAxisTitle("Si X", 1);
0067       hEnDensity_[i][j]->setAxisTitle("Si Y", 2);
0068
0069       sprintf(histo, "ES Num of RecHits Z %d P %d", iz, j + 1);
0070       hRecNHit_[i][j] = iBooker.book1DD(histo, histo, 60, 0, 1920);
0071       hRecNHit_[i][j]->setAxisTitle("# of RecHits", 1);
0072       hRecNHit_[i][j]->setAxisTitle("Num of Events", 2);
0073
0074       sprintf(histo, "ES Num of Good RecHits Z %d P %d", iz, j + 1);
0075       hGoodRecNHit_[i][j] = iBooker.book1DD(histo, histo, 60, 0, 1920);
0076       hGoodRecNHit_[i][j]->setAxisTitle("# of good RecHits", 1);
0077       hGoodRecNHit_[i][j]->setAxisTitle("Num of Events", 2);
0078
0079       sprintf(histo, "ES RecHit Energy Z %d P %d", iz, j + 1);
0080       hEng_[i][j] = iBooker.book1DD(histo, histo, 50, 0, 0.001);
0081       hEng_[i][j]->setAxisTitle("RecHit Energy", 1);
0082       hEng_[i][j]->setAxisTitle("Num of ReHits", 2);
0083
0084       sprintf(histo, "ES Event Energy Z %d P %d", iz, j + 1);
0085       hEvEng_[i][j] = iBooker.book1DD(histo, histo, 50, 0, 0.1);
0086       hEvEng_[i][j]->setAxisTitle("Event Energy", 1);
0087       hEvEng_[i][j]->setAxisTitle("Num of Events", 2);
0088
0089       // histograms with selected hits
0090       sprintf(histo, "ES RecHit Energy with selected hits Z %d P %d", iz, j + 1);
0091       hSelEng_[i][j] = iBooker.book1DD(histo, histo, 50, 0, 0.001);
0092       hSelEng_[i][j]->setAxisTitle("RecHit Energy", 1);
0093       hSelEng_[i][j]->setAxisTitle("Num of ReHits", 2);
0094
0095       sprintf(histo, "ES Occupancy with selected hits Z %d P %d", iz, j + 1);
0096       hSelOCC_[i][j] = iBooker.book2D(histo, histo, 40, 0.5, 40.5, 40, 0.5, 40.5);
0097       hSelOCC_[i][j]->setAxisTitle("Si X", 1);
0098       hSelOCC_[i][j]->setAxisTitle("Si Y", 2);
0099
0100       sprintf(histo, "ES Energy Density with selected hits Z %d P %d", iz, j + 1);
0101       hSelEnDensity_[i][j] = iBooker.book2D(histo, histo, 40, 0.5, 40.5, 40, 0.5, 40.5);
0102       hSelEnDensity_[i][j]->setAxisTitle("Si X", 1);
0103       hSelEnDensity_[i][j]->setAxisTitle("Si Y", 2);
0104     }
0105
0106   hE1E2_[0] = iBooker.book2D("ES+ EP1 vs EP2", "ES+ EP1 vs EP2", 50, 0, 0.1, 50, 0, 0.1);
0107   hE1E2_[1] = iBooker.book2D("ES- EP1 vs EP2", "ES- EP1 vs EP2", 50, 0, 0.1, 50, 0, 0.1);
0108 }
0109
0110 void ESOccupancyTask::analyze(const edm::Event& e, const edm::EventSetup& iSetup) {
0111   runNum_ = e.id().run();
0112   eCount_++;
0113
0114   // RecHits
0115   int zside, plane, ix, iy;
0116   int sum_RecHits[2][2], sum_GoodRecHits[2][2];
0117   float sum_Energy[2][2];
0118
0119   for (int i = 0; i < 2; ++i)
0120     for (int j = 0; j < 2; ++j) {
0121       sum_RecHits[i][j] = 0;
0122       sum_GoodRecHits[i][j] = 0;
0123       sum_Energy[i][j] = 0;
0124     }
0125
0126   Handle<ESRecHitCollection> ESRecHit;
0127   if (e.getByToken(rechittoken_, ESRecHit)) {
0128     for (ESRecHitCollection::const_iterator hitItr = ESRecHit->begin(); hitItr != ESRecHit->end(); ++hitItr) {
0129       ESDetId id = ESDetId(hitItr->id());
0130
0131       zside = id.zside();
0132       plane = id.plane();
0133       ix = id.six();
0134       iy = id.siy();
0135
0136       int i = (zside == 1) ? 0 : 1;
0137       int j = plane - 1;
0138
0139       sum_RecHits[i][j]++;
0140       sum_Energy[i][j] += hitItr->energy();
0141       hRecOCC_[i][j]->Fill(ix, iy);
0142       if (hitItr->energy() != 0) {
0143         hEng_[i][j]->Fill(hitItr->energy());
0144         hEnDensity_[i][j]->Fill(ix, iy, hitItr->energy());
0145
0146         if (hitItr->recoFlag() == 14 || hitItr->recoFlag() == 1 ||
0147             (hitItr->recoFlag() <= 10 && hitItr->recoFlag() >= 5))
0148           continue;
0149         sum_GoodRecHits[i][j]++;
0150         hSelEng_[i][j]->Fill(hitItr->energy());
0151         hSelEnDensity_[i][j]->Fill(ix, iy, hitItr->energy());
0152         hSelOCC_[i][j]->Fill(ix, iy);
0153       }
0154     }
0155   } else {
0156     LogWarning("ESOccupancyTask") << "RecHitCollection not available";
0157   }
0158
0159   //Fill histograms after a event
0160   for (int i = 0; i < 2; ++i)
0161     for (int j = 0; j < 2; ++j) {
0162       hRecNHit_[i][j]->Fill(sum_RecHits[i][j]);
0163       hGoodRecNHit_[i][j]->Fill(sum_GoodRecHits[i][j]);
0164       hEvEng_[i][j]->Fill(sum_Energy[i][j]);
0165
0166       //Save eCount_ for Scaling
0167       hRecOCC_[i][j]->setBinContent(40, 40, eCount_);
0168       hEnDensity_[i][j]->setBinContent(40, 40, eCount_);
0169
0170       hSelOCC_[i][j]->setBinContent(40, 40, eCount_);
0171       hSelEnDensity_[i][j]->setBinContent(40, 40, eCount_);
0172     }
0173
0174   hE1E2_[0]->Fill(sum_Energy[0][0], sum_Energy[0][1]);
0175   hE1E2_[1]->Fill(sum_Energy[1][0], sum_Energy[1][1]);
0176 }
0177
0178 //define this as a plug-in
0179 DEFINE_FWK_MODULE(ESOccupancyTask);