EcalPreshowerMonitorModule/src/ESTimingTask.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/EcalRecHit/interface/EcalRecHitCollections.h"
0014 #include "CondFormats/DataRecord/interface/ESGainRcd.h"
0015 #include "DQMServices/Core/interface/DQMStore.h"
0016 #include "DQM/EcalPreshowerMonitorModule/interface/ESTimingTask.h"
0017
0018 #include "TMath.h"
0019 #include "TGraph.h"
0020
0021 using namespace cms;
0022 using namespace edm;
0023 using namespace std;
0024
0025 // fit function
0026 double fitf(double* x, double* par) {
0027   double wc = par[2];
0028   double n = par[3];  // n-1 (in fact)
0029   double v1 = pow(wc / n * (x[0] - par[1]), n);
0030   double v2 = TMath::Exp(n - wc * (x[0] - par[1]));
0031   double v = par[0] * v1 * v2;
0032
0033   if (x[0] < par[1])
0034     v = 0;
0035
0036   return v;
0037 }
0038
0039 ESTimingTask::ESTimingTask(const edm::ParameterSet& ps) {
0040   digilabel_ = consumes<ESDigiCollection>(ps.getParameter<InputTag>("DigiLabel"));
0041   prefixME_ = ps.getUntrackedParameter<string>("prefixME", "EcalPreshower");
0042   esgainToken_ = esConsumes();
0043   eCount_ = 0;
0044
0045   fit_ = new TF1("fitShape", fitf, -200, 200, 4);
0046   fit_->SetParameters(50, 10, 0, 0);
0047
0048   //Histogram init
0049   for (int i = 0; i < 2; ++i)
0050     for (int j = 0; j < 2; ++j)
0051       hTiming_[i][j] = nullptr;
0052
0053   htESP_ = new TH1F("htESP", "Timing ES+", 81, -20.5, 20.5);
0054   htESM_ = new TH1F("htESM", "Timing ES-", 81, -20.5, 20.5);
0055 }
0056
0057 void ESTimingTask::bookHistograms(DQMStore::IBooker& iBooker, Run const&, EventSetup const&) {
0058   iBooker.setCurrentFolder(prefixME_ + "/ESTimingTask");
0059
0060   //Booking Histograms
0061   //Notice: Change ESRenderPlugin under DQM/RenderPlugins/src if you change this histogram name.
0062   char histo[200];
0063   for (int i = 0; i < 2; ++i)
0064     for (int j = 0; j < 2; ++j) {
0065       int iz = (i == 0) ? 1 : -1;
0066       sprintf(histo, "ES Timing Z %d P %d", iz, j + 1);
0067       hTiming_[i][j] = iBooker.book1D(histo, histo, 81, -20.5, 20.5);
0068       hTiming_[i][j]->setAxisTitle("ES Timing (ns)", 1);
0069     }
0070
0071   sprintf(histo, "ES 2D Timing");
0072   h2DTiming_ = iBooker.book2D(histo, histo, 81, -20.5, 20.5, 81, -20.5, 20.5);
0073   h2DTiming_->setAxisTitle("ES- Timing (ns)", 1);
0074   h2DTiming_->setAxisTitle("ES+ Timing (ns)", 2);
0075 }
0076
0077 ESTimingTask::~ESTimingTask() {
0078   delete fit_;
0079   delete htESP_;
0080   delete htESM_;
0081 }
0082
0083 void ESTimingTask::analyze(const edm::Event& e, const edm::EventSetup& iSetup) {
0084   set(iSetup);
0085
0086   runNum_ = e.id().run();
0087   eCount_++;
0088
0089   htESP_->Reset();
0090   htESM_->Reset();
0091
0092   //Digis
0093   int zside, plane, ix, iy, is;
0094   double adc[3];
0095   //  double para[10];
0096   //double tx[3] = {-5., 20., 45.};
0097   Handle<ESDigiCollection> digis;
0098   if (e.getByToken(digilabel_, digis)) {
0099     for (ESDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
0100       ESDataFrame dataframe = (*digiItr);
0101       ESDetId id = dataframe.id();
0102
0103       zside = id.zside();
0104       plane = id.plane();
0105       ix = id.six();
0106       iy = id.siy();
0107       is = id.strip();
0108
0109       //if (zside==1 && plane==1 && ix==15 && iy==6) continue;
0110       if (zside == 1 && plane == 1 && ix == 7 && iy == 28)
0111         continue;
0112       if (zside == 1 && plane == 1 && ix == 24 && iy == 9 && is == 21)
0113         continue;
0114       if (zside == -1 && plane == 2 && ix == 35 && iy == 17 && is == 23)
0115         continue;
0116
0117       int i = (zside == 1) ? 0 : 1;
0118       int j = plane - 1;
0119
0120       for (int k = 0; k < dataframe.size(); ++k)
0121         adc[k] = dataframe.sample(k).adc();
0122
0123       double status = 0;
0124       if (adc[1] < 200)
0125         status = 1;
0126       if (fabs(adc[0]) > 10)
0127         status = 1;
0128       if (adc[1] < 0 || adc[2] < 0)
0129         status = 1;
0130       if (adc[0] > adc[1] || adc[0] > adc[2])
0131         status = 1;
0132       if (adc[2] > adc[1])
0133         status = 1;
0134
0135       if (int(status) == 0) {
0136         double A1 = adc[1];
0137         double A2 = adc[2];
0138         double DeltaT = 25.;
0139         double aaa = (A2 > 0 && A1 > 0) ? log(A2 / A1) / n_ : 20.;  // if A1=0, t0=20
0140         double bbb = wc_ / n_ * DeltaT;
0141         double ccc = exp(aaa + bbb);
0142
0143         double t0 = (2. - ccc) / (1. - ccc) * DeltaT - 5;
0144         hTiming_[i][j]->Fill(t0);
0145         //cout<<"t0 : "<<t0<<endl;
0146         /*
0147     TGraph *gr = new TGraph(3, tx, adc);
0148     fit_->SetParameters(50, 10, wc_, n_);
0149     fit_->FixParameter(2, wc_);
0150     fit_->FixParameter(3, n_);
0151     fit_->Print();
0152     gr->Fit("fitShape", "MQ");
0153     fit_->GetParameters(para);
0154     delete gr;
0155     //hTiming_[i][j]->Fill(para[1]);
0156     */
0157         //cout<<"ADC : "<<zside<<" "<<plane<<" "<<ix<<" "<<iy<<" "<<is<<" "<<adc[0]<<" "<<adc[1]<<" "<<adc[2]<<" "<<para[1]<<" "<<wc_<<" "<<n_<<endl;
0158
0159         if (zside == 1)
0160           htESP_->Fill(t0);
0161         else if (zside == -1)
0162           htESM_->Fill(t0);
0163       }
0164     }
0165   } else {
0166     LogWarning("ESTimingTask") << "DigiCollection not available";
0167   }
0168
0169   if (htESP_->GetEntries() > 0 && htESM_->GetEntries() > 0)
0170     h2DTiming_->Fill(htESM_->GetMean(), htESP_->GetMean());
0171 }
0172
0173 void ESTimingTask::set(const edm::EventSetup& es) {
0174   const ESGain* gain = &es.getData(esgainToken_);
0175
0176   int ESGain = (int)gain->getESGain();
0177
0178   if (ESGain == 1) {  // LG
0179     wc_ = 0.0837264;
0180     n_ = 2.016;
0181   } else {  // HG
0182     wc_ = 0.07291;
0183     n_ = 1.798;
0184   }
0185
0186   //cout<<"gain : "<<ESGain<<endl;
0187   //cout<<wc_<<" "<<n_<<endl;
0188 }
0189
0190 //define this as a plug-in
0191 DEFINE_FWK_MODULE(ESTimingTask);