00001 #include <memory>
00002 #include <fstream>
00003 #include <iostream>
00004
00005 #include "FWCore/ServiceRegistry/interface/Service.h"
00006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00007 #include "FWCore/Framework/interface/MakerMacros.h"
00008 #include "FWCore/Framework/interface/Frameworkfwd.h"
00009 #include "FWCore/ParameterSet/interface/FileInPath.h"
00010 #include "DataFormats/DetId/interface/DetId.h"
00011 #include "DataFormats/EcalDetId/interface/ESDetId.h"
00012 #include "DataFormats/EcalDigi/interface/ESDataFrame.h"
00013 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00014 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
00015 #include "CondFormats/DataRecord/interface/ESGainRcd.h"
00016 #include "DQMServices/Core/interface/MonitorElement.h"
00017 #include "DQMServices/Core/interface/DQMStore.h"
00018 #include "DQM/EcalPreshowerMonitorModule/interface/ESTimingTask.h"
00019
00020 #include "TStyle.h"
00021 #include "TH2F.h"
00022 #include "TMath.h"
00023 #include "TGraph.h"
00024
00025 using namespace cms;
00026 using namespace edm;
00027 using namespace std;
00028
00029
00030 Double_t fitf(Double_t *x, Double_t *par) {
00031
00032 Double_t wc = par[2];
00033 Double_t n = par[3];
00034 Double_t v1 = pow(wc/n*(x[0]-par[1]), n);
00035 Double_t v2 = TMath::Exp(n-wc*(x[0]-par[1]));
00036 Double_t v = par[0]*v1*v2;
00037
00038 if (x[0] < par[1]) v = 0;
00039
00040 return v;
00041 }
00042
00043 ESTimingTask::ESTimingTask(const edm::ParameterSet& ps) {
00044
00045 rechitlabel_ = ps.getParameter<InputTag>("RecHitLabel");
00046 digilabel_ = ps.getParameter<InputTag>("DigiLabel");
00047 prefixME_ = ps.getUntrackedParameter<string>("prefixME", "EcalPreshower");
00048
00049 dqmStore_ = Service<DQMStore>().operator->();
00050 eCount_ = 0;
00051
00052
00053 for (int i = 0; i < 2; ++i)
00054 for (int j = 0; j < 2; ++j)
00055 hTiming_[i][j] = 0;
00056
00057 fit_ = new TF1("fit", fitf, -200, 200, 4);
00058 fit_->SetParameters(50, 10, 0, 0);
00059
00060 dqmStore_->setCurrentFolder(prefixME_ + "/ESTimingTask");
00061
00062
00063
00064 char histo[200];
00065 for (int i=0 ; i<2; ++i)
00066 for (int j=0 ; j<2; ++j) {
00067 int iz = (i==0)? 1:-1;
00068 sprintf(histo, "ES Timing Z %d P %d", iz, j+1);
00069 hTiming_[i][j] = dqmStore_->book1D(histo, histo, 81, -20.5, 20.5);
00070 hTiming_[i][j]->setAxisTitle("ES Timing (ns)", 1);
00071 }
00072
00073 sprintf(histo, "ES 2D Timing");
00074 h2DTiming_ = dqmStore_->book2D(histo, histo, 81, -20.5, 20.5, 81, -20.5, 20.5);
00075 h2DTiming_->setAxisTitle("ES- Timing (ns)", 1);
00076 h2DTiming_->setAxisTitle("ES+ Timing (ns)", 2);
00077
00078 htESP_ = new TH1F("htESP", "Timing ES+", 81, -20.5, 20.5);
00079 htESM_ = new TH1F("htESM", "Timing ES-", 81, -20.5, 20.5);
00080 }
00081
00082 ESTimingTask::~ESTimingTask() {
00083 delete htESP_;
00084 delete htESM_;
00085 }
00086
00087 void ESTimingTask::beginJob(void) {
00088 }
00089
00090 void ESTimingTask::endJob() {
00091 }
00092
00093 void ESTimingTask::analyze(const edm::Event& e, const edm::EventSetup& iSetup) {
00094
00095 set(iSetup);
00096
00097 runNum_ = e.id().run();
00098 eCount_++;
00099
00100 htESP_->Reset();
00101 htESM_->Reset();
00102
00103
00104 int zside, plane, ix, iy, is;
00105 double adc[3], para[10];
00106 double tx[3] = {-5., 20., 45.};
00107 Handle<ESDigiCollection> digis;
00108 if ( e.getByLabel(digilabel_, digis) ) {
00109
00110 for (ESDigiCollection::const_iterator digiItr = digis->begin(); digiItr != digis->end(); ++digiItr) {
00111
00112 ESDataFrame dataframe = (*digiItr);
00113 ESDetId id = dataframe.id();
00114
00115 zside = id.zside();
00116 plane = id.plane();
00117 ix = id.six();
00118 iy = id.siy();
00119 is = id.strip();
00120
00121
00122 if (zside==1 && plane==1 && ix==7 && iy==28) continue;
00123 if (zside==1 && plane==1 && ix==24 && iy==9 && is==21) continue;
00124 if (zside==-1 && plane==2 && ix==35 && iy==17 && is==23) continue;
00125
00126 int i = (zside==1)? 0:1;
00127 int j = plane-1;
00128
00129 for (int k=0; k<dataframe.size(); ++k)
00130 adc[k] = dataframe.sample(k).adc();
00131
00132 double status = 0;
00133 if (adc[1] < 200) status = 1;
00134 if (fabs(adc[0]) > 10) status = 1;
00135 if (adc[1] < 0 || adc[2] < 0) status = 1;
00136 if (adc[0] > adc[1] || adc[0] > adc[2]) status = 1;
00137 if (adc[2] > adc[1]) status = 1;
00138
00139 if (int(status) == 0) {
00140 TGraph *gr = new TGraph(3, tx, adc);
00141 fit_->SetParameters(50, 10, wc_, n_);
00142 fit_->FixParameter(2, wc_);
00143 fit_->FixParameter(3, n_);
00144 gr->Fit("fit", "MQ");
00145 fit_->GetParameters(para);
00146 delete gr;
00147
00148 hTiming_[i][j]->Fill(para[1]);
00149
00150 if (zside == 1) htESP_->Fill(para[1]);
00151 else if (zside == -1) htESM_->Fill(para[1]);
00152 }
00153
00154 }
00155 } else {
00156 LogWarning("ESTimingTask") << digilabel_ << " not available";
00157 }
00158
00159 if (htESP_->GetEntries() > 0 && htESM_->GetEntries() > 0)
00160 h2DTiming_->Fill(htESM_->GetMean(), htESP_->GetMean());
00161
00162 }
00163
00164 void ESTimingTask::set(const edm::EventSetup& es) {
00165
00166 es.get<ESGainRcd>().get(esgain_);
00167 const ESGain *gain = esgain_.product();
00168
00169 double ESGain = gain->getESGain();
00170
00171 if (ESGain == 1) {
00172 wc_ = 0.0837264;
00173 n_ = 2.016;
00174 } else {
00175 wc_ = 0.07291;
00176 n_ = 1.798;
00177 }
00178
00179 }
00180
00181
00182 DEFINE_FWK_MODULE(ESTimingTask);