Reference Manual

Go to the documentation of this file.
00001 /*
00002  * \file EETimingTask.cc
00003  *
00004  * $Date: 2012/04/27 13:46:16 $
00005  * $Revision: 1.85 $
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 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerEvmReadoutRecord.h"
00026 
00027 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
00028 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
00029 
00030 #include "DQM/EcalCommon/interface/Numbers.h"
00031 
00032 #include "DQM/EcalEndcapMonitorTasks/interface/EETimingTask.h"
00033 
00034 EETimingTask::EETimingTask(const edm::ParameterSet& ps){
00035 
00036   init_ = false;
00037 
00038   initCaloGeometry_ = 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>("energyThreshold", 3.);
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 < 18; i++) {
00058     meTime_[i] = 0;
00059     meTimeMap_[i] = 0;
00060     meTimeAmpli_[i] = 0;
00061   }
00062 
00063   for (int i = 0; i < 2; i++) {
00064     meTimeAmpliSummary_[i] = 0;
00065     meTimeSummary1D_[i] = 0;
00066     meTimeSummaryMap_[i] = 0;
00067   }
00068 
00069   meTimeDelta_ = 0;
00070   meTimeDelta2D_ = 0;
00071 
00072   stableBeamsDeclared_ = false;
00073 
00074 }
00075 
00076 EETimingTask::~EETimingTask(){
00077 
00078 }
00079 
00080 void EETimingTask::beginJob(void){
00081 
00082   ievt_ = 0;
00083 
00084   if ( dqmStore_ ) {
00085     dqmStore_->setCurrentFolder(prefixME_ + "/EETimingTask");
00086     dqmStore_->rmdir(prefixME_ + "/EETimingTask");
00087   }
00088 
00089 }
00090 
00091 void EETimingTask::beginRun(const edm::Run& r, const edm::EventSetup& c) {
00092 
00093   Numbers::initGeometry(c, true);
00094 
00095   if( !initCaloGeometry_ ) {
00096     c.get<CaloGeometryRecord>().get(pGeometry_);
00097     initCaloGeometry_ = true;
00098   }
00099 
00100   if ( ! mergeRuns_ ) this->reset();
00101 
00102   stableBeamsDeclared_ = false;
00103 
00104 }
00105 
00106 void EETimingTask::endRun(const edm::Run& r, const edm::EventSetup& c) {
00107 
00108 }
00109 
00110 void EETimingTask::reset(void) {
00111 
00112   for (int i = 0; i < 18; i++) {
00113     if ( meTime_[i] ) meTime_[i]->Reset();
00114     if ( meTimeMap_[i] ) meTimeMap_[i]->Reset();
00115     if ( meTimeAmpli_[i] ) meTimeAmpli_[i]->Reset();
00116   }
00117 
00118   for (int i = 0; i < 2; i++) {
00119     if ( meTimeAmpliSummary_[i] ) meTimeAmpliSummary_[i]->Reset();
00120     if ( meTimeSummary1D_[i] ) meTimeSummary1D_[i]->Reset();
00121     if ( meTimeSummaryMap_[i] ) meTimeSummaryMap_[i]->Reset();
00122   }
00123 
00124   if ( meTimeDelta_ ) meTimeDelta_->Reset();
00125   if ( meTimeDelta2D_ ) meTimeDelta2D_->Reset();
00126 
00127 }
00128 
00129 void EETimingTask::setup(void){
00130 
00131   init_ = true;
00132 
00133   std::string name;
00134 
00135   //for timing vs amplitude plots
00136   const int nbinsE = 25;
00137   const float minlogE = -0.5;
00138   const float maxlogE = 2.;
00139   float binEdgesE[nbinsE + 1];
00140   for(int i = 0; i <= nbinsE; i++)
00141     binEdgesE[i] = std::pow((float)10., minlogE + (maxlogE - minlogE) / nbinsE * i);
00142 
00143   const int nbinsT = 200;
00144   const float minT = -50.;
00145   const float maxT = 50.;
00146   float binEdgesT[nbinsT + 1];
00147   for(int i = 0; i <= nbinsT; i++)
00148     binEdgesT[i] = minT + (maxT - minT) / nbinsT * i;
00149 
00150   if ( dqmStore_ ) {
00151     dqmStore_->setCurrentFolder(prefixME_ + "/EETimingTask");
00152 
00153     for (int i = 0; i < 18; i++) {
00154       name = "EETMT timing 1D " + Numbers::sEE(i+1);
00155       meTime_[i] = dqmStore_->book1D(name, name, 50, -25., 25.);
00156       meTime_[i]->setAxisTitle("time (ns)", 1);
00157       dqmStore_->tag(meTime_[i], i+1);
00158 
00159       name = "EETMT timing " + Numbers::sEE(i+1);
00160       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");
00161       meTimeMap_[i]->setAxisTitle("ix", 1);
00162       if ( i+1 >= 1 && i+1 <= 9 ) meTimeMap_[i]->setAxisTitle("101-ix", 1);
00163       meTimeMap_[i]->setAxisTitle("iy", 2);
00164       meTimeMap_[i]->setAxisTitle("time (ns)", 3);
00165       dqmStore_->tag(meTimeMap_[i], i+1);
00166 
00167       name = "EETMT timing vs amplitude " + Numbers::sEE(i+1);
00168       meTimeAmpli_[i] = dqmStore_->book2D(name, name, nbinsE, binEdgesE, nbinsT, binEdgesT);
00169       meTimeAmpli_[i]->setAxisTitle("energy (GeV)", 1);
00170       meTimeAmpli_[i]->setAxisTitle("time (ns)", 2);
00171       dqmStore_->tag(meTimeAmpli_[i], i+1);
00172     }
00173 
00174     name = "EETMT timing vs amplitude summary EE -";
00175     meTimeAmpliSummary_[0] = dqmStore_->book2D(name, name, nbinsE, binEdgesE, nbinsT, binEdgesT);
00176     meTimeAmpliSummary_[0]->setAxisTitle("energy (GeV)", 1);
00177     meTimeAmpliSummary_[0]->setAxisTitle("time (ns)", 2);
00178 
00179     name = "EETMT timing vs amplitude summary EE +";
00180     meTimeAmpliSummary_[1] = dqmStore_->book2D(name, name, nbinsE, binEdgesE, nbinsT, binEdgesT);
00181     meTimeAmpliSummary_[1]->setAxisTitle("energy (GeV)", 1);
00182     meTimeAmpliSummary_[1]->setAxisTitle("time (ns)", 2);
00183 
00184     name = "EETMT timing 1D summary EE -";
00185     meTimeSummary1D_[0] = dqmStore_->book1D(name, name, 50, -25., 25.);
00186     meTimeSummary1D_[0]->setAxisTitle("time (ns)", 1);
00187 
00188     name = "EETMT timing 1D summary EE +";
00189     meTimeSummary1D_[1] = dqmStore_->book1D(name, name, 50, -25., 25.);
00190     meTimeSummary1D_[1]->setAxisTitle("time (ns)", 1);
00191 
00192     name = "EETMT timing map EE -";
00193     meTimeSummaryMap_[0] = dqmStore_->bookProfile2D(name, name, 20, 0., 100., 20, 0., 100., -20.+shiftProf2D, 20.+shiftProf2D, "s");
00194     meTimeSummaryMap_[0]->setAxisTitle("ix'", 1);
00195     meTimeSummaryMap_[0]->setAxisTitle("101-iy'", 2);
00196     meTimeSummaryMap_[0]->setAxisTitle("time (ns)", 3);
00197 
00198     name = "EETMT timing map EE +";
00199     meTimeSummaryMap_[1] = dqmStore_->bookProfile2D(name, name, 20, 0., 100., 20, 0., 100., -20.+shiftProf2D, 20.+shiftProf2D, "s");
00200     meTimeSummaryMap_[1]->setAxisTitle("ix'", 1);
00201     meTimeSummaryMap_[1]->setAxisTitle("iy'", 2);
00202     meTimeSummaryMap_[1]->setAxisTitle("time (ns)", 3);
00203 
00204     name = "EETMT timing EE+ - EE-";
00205     meTimeDelta_ = dqmStore_->book1D(name, name, 100, -3., 3.);
00206     meTimeDelta_->setAxisTitle("time (ns)", 1);
00207 
00208     name = "EETMT timing EE+ vs EE-";
00209     meTimeDelta2D_ = dqmStore_->book2D(name, name, 50, -25., 25., 50, -25., 25.);
00210     meTimeDelta2D_->setAxisTitle("EE+ average time (ns)", 1);
00211     meTimeDelta2D_->setAxisTitle("EE- average time (ns)", 2);
00212 
00213   }
00214 
00215 }
00216 
00217 void EETimingTask::cleanup(void){
00218 
00219   if ( ! init_ ) return;
00220 
00221   if ( dqmStore_ ) {
00222     dqmStore_->setCurrentFolder(prefixME_ + "/EETimingTask");
00223 
00224     for ( int i = 0; i < 18; i++ ) {
00225       if ( meTime_[i] ) dqmStore_->removeElement( meTime_[i]->getName() );
00226       meTime_[i] = 0;
00227 
00228       if ( meTimeMap_[i] ) dqmStore_->removeElement( meTimeMap_[i]->getName() );
00229       meTimeMap_[i] = 0;
00230 
00231       if ( meTimeAmpli_[i] ) dqmStore_->removeElement( meTimeAmpli_[i]->getName() );
00232       meTimeAmpli_[i] = 0;
00233     }
00234 
00235     for (int i = 0; i < 2; i++) {
00236       if ( meTimeAmpliSummary_[i] ) dqmStore_->removeElement( meTimeAmpliSummary_[i]->getName() );
00237       meTimeAmpliSummary_[i] = 0;
00238 
00239       if ( meTimeSummary1D_[i] ) dqmStore_->removeElement( meTimeSummary1D_[i]->getName() );
00240       meTimeSummary1D_[i] = 0;
00241 
00242       if ( meTimeSummaryMap_[i] ) dqmStore_->removeElement( meTimeSummaryMap_[i]->getName() );
00243       meTimeSummaryMap_[i] = 0;
00244 
00245     }
00246 
00247     if ( meTimeDelta_ ) dqmStore_->removeElement( meTimeDelta_->getName() );
00248     meTimeDelta_ = 0;
00249 
00250     if ( meTimeDelta2D_ ) dqmStore_->removeElement( meTimeDelta2D_->getName() );
00251     meTimeDelta2D_ = 0;
00252 
00253   }
00254 
00255   init_ = false;
00256 
00257 }
00258 
00259 void EETimingTask::endJob(void){
00260 
00261   edm::LogInfo("EETimingTask") << "analyzed " << ievt_ << " events";
00262 
00263   if ( enableCleanup_ ) this->cleanup();
00264 
00265 }
00266 
00267 void EETimingTask::analyze(const edm::Event& e, const edm::EventSetup& c){
00268 
00269   const unsigned STABLE_BEAMS = 11;
00270 
00271   bool isData = true;
00272   bool enable = false;
00273   int runType[18];
00274   for (int i=0; i<18; i++) runType[i] = -1;
00275 
00276   edm::Handle<EcalRawDataCollection> dcchs;
00277 
00278   if ( e.getByLabel(EcalRawDataCollection_, dcchs) ) {
00279 
00280     for ( EcalRawDataCollection::const_iterator dcchItr = dcchs->begin(); dcchItr != dcchs->end(); ++dcchItr ) {
00281 
00282       if ( Numbers::subDet( *dcchItr ) != EcalEndcap ) continue;
00283 
00284       int ism = Numbers::iSM( *dcchItr, EcalEndcap );
00285 
00286       runType[ism-1] = dcchItr->getRunType();
00287 
00288       if ( dcchItr->getRunType() == EcalDCCHeaderBlock::COSMIC ||
00289            dcchItr->getRunType() == EcalDCCHeaderBlock::MTCC ||
00290            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00291            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00292            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00293            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_LOCAL ) enable = true;
00294 
00295     }
00296 
00297   } else {
00298 
00299     isData = false; enable = true;
00300     edm::LogWarning("EETimingTask") << EcalRawDataCollection_ << " not available";
00301 
00302   }
00303 
00304   if ( ! enable ) return;
00305 
00306   if ( ! init_ ) this->setup();
00307 
00308   ievt_++;
00309 
00310   // resetting plots when stable beam is declared
00311   if( useBeamStatus_ && !stableBeamsDeclared_ ) {
00312     edm::Handle<L1GlobalTriggerEvmReadoutRecord> gtRecord;
00313     if( e.getByLabel(L1GtEvmReadoutRecord_, gtRecord) ) {
00314 
00315       unsigned lhcBeamMode = gtRecord->gtfeWord().beamMode();
00316 
00317       if( lhcBeamMode == STABLE_BEAMS ){
00318 
00319         reset();
00320 
00321         stableBeamsDeclared_ = true;
00322 
00323       }
00324     }
00325   }
00326 
00327   float sumTime_hithr[2] = {0.,0.};
00328   int n_hithr[2] = {0,0};
00329 
00330   edm::ESHandle<EcalSeverityLevelAlgo> sevlv;
00331   c.get<EcalSeverityLevelAlgoRcd>().get(sevlv);
00332 
00333   edm::Handle<EcalRecHitCollection> hits;
00334 
00335   if ( e.getByLabel(EcalRecHitCollection_, hits) ) {
00336 
00337     int neh = hits->size();
00338     LogDebug("EETimingTask") << "event " << ievt_ << " hits collection size " << neh;
00339 
00340     for ( EcalRecHitCollection::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr ) {
00341 
00342       EEDetId id = hitItr->id();
00343 
00344       int ix = id.ix();
00345       int iy = id.iy();
00346       int iz = ( id.positiveZ() ) ? 1 : 0;
00347 
00348       int ism = Numbers::iSM( id );
00349 
00350       if ( ism >= 1 && ism <= 9 ) ix = 101 - ix;
00351 
00352       float xix = ix - 0.5;
00353       float xiy = iy - 0.5;
00354 
00355       if ( isData ) {
00356 
00357         if ( ! ( runType[ism-1] == EcalDCCHeaderBlock::COSMIC ||
00358                  runType[ism-1] == EcalDCCHeaderBlock::MTCC ||
00359                  runType[ism-1] == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00360                  runType[ism-1] == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00361                  runType[ism-1] == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00362                  runType[ism-1] == EcalDCCHeaderBlock::PHYSICS_LOCAL ) ) continue;
00363 
00364       }
00365 
00366       MonitorElement* meTime = 0;
00367       MonitorElement* meTimeMap = 0;
00368       MonitorElement* meTimeAmpli = 0;
00369 
00370       meTime = meTime_[ism-1];
00371       meTimeMap = meTimeMap_[ism-1];
00372       meTimeAmpli = meTimeAmpli_[ism-1];
00373 
00374       float xval = hitItr->energy();
00375       float yval = hitItr->time();
00376 
00377       uint32_t flag = hitItr->recoFlag();
00378 
00379       uint32_t sev = sevlv->severityLevel(id, *hits );
00380 
00381       if ( (flag == EcalRecHit::kGood || flag == EcalRecHit::kOutOfTime) && sev != EcalSeverityLevel::kWeird ) {
00382         if ( meTimeAmpli ) meTimeAmpli->Fill(xval, yval);
00383         if ( meTimeAmpliSummary_[iz] ) meTimeAmpliSummary_[iz]->Fill(xval, yval);
00384         if ( hitItr->energy() > energyThreshold_ ) {
00385           if ( meTimeMap ) meTimeMap->Fill(xix, xiy, yval+shiftProf2D);
00386           if ( meTime ) meTime->Fill(yval);
00387           if ( meTimeSummary1D_[iz] ) meTimeSummary1D_[iz]->Fill(yval);
00388 
00389           if ( meTimeSummaryMap_[iz] ) meTimeSummaryMap_[iz]->Fill(id.ix()-0.5, xiy, yval+shiftProf2D);
00390 
00391           sumTime_hithr[iz] += yval;
00392           n_hithr[iz]++;
00393         }
00394       } // good rh for timing
00395     } // loop over rh
00396 
00397     if (n_hithr[0] > 0 && n_hithr[1] > 0 ) {
00398       if ( meTimeDelta_ ) meTimeDelta_->Fill( sumTime_hithr[1]/n_hithr[1] - sumTime_hithr[0]/n_hithr[0] );
00399       if ( meTimeDelta2D_ ) meTimeDelta2D_->Fill( sumTime_hithr[1]/n_hithr[1], sumTime_hithr[0]/n_hithr[0] );
00400     }
00401 
00402   } else {
00403 
00404     edm::LogWarning("EETimingTask") << EcalRecHitCollection_ << " not available";
00405 
00406   }
00407 
00408 }
00409
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_14/src/DQM/EcalEndcapMonitorTasks/src/EETimingTask.cc
