CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_0/src/DQM/EcalBarrelMonitorTasks/src/EBTimingTask.cc

Go to the documentation of this file.
00001 /*
00002  * \file EBTimingTask.cc
00003  *
00004  * $Date: 2011/09/15 21:54:51 $
00005  * $Revision: 1.73 $
00006  * \author G. Della Ricca
00007  *
00008 */
00009 
00010 #include <iostream>
00011 #include <fstream>
00012 #include <string>
00013 #include <cmath>
00014 
00015 #include "FWCore/ServiceRegistry/interface/Service.h"
00016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00017 #include "FWCore/Framework/interface/ESHandle.h"
00018 
00019 #include "DQMServices/Core/interface/MonitorElement.h"
00020 
00021 #include "DQMServices/Core/interface/DQMStore.h"
00022 
00023 #include "DataFormats/EcalRawData/interface/EcalRawDataCollections.h"
00024 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00025 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
00026 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00027 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerEvmReadoutRecord.h"
00028 
00029 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00030 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00031 
00032 #include "DQM/EcalCommon/interface/Numbers.h"
00033 
00034 #include "DQM/EcalBarrelMonitorTasks/interface/EBTimingTask.h"
00035 
00036 EBTimingTask::EBTimingTask(const edm::ParameterSet& ps){
00037 
00038   init_ = false;
00039 
00040   dqmStore_ = edm::Service<DQMStore>().operator->();
00041 
00042   prefixME_ = ps.getUntrackedParameter<std::string>("prefixME", "");
00043 
00044   enableCleanup_ = ps.getUntrackedParameter<bool>("enableCleanup", false);
00045 
00046   mergeRuns_ = ps.getUntrackedParameter<bool>("mergeRuns", false);
00047 
00048   energyThreshold_ = ps.getUntrackedParameter<double>("energyTreshold",1.0);
00049 
00050   EcalRawDataCollection_ = ps.getParameter<edm::InputTag>("EcalRawDataCollection");
00051   EcalRecHitCollection_ = ps.getParameter<edm::InputTag>("EcalRecHitCollection");
00052 
00053   L1GtEvmReadoutRecord_ = ps.getParameter<edm::InputTag>("L1GtEvmReadoutRecord");
00054 
00055   useBeamStatus_ = ps.getUntrackedParameter<bool>("useBeamStatus", false);
00056 
00057   for (int i = 0; i < 36; i++) {
00058     meTime_[i] = 0;
00059     meTimeMap_[i] = 0;
00060     meTimeAmpli_[i] = 0;
00061   }
00062 
00063   meTimeAmpliSummary_ = 0;
00064   meTimeSummary1D_ = 0;
00065   meTimeSummaryMap_ = 0;
00066 
00067   stableBeamsDeclared_ = false;
00068 
00069 }
00070 
00071 EBTimingTask::~EBTimingTask(){
00072 
00073 }
00074 
00075 void EBTimingTask::beginJob(void){
00076 
00077   ievt_ = 0;
00078 
00079   if ( dqmStore_ ) {
00080     dqmStore_->setCurrentFolder(prefixME_ + "/EBTimingTask");
00081     dqmStore_->rmdir(prefixME_ + "/EBTimingTask");
00082   }
00083 
00084 }
00085 
00086 void EBTimingTask::beginRun(const edm::Run& r, const edm::EventSetup& c) {
00087 
00088   Numbers::initGeometry(c, false);
00089 
00090   if ( ! mergeRuns_ ) this->reset();
00091 
00092   stableBeamsDeclared_ = false;
00093 
00094 }
00095 
00096 void EBTimingTask::endRun(const edm::Run& r, const edm::EventSetup& c) {
00097 
00098 }
00099 
00100 void EBTimingTask::reset(void) {
00101 
00102   for (int i = 0; i < 36; i++) {
00103     if ( meTime_[i] ) meTime_[i]->Reset();
00104     if ( meTimeMap_[i] ) meTimeMap_[i]->Reset();
00105     if ( meTimeAmpli_[i] ) meTimeAmpli_[i]->Reset();
00106   }
00107 
00108   if ( meTimeAmpliSummary_ ) meTimeAmpliSummary_->Reset();
00109   if ( meTimeSummary1D_ ) meTimeSummary1D_->Reset();
00110   if ( meTimeSummaryMap_ ) meTimeSummaryMap_->Reset();
00111 
00112 }
00113 
00114 void EBTimingTask::setup(void){
00115 
00116   init_ = true;
00117 
00118   std::string name;
00119 
00120   // for timing vs amplitude plots
00121   const int nbinsE = 25;
00122   const float minlogE = -0.5;
00123   const float maxlogE = 2.;
00124   float binEdgesE[nbinsE + 1];
00125   for(int i = 0; i <= nbinsE; i++)
00126     binEdgesE[i] = std::pow((float)10., minlogE + (maxlogE - minlogE) / nbinsE * i);
00127 
00128   const int nbinsT = 200;
00129   const float minT = -50.;
00130   const float maxT = 50.;
00131   float binEdgesT[nbinsT + 1];
00132   for(int i = 0; i <= nbinsT; i++)
00133     binEdgesT[i] = minT + (maxT - minT) / nbinsT * i;
00134 
00135   if ( dqmStore_ ) {
00136     dqmStore_->setCurrentFolder(prefixME_ + "/EBTimingTask");
00137 
00138     for (int i = 0; i < 36; i++) {
00139       name = "EBTMT timing 1D " + Numbers::sEB(i+1);
00140       meTime_[i] = dqmStore_->book1D(name, name, 50, -25., 25.);
00141       meTime_[i]->setAxisTitle("time (ns)", 1);
00142       dqmStore_->tag(meTime_[i], i+1);
00143 
00144       name = "EBTMT timing " + Numbers::sEB(i+1);
00145       meTimeMap_[i] = dqmStore_->bookProfile2D(name, name, 85, 0., 85., 20, 0., 20., -20.+shiftProf2D_, 20.+shiftProf2D_, "s");
00146       meTimeMap_[i]->setAxisTitle("ieta", 1);
00147       meTimeMap_[i]->setAxisTitle("iphi", 2);
00148       meTimeMap_[i]->setAxisTitle("time (ns)", 3);
00149       dqmStore_->tag(meTimeMap_[i], i+1);
00150 
00151       name = "EBTMT timing vs amplitude " + Numbers::sEB(i+1);
00152       meTimeAmpli_[i] = dqmStore_->book2D(name, name, nbinsE, binEdgesE, nbinsT, binEdgesT);
00153       meTimeAmpli_[i]->setAxisTitle("energy (GeV)", 1);
00154       meTimeAmpli_[i]->setAxisTitle("time (ns)", 2);
00155       dqmStore_->tag(meTimeAmpli_[i], i+1);
00156     }
00157 
00158     name = "EBTMT timing vs amplitude summary";
00159     meTimeAmpliSummary_ = dqmStore_->book2D(name, name, nbinsE, binEdgesE, nbinsT, binEdgesT);
00160     meTimeAmpliSummary_->setAxisTitle("energy (GeV)", 1);
00161     meTimeAmpliSummary_->setAxisTitle("time (ns)", 2);
00162 
00163     name = "EBTMT timing 1D summary";
00164     meTimeSummary1D_ = dqmStore_->book1D(name, name, 50, -25., 25.);
00165     meTimeSummary1D_->setAxisTitle("time (ns)", 1);
00166 
00167     name = "EBTMT timing map";
00168     meTimeSummaryMap_ = dqmStore_->bookProfile2D(name, name, 72, 0., 360., 34, -85, 85, -20.+shiftProf2D_, 20.+shiftProf2D_, "s");
00169     meTimeSummaryMap_->setAxisTitle("jphi", 1);
00170     meTimeSummaryMap_->setAxisTitle("jeta", 2);
00171     meTimeSummaryMap_->setAxisTitle("time (ns)", 3);
00172 
00173   }
00174 
00175 }
00176 
00177 void EBTimingTask::cleanup(void){
00178 
00179   if ( ! init_ ) return;
00180 
00181   if ( dqmStore_ ) {
00182     dqmStore_->setCurrentFolder(prefixME_ + "/EBTimingTask");
00183 
00184     for ( int i = 0; i < 36; i++ ) {
00185       if ( meTime_[i] ) dqmStore_->removeElement( meTime_[i]->getName() );
00186       meTime_[i] = 0;
00187 
00188       if ( meTimeMap_[i] ) dqmStore_->removeElement( meTimeMap_[i]->getName() );
00189       meTimeMap_[i] = 0;
00190 
00191       if ( meTimeAmpli_[i] ) dqmStore_->removeElement( meTimeAmpli_[i]->getName() );
00192       meTimeAmpli_[i] = 0;
00193     }
00194 
00195     if ( meTimeAmpliSummary_ ) dqmStore_->removeElement( meTimeAmpliSummary_->getName() );
00196     meTimeAmpliSummary_ = 0;
00197 
00198     if ( meTimeSummary1D_ ) dqmStore_->removeElement( meTimeSummary1D_->getName() );
00199     meTimeSummary1D_ = 0;
00200 
00201     if ( meTimeSummaryMap_ ) dqmStore_->removeElement( meTimeSummaryMap_->getName() );
00202     meTimeSummaryMap_ = 0;
00203 
00204   }
00205 
00206   init_ = false;
00207 
00208 }
00209 
00210 void EBTimingTask::endJob(void){
00211 
00212   edm::LogInfo("EBTimingTask") << "analyzed " << ievt_ << " events";
00213 
00214   if ( enableCleanup_ ) this->cleanup();
00215 
00216 }
00217 
00218 void EBTimingTask::analyze(const edm::Event& e, const edm::EventSetup& c){
00219 
00220   const unsigned STABLE_BEAMS = 11;
00221 
00222   bool isData = true;
00223   bool enable = false;
00224   int runType[36];
00225   for (int i=0; i<36; i++) runType[i] = -1;
00226 
00227   edm::Handle<EcalRawDataCollection> dcchs;
00228 
00229   if ( e.getByLabel(EcalRawDataCollection_, dcchs) ) {
00230 
00231     for ( EcalRawDataCollection::const_iterator dcchItr = dcchs->begin(); dcchItr != dcchs->end(); ++dcchItr ) {
00232 
00233       if ( Numbers::subDet( *dcchItr ) != EcalBarrel ) continue;
00234 
00235       int ism = Numbers::iSM( *dcchItr, EcalBarrel );
00236 
00237       runType[ism-1] = dcchItr->getRunType();
00238 
00239       if ( dcchItr->getRunType() == EcalDCCHeaderBlock::COSMIC ||
00240            dcchItr->getRunType() == EcalDCCHeaderBlock::MTCC ||
00241            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00242            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00243            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00244            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_LOCAL ) enable = true;
00245 
00246     }
00247 
00248   } else {
00249 
00250     isData = false; enable = true;
00251     edm::LogWarning("EBTimingTask") << EcalRawDataCollection_ << " not available";
00252 
00253   }
00254 
00255   if ( ! enable ) return;
00256 
00257   if ( ! init_ ) this->setup();
00258 
00259   ievt_++;
00260 
00261   // resetting plots when stable beam is declared
00262   if( useBeamStatus_ && !stableBeamsDeclared_ ) {
00263     edm::Handle<L1GlobalTriggerEvmReadoutRecord> gtRecord;
00264     if( e.getByLabel(L1GtEvmReadoutRecord_, gtRecord) ) {
00265 
00266       unsigned lhcBeamMode = gtRecord->gtfeWord().beamMode();
00267 
00268       if( lhcBeamMode == STABLE_BEAMS ){
00269 
00270         reset();
00271 
00272         stableBeamsDeclared_ = true;
00273 
00274       }
00275     }
00276   }
00277 
00278   edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00279   c.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00280 
00281   edm::Handle<EcalRecHitCollection> hits;
00282 
00283   if ( e.getByLabel(EcalRecHitCollection_, hits) ) {
00284 
00285     int neh = hits->size();
00286     LogDebug("EBTimingTask") << "event " << ievt_ << " hits collection size " << neh;
00287 
00288     for ( EcalRecHitCollection::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr ) {
00289 
00290       EBDetId id = hitItr->id();
00291 
00292       int ic = id.ic();
00293       int ie = (ic-1)/20 + 1;
00294       int ip = (ic-1)%20 + 1;
00295 
00296       int ism = Numbers::iSM( id );
00297 
00298       float xie = ie - 0.5;
00299       float xip = ip - 0.5;
00300 
00301       if ( isData ) {
00302 
00303         if ( ! ( runType[ism-1] == EcalDCCHeaderBlock::COSMIC ||
00304                  runType[ism-1] == EcalDCCHeaderBlock::MTCC ||
00305                  runType[ism-1] == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00306                  runType[ism-1] == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00307                  runType[ism-1] == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00308                  runType[ism-1] == EcalDCCHeaderBlock::PHYSICS_LOCAL ) ) continue;
00309 
00310       }
00311 
00312       MonitorElement* meTime = 0;
00313       MonitorElement* meTimeMap = 0;
00314       MonitorElement* meTimeAmpli = 0;
00315 
00316       meTime = meTime_[ism-1];
00317       meTimeMap = meTimeMap_[ism-1];
00318       meTimeAmpli = meTimeAmpli_[ism-1];
00319 
00320       float xval = hitItr->energy();
00321       float yval = hitItr->time();
00322 
00323       uint32_t flag = hitItr->recoFlag();
00324 
00325       uint32_t sev = sevlv->severityLevel(id, *hits);
00326 
00327       if ( (flag == EcalRecHit::kGood || flag == EcalRecHit::kOutOfTime) && sev != EcalSeverityLevel::kWeird ) {
00328         if ( meTimeAmpli ) meTimeAmpli->Fill(xval, yval);
00329         if ( meTimeAmpliSummary_ ) meTimeAmpliSummary_->Fill(xval, yval);
00330 
00331         if ( xval > energyThreshold_ ) {
00332           if ( meTime ) meTime->Fill(yval);
00333           if ( meTimeMap ) meTimeMap->Fill(xie, xip, yval+shiftProf2D_);
00334           if ( meTimeSummary1D_ ) meTimeSummary1D_->Fill(yval);
00335 
00336           float xebeta = id.ieta() - 0.5 * id.zside();
00337           float xebphi = id.iphi() - 0.5;
00338           if ( meTimeSummaryMap_ ) meTimeSummaryMap_->Fill(xebphi, xebeta, yval+shiftProf2D_);
00339         }
00340 
00341       }
00342     }
00343 
00344   } else {
00345 
00346     edm::LogWarning("EBTimingTask") << EcalRecHitCollection_ << " not available";
00347 
00348   }
00349 
00350 }
00351