Reference Manual

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