CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch13/src/DQM/EcalPreshowerMonitorModule/src/ESTimingTask.cc

Go to the documentation of this file.
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 // fit function
00030 Double_t fitf(Double_t *x, Double_t *par) {
00031 
00032   Double_t wc = par[2];
00033   Double_t n  = par[3]; // n-1 (in fact)
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   //Histogram init  
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   //Booking Histograms
00063   //Notice: Change ESRenderPlugin under DQM/RenderPlugins/src if you change this histogram name.
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    //Digis
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        //if (zside==1 && plane==1 && ix==15 && iy==6) continue;       
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          //cout<<"ADC : "<<zside<<" "<<plane<<" "<<ix<<" "<<iy<<" "<<is<<" "<<adc[0]<<" "<<adc[1]<<" "<<adc[2]<<" "<<para[1]<<endl;
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) { // LG
00172     wc_ = 0.0837264;
00173     n_  = 2.016;
00174   } else { // HG
00175     wc_ = 0.07291;
00176     n_  = 1.798; 
00177   }
00178 
00179 }
00180 
00181 //define this as a plug-in
00182 DEFINE_FWK_MODULE(ESTimingTask);