Reference Manual

Go to the documentation of this file.
00001 /*
00002  * \file EBClusterTask.cc
00003  *
00004  * $Date: 2010/11/10 16:59:50 $
00005  * $Revision: 1.91 $
00006  * \author G. Della Ricca
00007  * \author E. Di Marco
00008  *
00009 */
00010 
00011 #include <iostream>
00012 #include <fstream>
00013 #include <vector>
00014 #include <math.h>
00015 
00016 #include "FWCore/ServiceRegistry/interface/Service.h"
00017 #include "FWCore/MessageLogger/interface/MessageLogger.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/EgammaReco/interface/BasicCluster.h"
00025 #include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
00026 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
00027 #include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
00028 #include "FWCore/Framework/interface/ESHandle.h"
00029 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00030 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00031 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
00032 #include "DataFormats/Math/interface/Point3D.h"
00033 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
00034 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
00035 
00036 #include "DQM/EcalCommon/interface/Numbers.h"
00037 
00038 #include "DQM/EcalBarrelMonitorTasks/interface/EBClusterTask.h"
00039 
00040 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00041 
00042 #include "TLorentzVector.h"
00043 
00044 EBClusterTask::EBClusterTask(const edm::ParameterSet& ps){
00045 
00046   init_ = false;
00047 
00048   dqmStore_ = edm::Service<DQMStore>().operator->();
00049 
00050   prefixME_ = ps.getUntrackedParameter<std::string>("prefixME", "");
00051 
00052   enableCleanup_ = ps.getUntrackedParameter<bool>("enableCleanup", false);
00053 
00054   mergeRuns_ = ps.getUntrackedParameter<bool>("mergeRuns", false);
00055 
00056   // parameters...
00057   EcalRawDataCollection_ = ps.getParameter<edm::InputTag>("EcalRawDataCollection");
00058   BasicClusterCollection_ = ps.getParameter<edm::InputTag>("BasicClusterCollection");
00059   SuperClusterCollection_ = ps.getParameter<edm::InputTag>("SuperClusterCollection");
00060   EcalRecHitCollection_ = ps.getParameter<edm::InputTag>("EcalRecHitCollection");
00061 
00062   // histograms...
00063   meBCEne_ = 0;
00064   meBCNum_ = 0;
00065   meBCSiz_ = 0;
00066 
00067   meBCEneMap_ = 0;
00068   meBCNumMap_ = 0;
00069   meBCETMap_  = 0;
00070   meBCSizMap_ = 0;
00071 
00072   meBCEneMapProjEta_ = 0;
00073   meBCNumMapProjEta_ = 0;
00074   meBCETMapProjEta_  = 0;
00075   meBCSizMapProjEta_ = 0;
00076 
00077   meBCEneMapProjPhi_ = 0;
00078   meBCNumMapProjPhi_ = 0;
00079   meBCETMapProjPhi_  = 0;
00080   meBCSizMapProjPhi_ = 0;
00081 
00082   meSCEne_ = 0;
00083   meSCNum_ = 0;
00084   meSCSiz_ = 0;
00085 
00086   meSCCrystalSiz_ = 0;
00087   meSCSeedEne_ = 0;
00088   meSCEne2_ = 0;
00089   meSCEneVsEMax_ = 0;
00090   meSCEneLowScale_ = 0;
00091   meSCSeedMapOcc_ = 0;
00092   meSCMapSingleCrystal_ = 0;
00093 
00094   mes1s9_  = 0;
00095   mes1s9thr_  = 0;
00096   mes9s25_  = 0;
00097 
00098   meInvMassPi0_ = 0;
00099   meInvMassJPsi_ = 0;
00100   meInvMassZ0_ = 0;
00101   meInvMassHigh_ = 0;
00102 
00103   meInvMassPi0Sel_ = 0;
00104   meInvMassJPsiSel_ = 0;
00105   meInvMassZ0Sel_ = 0;
00106   meInvMassHighSel_ = 0;
00107 
00108   thrS4S9_ = 0.85;
00109   thrClusEt_ = 0.200;
00110   thrCandEt_ = 0.650;
00111 
00112 }
00113 
00114 EBClusterTask::~EBClusterTask(){
00115 
00116 }
00117 
00118 void EBClusterTask::beginJob(void){
00119 
00120   ievt_ = 0;
00121 
00122   if ( dqmStore_ ) {
00123     dqmStore_->setCurrentFolder(prefixME_ + "/EBClusterTask");
00124     dqmStore_->rmdir(prefixME_ + "/EBClusterTask");
00125   }
00126 
00127 }
00128 
00129 void EBClusterTask::beginRun(const edm::Run& r, const edm::EventSetup& c) {
00130 
00131   Numbers::initGeometry(c, false);
00132 
00133   if ( ! mergeRuns_ ) this->reset();
00134 
00135 }
00136 
00137 void EBClusterTask::endRun(const edm::Run& r, const edm::EventSetup& c) {
00138 
00139 }
00140 
00141 void EBClusterTask::reset(void) {
00142 
00143   if ( meBCEne_ ) meBCEne_->Reset();
00144 
00145   if ( meBCNum_ ) meBCNum_->Reset();
00146 
00147   if ( meBCSiz_ ) meBCSiz_->Reset();
00148 
00149   if ( meBCEneMap_ ) meBCEneMap_->Reset();
00150 
00151   if ( meBCNumMap_ ) meBCNumMap_->Reset();
00152 
00153   if ( meBCETMap_ ) meBCETMap_->Reset();
00154 
00155   if ( meBCSizMap_ ) meBCSizMap_->Reset();
00156 
00157   if ( meBCEneMapProjEta_ ) meBCEneMapProjEta_->Reset();
00158 
00159   if ( meBCEneMapProjPhi_ ) meBCEneMapProjPhi_->Reset();
00160 
00161   if ( meBCNumMapProjEta_ ) meBCNumMapProjEta_->Reset();
00162 
00163   if ( meBCNumMapProjPhi_ ) meBCNumMapProjPhi_->Reset();
00164 
00165   if ( meBCETMapProjEta_ ) meBCETMapProjEta_->Reset();
00166 
00167   if ( meBCETMapProjPhi_ ) meBCETMapProjPhi_->Reset();
00168 
00169   if ( meBCSizMapProjEta_ ) meBCSizMapProjEta_->Reset();
00170 
00171   if ( meBCSizMapProjPhi_ ) meBCSizMapProjPhi_->Reset();
00172 
00173   if ( meSCEne_ ) meSCEne_->Reset();
00174 
00175   if ( meSCNum_ ) meSCNum_->Reset();
00176 
00177   if ( meSCSiz_ ) meSCSiz_->Reset();
00178 
00179   if ( meSCCrystalSiz_ ) meSCCrystalSiz_->Reset();
00180 
00181   if ( meSCSeedEne_ ) meSCSeedEne_->Reset();
00182 
00183   if ( meSCEne2_ ) meSCEne2_->Reset();
00184 
00185   if ( meSCEneVsEMax_ ) meSCEneVsEMax_->Reset();
00186 
00187   if ( meSCEneLowScale_ ) meSCEneLowScale_->Reset();
00188 
00189   if ( meSCSeedMapOcc_ ) meSCSeedMapOcc_->Reset();
00190 
00191   if ( meSCMapSingleCrystal_ ) meSCMapSingleCrystal_->Reset();
00192 
00193   if ( mes1s9_ ) mes1s9_->Reset();
00194 
00195   if ( mes1s9thr_ ) mes1s9thr_->Reset();
00196 
00197   if ( mes9s25_ ) mes9s25_->Reset();
00198 
00199   if ( meInvMassPi0_ ) meInvMassPi0_->Reset();
00200 
00201   if ( meInvMassJPsi_ ) meInvMassJPsi_->Reset();
00202 
00203   if ( meInvMassZ0_ ) meInvMassZ0_->Reset();
00204 
00205   if ( meInvMassHigh_ ) meInvMassHigh_->Reset();
00206 
00207   if ( meInvMassPi0Sel_ ) meInvMassPi0Sel_->Reset();
00208 
00209   if ( meInvMassJPsiSel_ ) meInvMassJPsiSel_->Reset();
00210 
00211   if ( meInvMassZ0Sel_ ) meInvMassZ0Sel_->Reset();
00212 
00213   if ( meInvMassHighSel_ ) meInvMassHighSel_->Reset();
00214 
00215 }
00216 
00217 void EBClusterTask::setup(void){
00218 
00219   init_ = true;
00220 
00221   char histo[200];
00222 
00223   if ( dqmStore_ ) {
00224     dqmStore_->setCurrentFolder(prefixME_ + "/EBClusterTask");
00225 
00226     sprintf(histo, "EBCLT BC energy");
00227     meBCEne_ = dqmStore_->book1D(histo, histo, 100, 0., 150.);
00228     meBCEne_->setAxisTitle("energy (GeV)", 1);
00229 
00230     sprintf(histo, "EBCLT BC number");
00231     meBCNum_ = dqmStore_->book1D(histo, histo, 100, 0., 100.);
00232     meBCNum_->setAxisTitle("number of clusters", 1);
00233 
00234     sprintf(histo, "EBCLT BC size");
00235     meBCSiz_ = dqmStore_->book1D(histo, histo, 100, 0., 100.);
00236     meBCSiz_->setAxisTitle("cluster size", 1);
00237 
00238     sprintf(histo, "EBCLT BC energy map");
00239     meBCEneMap_ = dqmStore_->bookProfile2D(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 34, -1.479, 1.479, 100, 0., 500., "s");
00240     meBCEneMap_->setAxisTitle("phi", 1);
00241     meBCEneMap_->setAxisTitle("eta", 2);
00242 
00243     sprintf(histo, "EBCLT BC number map");
00244     meBCNumMap_ = dqmStore_->book2D(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 34, -1.479, 1.479);
00245     meBCNumMap_->setAxisTitle("phi", 1);
00246     meBCNumMap_->setAxisTitle("eta", 2);
00247 
00248     sprintf(histo, "EBCLT BC ET map");
00249     meBCETMap_ = dqmStore_->bookProfile2D(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 34, -1.479, 1.479, 100, 0., 500., "s");
00250     meBCETMap_->setAxisTitle("phi", 1);
00251     meBCETMap_->setAxisTitle("eta", 2);
00252 
00253     sprintf(histo, "EBCLT BC size map");
00254     meBCSizMap_ = dqmStore_->bookProfile2D(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 34, -1.479, 1.479, 100, 0., 100., "s");
00255     meBCSizMap_->setAxisTitle("phi", 1);
00256     meBCSizMap_->setAxisTitle("eta", 2);
00257 
00258     sprintf(histo, "EBCLT BC energy projection eta");
00259     meBCEneMapProjEta_ = dqmStore_->bookProfile(histo, histo, 34, -1.479, 1.479, 100, 0., 500., "s");
00260     meBCEneMapProjEta_->setAxisTitle("eta", 1);
00261     meBCEneMapProjEta_->setAxisTitle("energy (GeV)", 2);
00262 
00263     sprintf(histo, "EBCLT BC energy projection phi");
00264     meBCEneMapProjPhi_ = dqmStore_->bookProfile(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 100, 0., 500., "s");
00265     meBCEneMapProjPhi_->setAxisTitle("phi", 1);
00266     meBCEneMapProjPhi_->setAxisTitle("energy (GeV)", 2);
00267 
00268     sprintf(histo, "EBCLT BC number projection eta");
00269     meBCNumMapProjEta_ = dqmStore_->book1D(histo, histo, 34, -1.479, 1.479);
00270     meBCNumMapProjEta_->setAxisTitle("eta", 1);
00271     meBCNumMapProjEta_->setAxisTitle("number of clusters", 2);
00272 
00273     sprintf(histo, "EBCLT BC number projection phi");
00274     meBCNumMapProjPhi_ = dqmStore_->book1D(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9);
00275     meBCNumMapProjPhi_->setAxisTitle("phi", 1);
00276     meBCNumMapProjPhi_->setAxisTitle("number of clusters", 2);
00277 
00278     sprintf(histo, "EBCLT BC ET projection eta");
00279     meBCETMapProjEta_ = dqmStore_->bookProfile(histo, histo, 34, -1.479, 1.479, 100, 0., 500., "s");
00280     meBCETMapProjEta_->setAxisTitle("eta", 1);
00281     meBCETMapProjEta_->setAxisTitle("transverse energy (GeV)", 2);
00282 
00283     sprintf(histo, "EBCLT BC ET projection phi");
00284     meBCETMapProjPhi_ = dqmStore_->bookProfile(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 100, 0., 500., "s");
00285     meBCETMapProjPhi_->setAxisTitle("phi", 1);
00286     meBCETMapProjPhi_->setAxisTitle("transverse energy (GeV)", 2);
00287 
00288     sprintf(histo, "EBCLT BC size projection eta");
00289     meBCSizMapProjEta_ = dqmStore_->bookProfile(histo, histo, 34, -1.479, 1.479, 100, 0., 100., "s");
00290     meBCSizMapProjEta_->setAxisTitle("eta", 1);
00291     meBCSizMapProjEta_->setAxisTitle("cluster size", 2);
00292 
00293     sprintf(histo, "EBCLT BC size projection phi");
00294     meBCSizMapProjPhi_ = dqmStore_->bookProfile(histo, histo, 72, -M_PI*(9+1.5)/9, M_PI*(9-1.5)/9, 100, 0., 100., "s");
00295     meBCSizMapProjPhi_->setAxisTitle("phi", 1);
00296     meBCSizMapProjPhi_->setAxisTitle("cluster size", 2);
00297 
00298     sprintf(histo, "EBCLT SC energy");
00299     meSCEne_ = dqmStore_->book1D(histo, histo, 100, 0., 150.);
00300     meSCEne_->setAxisTitle("energy (GeV)", 1);
00301 
00302     sprintf(histo, "EBCLT SC number");
00303     meSCNum_ = dqmStore_->book1D(histo, histo, 50, 0., 50.);
00304     meSCNum_->setAxisTitle("number of clusters", 1);
00305 
00306     sprintf(histo, "EBCLT SC size");
00307     meSCSiz_ = dqmStore_->book1D(histo, histo, 50, 0., 50.);
00308     meSCSiz_->setAxisTitle("cluster size", 1);
00309 
00310     sprintf(histo, "EBCLT SC size (crystal)");
00311     meSCCrystalSiz_ = dqmStore_->book1D(histo, histo, 150, 0, 150);
00312     meSCCrystalSiz_->setAxisTitle("cluster size in crystals", 1);
00313 
00314     sprintf(histo, "EBCLT SC seed crystal energy");
00315     meSCSeedEne_ = dqmStore_->book1D(histo, histo, 100, 0., 10.);
00316     meSCSeedEne_->setAxisTitle("seed crystal energy (GeV)", 1);
00317 
00318     sprintf(histo, "EBCLT SC e2");
00319     meSCEne2_ = dqmStore_->book1D(histo, histo, 100, 0., 10.);
00320     meSCEne2_->setAxisTitle("seed + highest neighbor crystal energy (GeV)", 1);
00321 
00322     sprintf(histo, "EBCLT SC energy vs seed crystal energy");
00323     meSCEneVsEMax_ = dqmStore_->book2D(histo, histo, 50, 0., 10., 50, 0., 10.);
00324     meSCEneVsEMax_->setAxisTitle("seed crystal energy (GeV)", 1);
00325     meSCEneVsEMax_->setAxisTitle("cluster energy (GeV)", 2);
00326 
00327     sprintf(histo, "EBCLT SC energy (low scale)");
00328     meSCEneLowScale_ = dqmStore_->book1D(histo, histo, 200, 0., 10.);
00329     meSCEneLowScale_->setAxisTitle("cluster energy (GeV)", 1);
00330 
00331     sprintf(histo, "EBCLT SC seed occupancy map");
00332     meSCSeedMapOcc_ = dqmStore_->book2D(histo, histo, 72, 0., 360., 34, -85, 85);
00333     meSCSeedMapOcc_->setAxisTitle("jphi", 1);
00334     meSCSeedMapOcc_->setAxisTitle("jeta", 2);
00335 
00336     sprintf(histo, "EBCLT SC single crystal cluster seed occupancy map");
00337     meSCMapSingleCrystal_ = dqmStore_->book2D(histo, histo, 72, 0., 360., 34, -85, 85);
00338     meSCMapSingleCrystal_->setAxisTitle("jphi", 1);
00339     meSCMapSingleCrystal_->setAxisTitle("jeta", 2);
00340 
00341     sprintf(histo, "EBCLT s1s9");
00342     mes1s9_ = dqmStore_->book1D(histo, histo, 50, 0., 1.5);
00343     mes1s9_->setAxisTitle("s1/s9", 1);
00344 
00345     sprintf(histo, "EBCLT s1s9 thr");
00346     mes1s9thr_ = dqmStore_->book1D(histo, histo, 50, 0., 1.5);
00347     mes1s9thr_->setAxisTitle("s1/s9", 1);
00348 
00349     sprintf(histo, "EBCLT s9s25");
00350     mes9s25_ = dqmStore_->book1D(histo, histo, 75, 0., 1.5);
00351     mes9s25_->setAxisTitle("s9/s25", 1);
00352 
00353     sprintf(histo, "EBCLT dicluster invariant mass Pi0");
00354     meInvMassPi0_ = dqmStore_->book1D(histo, histo, 50, 0.0, 0.500);
00355     meInvMassPi0_->setAxisTitle("mass (GeV)", 1);
00356 
00357     sprintf(histo, "EBCLT dicluster invariant mass JPsi");
00358     meInvMassJPsi_ = dqmStore_->book1D(histo, histo, 50, 2.9, 3.3);
00359     meInvMassJPsi_->setAxisTitle("mass (GeV)", 1);
00360 
00361     sprintf(histo, "EBCLT dicluster invariant mass Z0");
00362     meInvMassZ0_ = dqmStore_->book1D(histo, histo, 50, 40, 110);
00363     meInvMassZ0_->setAxisTitle("mass (GeV)", 1);
00364 
00365     sprintf(histo, "EBCLT dicluster invariant mass high");
00366     meInvMassHigh_ = dqmStore_->book1D(histo, histo, 500, 110, 3000);
00367     meInvMassHigh_->setAxisTitle("mass (GeV)", 1);
00368 
00369     sprintf(histo, "EBCLT dicluster invariant mass Pi0 sel");
00370     meInvMassPi0Sel_ = dqmStore_->book1D(histo, histo, 50, 0.00, 0.500);
00371     meInvMassPi0Sel_->setAxisTitle("mass (GeV)", 1);
00372 
00373     sprintf(histo, "EBCLT dicluster invariant mass JPsi sel");
00374     meInvMassJPsiSel_ = dqmStore_->book1D(histo, histo, 50, 2.9, 3.3);
00375     meInvMassJPsiSel_->setAxisTitle("mass (GeV)", 1);
00376 
00377     sprintf(histo, "EBCLT dicluster invariant mass Z0 sel");
00378     meInvMassZ0Sel_ = dqmStore_->book1D(histo, histo, 50, 40, 110);
00379     meInvMassZ0Sel_->setAxisTitle("mass (GeV)", 1);
00380 
00381     sprintf(histo, "EBCLT dicluster invariant mass high sel");
00382     meInvMassHighSel_ = dqmStore_->book1D(histo, histo, 500, 110, 3000);
00383     meInvMassHighSel_->setAxisTitle("mass (GeV)", 1);
00384 
00385   }
00386 
00387 }
00388 
00389 void EBClusterTask::cleanup(void){
00390 
00391   if ( ! init_ ) return;
00392 
00393   if ( dqmStore_ ) {
00394     dqmStore_->setCurrentFolder(prefixME_ + "/EBClusterTask");
00395 
00396     if ( meBCEne_ ) dqmStore_->removeElement( meBCEne_->getName() );
00397     meBCEne_ = 0;
00398 
00399     if ( meBCNum_ ) dqmStore_->removeElement( meBCNum_->getName() );
00400     meBCNum_ = 0;
00401 
00402     if ( meBCSiz_ ) dqmStore_->removeElement( meBCSiz_->getName() );
00403     meBCSiz_ = 0;
00404 
00405     if ( meBCEneMap_ ) dqmStore_->removeElement( meBCEneMap_->getName() );
00406     meBCEneMap_ = 0;
00407 
00408     if ( meBCNumMap_ ) dqmStore_->removeElement( meBCNumMap_->getName() );
00409     meBCNumMap_ = 0;
00410 
00411     if ( meBCETMap_ ) dqmStore_->removeElement( meBCETMap_->getName() );
00412     meBCETMap_ = 0;
00413 
00414     if ( meBCSizMap_ ) dqmStore_->removeElement( meBCSizMap_->getName() );
00415     meBCSizMap_ = 0;
00416 
00417     if ( meBCEneMapProjEta_ ) dqmStore_->removeElement( meBCEneMapProjEta_->getName() );
00418     meBCEneMapProjEta_ = 0;
00419 
00420     if ( meBCEneMapProjPhi_ ) dqmStore_->removeElement( meBCEneMapProjPhi_->getName() );
00421     meBCEneMapProjPhi_ = 0;
00422 
00423     if ( meBCNumMapProjEta_ ) dqmStore_->removeElement( meBCNumMapProjEta_->getName() );
00424     meBCNumMapProjEta_ = 0;
00425 
00426     if ( meBCNumMapProjPhi_ ) dqmStore_->removeElement( meBCNumMapProjPhi_->getName() );
00427     meBCNumMapProjPhi_ = 0;
00428 
00429     if ( meBCETMapProjEta_ ) dqmStore_->removeElement( meBCETMapProjEta_->getName() );
00430     meBCETMapProjEta_ = 0;
00431 
00432     if ( meBCETMapProjPhi_ ) dqmStore_->removeElement( meBCETMapProjPhi_->getName() );
00433     meBCETMapProjPhi_ = 0;
00434 
00435     if ( meBCSizMapProjEta_ ) dqmStore_->removeElement( meBCSizMapProjEta_->getName() );
00436     meBCSizMapProjEta_ = 0;
00437 
00438     if ( meBCSizMapProjPhi_ ) dqmStore_->removeElement( meBCSizMapProjPhi_->getName() );
00439     meBCSizMapProjPhi_ = 0;
00440 
00441     if ( meSCEne_ ) dqmStore_->removeElement( meSCEne_->getName() );
00442     meSCEne_ = 0;
00443 
00444     if ( meSCNum_ ) dqmStore_->removeElement( meSCNum_->getName() );
00445     meSCNum_ = 0;
00446 
00447     if ( meSCSiz_ ) dqmStore_->removeElement( meSCSiz_->getName() );
00448     meSCSiz_ = 0;
00449 
00450     if ( meSCCrystalSiz_ ) dqmStore_->removeElement( meSCCrystalSiz_->getName() );
00451     meSCCrystalSiz_ = 0;
00452 
00453     if ( meSCSeedEne_ ) dqmStore_->removeElement( meSCSeedEne_->getName() );
00454     meSCSeedEne_ = 0;
00455 
00456     if ( meSCEne2_ ) dqmStore_->removeElement( meSCEne2_->getName() );
00457     meSCEne2_ = 0;
00458 
00459     if ( meSCEneVsEMax_ ) dqmStore_->removeElement( meSCEneVsEMax_->getName() );
00460     meSCEneVsEMax_ = 0;
00461 
00462     if ( meSCEneLowScale_ ) dqmStore_->removeElement( meSCEneLowScale_->getName() );
00463     meSCEneLowScale_ = 0;
00464 
00465     if ( meSCSeedMapOcc_ ) dqmStore_->removeElement( meSCSeedMapOcc_->getName() );
00466     meSCSeedMapOcc_ = 0;
00467 
00468     if ( meSCMapSingleCrystal_ ) dqmStore_->removeElement( meSCMapSingleCrystal_->getName() );
00469     meSCMapSingleCrystal_ = 0;
00470 
00471     if ( mes1s9_ ) dqmStore_->removeElement( mes1s9_->getName() );
00472     mes1s9_ = 0;
00473 
00474     if ( mes1s9thr_ ) dqmStore_->removeElement( mes1s9thr_->getName() );
00475     mes1s9thr_ = 0;
00476 
00477     if ( mes9s25_ ) dqmStore_->removeElement( mes9s25_->getName() );
00478     mes9s25_ = 0;
00479 
00480     if ( meInvMassPi0_ ) dqmStore_->removeElement( meInvMassPi0_->getName() );
00481     meInvMassPi0_ = 0;
00482 
00483     if ( meInvMassJPsi_ ) dqmStore_->removeElement( meInvMassJPsi_->getName() );
00484     meInvMassJPsi_ = 0;
00485 
00486     if ( meInvMassZ0_ ) dqmStore_->removeElement( meInvMassZ0_->getName() );
00487     meInvMassZ0_ = 0;
00488 
00489     if ( meInvMassHigh_ ) dqmStore_->removeElement( meInvMassHigh_->getName() );
00490     meInvMassHigh_ = 0;
00491 
00492     if ( meInvMassPi0Sel_ ) dqmStore_->removeElement( meInvMassPi0Sel_->getName() );
00493     meInvMassPi0Sel_ = 0;
00494 
00495     if ( meInvMassJPsiSel_ ) dqmStore_->removeElement( meInvMassJPsiSel_->getName() );
00496     meInvMassJPsiSel_ = 0;
00497 
00498     if ( meInvMassZ0Sel_ ) dqmStore_->removeElement( meInvMassZ0Sel_->getName() );
00499     meInvMassZ0Sel_ = 0;
00500 
00501     if ( meInvMassHighSel_ ) dqmStore_->removeElement( meInvMassHighSel_->getName() );
00502     meInvMassHighSel_ = 0;
00503 
00504   }
00505 
00506   init_ = false;
00507 
00508 }
00509 
00510 void EBClusterTask::endJob(void){
00511 
00512   edm::LogInfo("EBClusterTask") << "analyzed " << ievt_ << " events";
00513 
00514   if ( enableCleanup_ ) this->cleanup();
00515 
00516 }
00517 
00518 void EBClusterTask::analyze(const edm::Event& e, const edm::EventSetup& c){
00519 
00520   bool enable = false;
00521 
00522   edm::Handle<EcalRawDataCollection> dcchs;
00523 
00524   if ( e.getByLabel(EcalRawDataCollection_, dcchs) ) {
00525 
00526     for ( EcalRawDataCollection::const_iterator dcchItr = dcchs->begin(); dcchItr != dcchs->end(); ++dcchItr ) {
00527 
00528       if ( Numbers::subDet( *dcchItr ) != EcalBarrel ) continue;
00529 
00530       if ( dcchItr->getRunType() == EcalDCCHeaderBlock::BEAMH4 ||
00531            dcchItr->getRunType() == EcalDCCHeaderBlock::BEAMH2 ) enable = true;
00532 
00533       if ( dcchItr->getRunType() == EcalDCCHeaderBlock::COSMIC ||
00534            dcchItr->getRunType() == EcalDCCHeaderBlock::MTCC ||
00535            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00536            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00537            dcchItr->getRunType() == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00538            dcchItr->getRunType() == EcalDCCHeaderBlock::PHYSICS_LOCAL ) enable = true;
00539 
00540       break;
00541 
00542     }
00543 
00544   } else {
00545 
00546     enable = true;
00547     edm::LogWarning("EBClusterTask") << EcalRawDataCollection_ << " not available";
00548 
00549   }
00550 
00551   if ( ! enable ) return;
00552 
00553   if ( ! init_ ) this->setup();
00554 
00555   ievt_++;
00556 
00557   // ECAL topology
00558   edm::ESHandle<CaloTopology> pTopology;
00559   c.get<CaloTopologyRecord>().get(pTopology);
00560   if ( !pTopology.isValid() ) {
00561     edm::LogWarning("EBClusterTask") << "Topology not valid";
00562     return;
00563   }
00564   const CaloTopology* topology = pTopology.product();
00565 
00566   // recHits
00567   edm::Handle< EcalRecHitCollection > pEBRecHits;
00568   e.getByLabel( EcalRecHitCollection_, pEBRecHits );
00569   if ( !pEBRecHits.isValid() ) {
00570     edm::LogWarning("EBClusterTask") << "RecHit collection " << EcalRecHitCollection_ << " not available.";
00571     return;
00572   }
00573   const EcalRecHitCollection* ebRecHits = pEBRecHits.product();
00574 
00575   reco::BasicClusterCollection bcSel;
00576 
00577   // --- Barrel Basic Clusters ---
00578   edm::Handle<reco::BasicClusterCollection> pBasicClusters;
00579   if ( e.getByLabel(BasicClusterCollection_, pBasicClusters) ) {
00580 
00581     int nbcc = pBasicClusters->size();
00582     if ( nbcc > 0 ) meBCNum_->Fill(float(nbcc));
00583 
00584     for ( reco::BasicClusterCollection::const_iterator bCluster = pBasicClusters->begin(); bCluster != pBasicClusters->end(); ++bCluster ) {
00585 
00586       meBCEne_->Fill(bCluster->energy());
00587       meBCSiz_->Fill(float(bCluster->size()));
00588 
00589       float xphi = bCluster->phi();
00590       if ( xphi > M_PI*(9-1.5)/9 ) xphi = xphi - M_PI*2;
00591 
00592       meBCEneMap_->Fill(xphi, bCluster->eta(), bCluster->energy());
00593       meBCEneMapProjEta_->Fill(bCluster->eta(), bCluster->energy());
00594       meBCEneMapProjPhi_->Fill(xphi, bCluster->energy());
00595 
00596       meBCNumMap_->Fill(xphi, bCluster->eta());
00597       meBCNumMapProjEta_->Fill(bCluster->eta());
00598       meBCNumMapProjPhi_->Fill(xphi);
00599 
00600       meBCSizMap_->Fill(xphi, bCluster->eta(), float(bCluster->size()));
00601       meBCSizMapProjEta_->Fill(bCluster->eta(), float(bCluster->size()));
00602       meBCSizMapProjPhi_->Fill(xphi, float(bCluster->size()));
00603 
00604       meBCETMap_->Fill(xphi, bCluster->eta(), float(bCluster->energy()) * sin(bCluster->position().theta()));
00605       meBCETMapProjEta_->Fill(bCluster->eta(), float(bCluster->energy()) * sin(bCluster->position().theta()));
00606       meBCETMapProjPhi_->Fill(xphi, float(bCluster->energy()) * sin(bCluster->position().theta()));
00607 
00608       float e2x2 = EcalClusterTools::e2x2( *bCluster, ebRecHits, topology );
00609       float e3x3 = EcalClusterTools::e3x3( *bCluster, ebRecHits, topology );
00610 
00611       // fill the selected cluster collection
00612       float pt = std::abs( bCluster->energy()*sin(bCluster->position().theta()) );
00613       if ( pt > thrClusEt_ && e2x2/e3x3 > thrS4S9_ ) bcSel.push_back(*bCluster);
00614     }
00615 
00616   } else {
00617     edm::LogWarning("EBClusterTask") << BasicClusterCollection_ << " not available";
00618   }
00619 
00620   for ( reco::BasicClusterCollection::const_iterator bc1 = bcSel.begin(); bc1 != bcSel.end(); ++bc1 ) {
00621     TLorentzVector bc1P;
00622     bc1P.SetPtEtaPhiE(std::abs(bc1->energy()*sin(bc1->position().theta())),
00623                       bc1->eta(), bc1->phi(), bc1->energy());
00624     for ( reco::BasicClusterCollection::const_iterator bc2 = bc1+1; bc2 != bcSel.end(); ++bc2 ) {
00625       TLorentzVector bc2P;
00626       bc2P.SetPtEtaPhiE(std::abs(bc2->energy()*sin(bc2->position().theta())),
00627                         bc2->eta(), bc2->phi(), bc2->energy());
00628 
00629       TLorentzVector candP = bc1P + bc2P;
00630 
00631       if ( candP.Pt() > thrCandEt_ ) {
00632         float mass = candP.M();
00633         if ( mass < 0.500 ) {
00634           meInvMassPi0Sel_->Fill( mass );
00635         } else if ( mass > 2.9 && mass < 3.3 ) {
00636           meInvMassJPsiSel_->Fill( mass );
00637         } else if ( mass > 40 && mass < 110 ) {
00638           meInvMassZ0Sel_->Fill( mass );
00639         } else if ( mass > 110 ) {
00640           meInvMassHighSel_->Fill( mass );
00641         }
00642 
00643       }
00644 
00645     }
00646   }
00647 
00648   // --- Barrel Super Clusters ---
00649   edm::Handle<reco::SuperClusterCollection> pSuperClusters;
00650   if ( e.getByLabel(SuperClusterCollection_, pSuperClusters) ) {
00651 
00652     int nscc = pSuperClusters->size();
00653     if ( nscc > 0 ) meSCNum_->Fill(float(nscc));
00654 
00655     TLorentzVector sc1_p(0,0,0,0);
00656     TLorentzVector sc2_p(0,0,0,0);
00657 
00658     for ( reco::SuperClusterCollection::const_iterator sCluster = pSuperClusters->begin(); sCluster != pSuperClusters->end(); ++sCluster ) {
00659 
00660       // energy, size
00661       meSCEne_->Fill( sCluster->energy() );
00662       meSCSiz_->Fill( float(sCluster->clustersSize()) );
00663 
00664       reco::CaloClusterPtr theSeed = sCluster->seed();
00665 
00666       // Find the seed rec hit
00667       std::vector< std::pair<DetId,float> > sIds = sCluster->hitsAndFractions();
00668 
00669       float eMax, e2nd;
00670       EcalRecHitCollection::const_iterator seedItr = ebRecHits->begin();
00671       EcalRecHitCollection::const_iterator secondItr = ebRecHits->begin();
00672 
00673       for(std::vector< std::pair<DetId,float> >::const_iterator idItr = sIds.begin(); idItr != sIds.end(); ++idItr) {
00674         DetId id = idItr->first;
00675         if(id.det() != DetId::Ecal) { continue; }
00676         EcalRecHitCollection::const_iterator hitItr = ebRecHits->find(id);
00677         if(hitItr == ebRecHits->end()) { continue; }
00678         if(hitItr->energy() > secondItr->energy()) { secondItr = hitItr; }
00679         if(hitItr->energy() > seedItr->energy()) { std::swap(seedItr,secondItr); }
00680       }
00681 
00682       eMax = seedItr->energy();
00683       e2nd = secondItr->energy();
00684       EBDetId seedId = (EBDetId) seedItr->id();
00685 
00686       float e3x3 = EcalClusterTools::e3x3( *theSeed, ebRecHits, topology );
00687       float e5x5 = EcalClusterTools::e5x5( *theSeed, ebRecHits, topology );
00688 
00689       meSCCrystalSiz_->Fill(sIds.size());
00690       meSCSeedEne_->Fill(eMax);
00691       meSCEne2_->Fill(eMax+e2nd);
00692       meSCEneVsEMax_->Fill(eMax,sCluster->energy());
00693       meSCEneLowScale_->Fill(sCluster->energy());
00694 
00695       // Prepare to fill maps
00696       int ebeta = seedId.ieta();
00697       int ebphi = seedId.iphi();
00698       float xebeta = ebeta - 0.5 * seedId.zside();
00699       float xebphi = ebphi - 0.5;
00700 
00701       meSCSeedMapOcc_->Fill(xebphi,xebeta);
00702 
00703       if(sIds.size() == 1) meSCMapSingleCrystal_->Fill(xebphi,xebeta);
00704 
00705       mes1s9_->Fill( eMax/e3x3 );
00706       if ( eMax > 3.0 ) mes1s9thr_->Fill( eMax/e3x3 );
00707       mes9s25_->Fill( e3x3/e5x5 );
00708 
00709       if ( nscc >= 2 ) {
00710         // look for the two most energetic super clusters
00711         if ( sCluster->energy() > sc1_p.Energy() ) {
00712           sc2_p=sc1_p;
00713           sc1_p.SetPtEtaPhiE(std::abs(sCluster->energy()*sin(sCluster->position().theta())),
00714                              sCluster->eta(), sCluster->phi(), sCluster->energy());
00715         } else if ( sCluster->energy() > sc2_p.Energy() ) {
00716           sc2_p.SetPtEtaPhiE(std::abs(sCluster->energy()*sin(sCluster->position().theta())),
00717                              sCluster->eta(), sCluster->phi(), sCluster->energy());
00718         }
00719       }
00720 
00721     }
00722     // Get the invariant mass of the two most energetic super clusters
00723     if ( nscc >= 2 ) {
00724       TLorentzVector sum = sc1_p+sc2_p;
00725       float mass = sum.M();
00726       if ( mass < 0.500 ) {
00727         meInvMassPi0_->Fill( mass );
00728       } else if ( mass > 2.9 && mass < 3.3 ) {
00729         meInvMassJPsi_->Fill( mass );
00730       } else if ( mass > 40 && mass < 110 ) {
00731         meInvMassZ0_->Fill( mass );
00732       } else if ( mass > 110 ) {
00733         meInvMassHigh_->Fill( mass );
00734       }
00735     }
00736 
00737   } else {
00738 
00739     edm::LogWarning("EBClusterTask") << SuperClusterCollection_ << " not available";
00740 
00741   }
00742 
00743 }
/data/refman/pasoursint/CMSSW_4_1_8_patch9/src/DQM/EcalBarrelMonitorTasks/src/EBClusterTask.cc
