CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/DQM/EcalEndcapMonitorTasks/src/EECosmicTask.cc

Go to the documentation of this file.
00001 /*
00002  * \file EECosmicTask.cc
00003  *
00004  * $Date: 2012/04/27 13:46:14 $
00005  * $Revision: 1.60 $
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 
00016 #include "DQMServices/Core/interface/MonitorElement.h"
00017 
00018 #include "DQMServices/Core/interface/DQMStore.h"
00019 
00020 #include "DataFormats/EcalRawData/interface/EcalRawDataCollections.h"
00021 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00022 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
00023 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00024 
00025 #include "DQM/EcalCommon/interface/Numbers.h"
00026 
00027 #include "DQM/EcalEndcapMonitorTasks/interface/EECosmicTask.h"
00028 
00029 EECosmicTask::EECosmicTask(const edm::ParameterSet& ps){
00030 
00031   init_ = false;
00032 
00033   dqmStore_ = edm::Service<DQMStore>().operator->();
00034 
00035   prefixME_ = ps.getUntrackedParameter<std::string>("prefixME", "");
00036 
00037   enableCleanup_ = ps.getUntrackedParameter<bool>("enableCleanup", false);
00038 
00039   mergeRuns_ = ps.getUntrackedParameter<bool>("mergeRuns", false);
00040 
00041   EcalRawDataCollection_ = ps.getParameter<edm::InputTag>("EcalRawDataCollection");
00042   EcalUncalibratedRecHitCollection_ = ps.getParameter<edm::InputTag>("EcalUncalibratedRecHitCollection");
00043   EcalRecHitCollection_ = ps.getParameter<edm::InputTag>("EcalRecHitCollection");
00044 
00045   threshold_ = 0.12500; // typical muon energy deposit is 250 MeV
00046 
00047   minJitter_ = -2.0;
00048   maxJitter_ =  1.5;
00049 
00050   for (int i = 0; i < 18; i++) {
00051     meSelMap_[i] = 0;
00052     meSpectrum_[0][i] = 0;
00053     meSpectrum_[1][i] = 0;
00054   }
00055 
00056 }
00057 
00058 EECosmicTask::~EECosmicTask(){
00059 
00060 }
00061 
00062 void EECosmicTask::beginJob(void){
00063 
00064   ievt_ = 0;
00065 
00066   if ( dqmStore_ ) {
00067     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask");
00068     dqmStore_->rmdir(prefixME_ + "/EECosmicTask");
00069   }
00070 
00071 }
00072 
00073 void EECosmicTask::beginRun(const edm::Run& r, const edm::EventSetup& c) {
00074 
00075   Numbers::initGeometry(c, false);
00076 
00077   if ( ! mergeRuns_ ) this->reset();
00078 
00079 }
00080 
00081 void EECosmicTask::endRun(const edm::Run& r, const edm::EventSetup& c) {
00082 
00083 }
00084 
00085 void EECosmicTask::reset(void) {
00086 
00087   for (int i = 0; i < 18; i++) {
00088     if ( meSelMap_[i] ) meSelMap_[i]->Reset();
00089     if ( meSpectrum_[0][i] ) meSpectrum_[0][i]->Reset();
00090     if ( meSpectrum_[1][i] ) meSpectrum_[1][i]->Reset();
00091   }
00092 
00093 }
00094 
00095 void EECosmicTask::setup(void){
00096 
00097   init_ = true;
00098 
00099   std::string name;
00100 
00101   if ( dqmStore_ ) {
00102     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask");
00103 
00104     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask/Sel");
00105     for (int i = 0; i < 18; i++) {
00106       name = "EECT energy sel " + Numbers::sEE(i+1);
00107       meSelMap_[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., 4096, 0., 4096., "s");
00108       meSelMap_[i]->setAxisTitle("ix", 1);
00109       if ( i+1 >= 1 && i+1 <= 9 ) meSelMap_[i]->setAxisTitle("101-ix", 1);
00110       meSelMap_[i]->setAxisTitle("iy", 2);
00111       meSelMap_[i]->setAxisTitle("energy (GeV)", 3);
00112     }
00113 
00114     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask/Spectrum");
00115     for (int i = 0; i < 18; i++) {
00116       name = "EECT 1x1 energy spectrum " + Numbers::sEE(i+1);
00117       meSpectrum_[0][i] = dqmStore_->book1D(name, name, 100, 0., 1.5);
00118       meSpectrum_[0][i]->setAxisTitle("energy (GeV)", 1);
00119 
00120       name = "EECT 3x3 energy spectrum " + Numbers::sEE(i+1);
00121       meSpectrum_[1][i] = dqmStore_->book1D(name, name, 100, 0., 1.5);
00122       meSpectrum_[1][i]->setAxisTitle("energy (GeV)", 1);
00123     }
00124 
00125   }
00126 
00127 }
00128 
00129 void EECosmicTask::cleanup(void){
00130 
00131   if ( ! init_ ) return;
00132 
00133   if ( dqmStore_ ) {
00134     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask");
00135 
00136     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask/Sel");
00137     for (int i = 0; i < 18; i++) {
00138       if ( meSelMap_[i] ) dqmStore_->removeElement( meSelMap_[i]->getName() );
00139       meSelMap_[i] = 0;
00140     }
00141 
00142     dqmStore_->setCurrentFolder(prefixME_ + "/EECosmicTask/Spectrum");
00143     for (int i = 0; i < 18; i++) {
00144       if ( meSpectrum_[0][i] ) dqmStore_->removeElement( meSpectrum_[0][i]->getName() );
00145       meSpectrum_[0][i] = 0;
00146       if ( meSpectrum_[1][i] ) dqmStore_->removeElement( meSpectrum_[1][i]->getName() );
00147       meSpectrum_[1][i] = 0;
00148     }
00149 
00150   }
00151 
00152   init_ = false;
00153 
00154 }
00155 
00156 void EECosmicTask::endJob(void){
00157 
00158   edm::LogInfo("EECosmicTask") << "analyzed " << ievt_ << " events";
00159 
00160   if ( enableCleanup_ ) this->cleanup();
00161 
00162 }
00163 
00164 void EECosmicTask::analyze(const edm::Event& e, const edm::EventSetup& c){
00165 
00166   bool isData = true;
00167   bool enable = false;
00168   int runType[18];
00169   for (int i=0; i<18; i++) runType[i] = -1;
00170 
00171   edm::Handle<EcalRawDataCollection> dcchs;
00172 
00173   if ( e.getByLabel(EcalRawDataCollection_, dcchs) ) {
00174 
00175     for ( EcalRawDataCollection::const_iterator dcchItr = dcchs->begin(); dcchItr != dcchs->end(); ++dcchItr ) {
00176 
00177       if ( Numbers::subDet( *dcchItr ) != EcalEndcap ) continue;
00178 
00179       int ism = Numbers::iSM( *dcchItr, EcalEndcap );
00180 
00181       runType[ism-1] = dcchItr->getRunType();
00182 
00183       if ( dcchItr->getRunType() == EcalDCCHeaderBlock::COSMIC ||
00184            dcchItr->getRunType() == EcalDCCHeaderBlock::MTCC ||
00185            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00186            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00187            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00188            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_LOCAL ) enable = true;
00189 
00190     }
00191 
00192   } else {
00193 
00194     isData = false; enable = true;
00195     edm::LogWarning("EECosmicTask") << EcalRawDataCollection_ << " not available";
00196 
00197   }
00198 
00199   if ( ! enable ) return;
00200 
00201   if ( ! init_ ) this->setup();
00202 
00203   ievt_++;
00204 
00205   edm::Handle<EcalRecHitCollection> hits;
00206 
00207   if ( e.getByLabel(EcalRecHitCollection_, hits) ) {
00208 
00209     int neeh = hits->size();
00210     LogDebug("EECosmicTask") << "event " << ievt_ << " hits collection size " << neeh;
00211 
00212     edm::Handle<EcalUncalibratedRecHitCollection> uhits;
00213 
00214     if ( ! e.getByLabel(EcalUncalibratedRecHitCollection_, uhits) ) {
00215       edm::LogWarning("EECosmicTask") << EcalUncalibratedRecHitCollection_ << " not available";
00216     }
00217 
00218     for ( EcalRecHitCollection::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr ) {
00219 
00220       EEDetId id = hitItr->id();
00221 
00222       int ix = id.ix();
00223       int iy = id.iy();
00224 
00225       int ism = Numbers::iSM( id );
00226 
00227       if ( ism >= 1 && ism <= 9 ) ix = 101 - ix;
00228 
00229       float xix = ix - 0.5;
00230       float xiy = iy - 0.5;
00231 
00232       int iz = 0;
00233 
00234       if ( ism >=  1 && ism <=  9 ) iz = -1;
00235       if ( ism >= 10 && ism <= 18 ) iz = +1;
00236 
00237       if ( isData ) {
00238 
00239         if ( ! ( runType[ism-1] == EcalDCCHeaderBlock::COSMIC ||
00240                  runType[ism-1] == EcalDCCHeaderBlock::MTCC ||
00241                  runType[ism-1] == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00242                  runType[ism-1] == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00243                  runType[ism-1] == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00244                  runType[ism-1] == EcalDCCHeaderBlock::PHYSICS_LOCAL ) ) continue;
00245 
00246       }
00247 
00248       float xval = hitItr->energy();
00249       if ( xval <= 0. ) xval = 0.0;
00250 
00251       // look for the seeds
00252       float e3x3 = 0.;
00253       bool isSeed = true;
00254 
00255       // evaluate 3x3 matrix around a seed
00256       for(int icry=0; icry<9; ++icry) {
00257         unsigned int row    = icry/3;
00258         unsigned int column = icry%3;
00259         int icryX = id.ix()+column-1;
00260         int icryY = id.iy()+row-1;
00261         if ( EEDetId::validDetId(icryX, icryY, iz) ) {
00262           EEDetId id3x3 = EEDetId(icryX, icryY, iz, EEDetId::XYMODE);
00263           if ( hits->find(id3x3) != hits->end() ) {
00264             float neighbourEnergy = hits->find(id3x3)->energy();
00265             e3x3 += neighbourEnergy;
00266             if ( neighbourEnergy > xval ) isSeed = false;
00267           }
00268         }
00269       }
00270 
00271       // find the jitter of the seed
00272       float jitter = -999.;
00273       if ( isSeed ) {
00274         if ( uhits.isValid() ) {
00275           if ( uhits->find(id) != uhits->end() ) {
00276             jitter = uhits->find(id)->jitter();
00277           }
00278         }
00279       }
00280 
00281       if ( isSeed && e3x3 >= threshold_ && jitter > minJitter_ && jitter < maxJitter_ ) {
00282         if ( meSelMap_[ism-1] ) meSelMap_[ism-1]->Fill(xix, xiy, e3x3);
00283       }
00284 
00285       if ( meSpectrum_[0][ism-1] ) meSpectrum_[0][ism-1]->Fill(xval);
00286 
00287       if ( isSeed && xval >= threshold_ && jitter > minJitter_ && jitter < maxJitter_ ) {
00288         if ( meSpectrum_[1][ism-1] ) meSpectrum_[1][ism-1]->Fill(e3x3);
00289       }
00290 
00291     }
00292 
00293   } else {
00294 
00295     edm::LogWarning("EECosmicTask") << EcalRecHitCollection_ << " not available";
00296 
00297   }
00298 
00299 }
00300