CMSSW Reference Manual

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