Reference Manual

Go to the documentation of this file.
00001 #include "../interface/ClusterTask.h"
00002 
00003 #include <cmath>
00004 
00005 #include "FWCore/Framework/interface/ESHandle.h"
00006 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
00007 #include "Geometry/CaloTopology/interface/CaloTopology.h"
00008 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00009 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00010 #include "Geometry/Records/interface/CaloTopologyRecord.h"
00011 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00012 #include "DataFormats/EcalRawData/interface/EcalDCCHeaderBlock.h"
00013 #include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
00014 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
00015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00016 
00017 namespace ecaldqm {
00018 
00019   ClusterTask::ClusterTask(const edm::ParameterSet &_params, const::edm::ParameterSet& _paths) :
00020     DQWorkerTask(_params, _paths, "ClusterTask"),
00021     topology_(0),
00022     ebGeometry_(0),
00023     eeGeometry_(0),
00024     ebHits_(0),
00025     eeHits_(0),
00026     ievt_(0),
00027     lowEMax_(0.),
00028     massCalcPrescale_(0)
00029   {
00030     collectionMask_ = 
00031       (0x1 << kRun) |
00032       (0x1 << kEBRecHit) |
00033       (0x1 << kEERecHit) |
00034       (0x1 << kEBBasicCluster) |
00035       (0x1 << kEEBasicCluster) |
00036       (0x1 << kEBSuperCluster) |
00037       (0x1 << kEESuperCluster);
00038 
00039     dependencies_.push_back(std::pair<Collections, Collections>(kEBSuperCluster, kEBRecHit));
00040     dependencies_.push_back(std::pair<Collections, Collections>(kEESuperCluster, kEERecHit));
00041 
00042     edm::ParameterSet const& taskParams(_params.getUntrackedParameterSet(name_));
00043 
00044     lowEMax_ = taskParams.getUntrackedParameter<double>("lowEMax");
00045     massCalcPrescale_ = taskParams.getUntrackedParameter<int>("massCalcPrescale");
00046   }
00047 
00048   ClusterTask::~ClusterTask()
00049   {
00050   }
00051 
00052   void
00053   ClusterTask::beginRun(const edm::Run &, const edm::EventSetup &_es)
00054   {
00055     edm::ESHandle<CaloTopology> topoHndl;
00056     _es.get<CaloTopologyRecord>().get(topoHndl);
00057     topology_ = topoHndl.product();
00058     if(!topology_)
00059       throw cms::Exception("EventSetup") << "CaloTopology missing" << std::endl;
00060 
00061     edm::ESHandle<CaloGeometry> geomHndl;
00062     _es.get<CaloGeometryRecord>().get(geomHndl);
00063     ebGeometry_ = geomHndl->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
00064     eeGeometry_ = geomHndl->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
00065     if(!ebGeometry_ || !eeGeometry_)
00066       throw cms::Exception("EventSetup") << "CaloSubdetectorGeometry missing" << std::endl;
00067 
00068     ievt_ = 0;
00069   }
00070 
00071   void
00072   ClusterTask::beginEvent(const edm::Event &, const edm::EventSetup &)
00073   {
00074     ebHits_ = 0;
00075     eeHits_ = 0;
00076 
00077     ievt_++;
00078   }
00079 
00080   void
00081   ClusterTask::bookMEs()
00082   {
00083     DQWorker::bookMEs();
00084 
00085     MEs_[kBCE]->setAxisTitle("energy (GeV)", 1);
00086     MEs_[kBCNum]->setAxisTitle("number of clusters", 1);
00087     MEs_[kBCSize]->setAxisTitle("number of clusters", 1);
00088     MEs_[kSCE]->setAxisTitle("energy (GeV)", 1);
00089     MEs_[kSCELow]->setAxisTitle("energy (GeV)", 1);
00090     MEs_[kSCSeedEnergy]->setAxisTitle("energy (GeV)", 1);
00091     MEs_[kSCClusterVsSeed]->setAxisTitle("seed crystal energy (GeV)", 1);
00092     MEs_[kSCClusterVsSeed]->setAxisTitle("cluster energy (GeV)", 2);
00093     MEs_[kSCNum]->setAxisTitle("number of clusters", 1);
00094     MEs_[kSCNBCs]->setAxisTitle("cluster size", 1);
00095     MEs_[kSCNcrystals]->setAxisTitle("cluster size", 1);
00096     MEs_[kSCR9]->setAxisTitle("R9", 1);
00097     MEs_[kPi0]->setAxisTitle("mass (GeV)", 1);
00098     MEs_[kJPsi]->setAxisTitle("mass (GeV)", 1);
00099     MEs_[kZ]->setAxisTitle("mass (GeV)", 1);
00100     MEs_[kHighMass]->setAxisTitle("mass (GeV)", 1);
00101 
00102     for(int i(0); i < 2; i++)
00103       MEs_[kSCELow]->getME(i)->getTH1()->GetXaxis()->SetLimits(0., lowEMax_);
00104   }
00105   
00106   bool
00107   ClusterTask::filterRunType(const std::vector<short>& _runType)
00108   {
00109     for(int iFED(0); iFED < 54; iFED++){
00110       if ( _runType[iFED] == EcalDCCHeaderBlock::COSMIC ||
00111            _runType[iFED] == EcalDCCHeaderBlock::MTCC ||
00112            _runType[iFED] == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
00113            _runType[iFED] == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
00114            _runType[iFED] == EcalDCCHeaderBlock::COSMICS_LOCAL ||
00115            _runType[iFED] == EcalDCCHeaderBlock::PHYSICS_LOCAL ) return true;
00116     }
00117 
00118     return false;
00119   }
00120 
00121   void 
00122   ClusterTask::runOnRecHits(const EcalRecHitCollection &_hits, Collections _collection)
00123   {
00124     switch(_collection){
00125     case kEBRecHit:
00126       ebHits_ = &_hits;
00127       break;
00128     case kEERecHit:
00129       eeHits_ = &_hits;
00130       break;
00131     default:
00132       break;
00133     }
00134   }
00135 
00136   void
00137   ClusterTask::runOnBasicClusters(const reco::BasicClusterCollection &_bcs, Collections _collection)
00138   {
00139     using namespace std;
00140 
00141     int nBC[] = {0, 0};
00142     bool isBarrel(_collection == kEBBasicCluster);
00143 
00144     vector<reco::BasicCluster const*> lowMassCands;
00145 
00146     for(reco::BasicClusterCollection::const_iterator bcItr(_bcs.begin()); bcItr != _bcs.end(); ++bcItr){
00147       const math::XYZPoint &position(bcItr->position());
00148 
00149       DetId id(bcItr->seed());
00150       if(id.null()){
00151         GlobalPoint gp(position.x(), position.y(), position.z());
00152         const CaloSubdetectorGeometry* subgeom(isBarrel ? ebGeometry_ : eeGeometry_);
00153 
00154         id = subgeom->getClosestCell(gp);
00155       }
00156 
00157       if(id.null() || (id.subdetId() == EcalBarrel && !isBarrel) || (id.subdetId() == EcalEndcap && isBarrel)) continue;
00158 
00159       float energy(bcItr->energy());
00160 
00161       MEs_[kBCE]->fill(id, energy);
00162 
00163       MEs_[kBCEMap]->fill(id, energy);
00164       MEs_[kBCEMapProjEta]->fill(id, energy);
00165       MEs_[kBCEMapProjPhi]->fill(id, energy);
00166 
00167       MEs_[kBCOccupancy]->fill(id);
00168       MEs_[kBCOccupancyProjEta]->fill(id);
00169       MEs_[kBCOccupancyProjPhi]->fill(id);
00170 
00171       float size(bcItr->size());
00172 
00173       MEs_[kBCSize]->fill(id, size);
00174 
00175       MEs_[kBCSizeMap]->fill(id, size);
00176       MEs_[kBCSizeMapProjEta]->fill(id, size);
00177       MEs_[kBCSizeMapProjPhi]->fill(id, size);
00178 
00179       int zside(position.z() > 0 ? 1 : 0);
00180       nBC[zside]++;
00181 
00182       if(ievt_ % massCalcPrescale_ != 0) continue;
00183 
00184       if(energy > 10.) continue;
00185 
00186       EcalRecHitCollection::const_iterator hitItr(isBarrel ? ebHits_->find(id) : eeHits_->find(id));
00187       if(hitItr == (isBarrel ? ebHits_->end() : eeHits_->end())) continue;
00188 
00189       // cuts here must be parametrized
00190       if(hitItr->energy() < 0.5) continue;
00191 
00192       if(hitItr->energy() / energy > 0.95) continue;
00193 
00194       lowMassCands.push_back(&(*bcItr));
00195     }
00196 
00197     if(isBarrel){
00198       MEs_[kBCNum]->fill((unsigned)BinService::kEB + 1, nBC[0] + nBC[1]);
00199     }else{
00200       MEs_[kBCNum]->fill((unsigned)BinService::kEEm + 1, nBC[0]);
00201       MEs_[kBCNum]->fill((unsigned)BinService::kEEp + 1, nBC[1]);
00202     }
00203 
00204     if(ievt_ % massCalcPrescale_ != 0) return;
00205 
00206     for(vector<reco::BasicCluster const*>::iterator bcItr1(lowMassCands.begin()); bcItr1 != lowMassCands.end(); ++bcItr1){
00207       reco::BasicCluster const& bc1(**bcItr1);
00208       float energy1(bc1.energy());
00209       float px1(energy1 * sin(bc1.position().theta()) * cos(bc1.phi()));
00210       float py1(energy1 * sin(bc1.position().theta()) * sin(bc1.phi()));
00211       float pz1(energy1 * cos(bc1.position().theta()));
00212 
00213       for(vector<reco::BasicCluster const*>::iterator bcItr2(lowMassCands.begin()); bcItr2 != lowMassCands.end(); ++bcItr2){
00214         if(*bcItr1 == *bcItr2) continue;
00215         reco::BasicCluster const& bc2(**bcItr2);
00216         float energy2(bc2.energy());
00217         float px2(energy2 * sin(bc2.position().theta()) * cos(bc2.phi()));
00218         float py2(energy2 * sin(bc2.position().theta()) * sin(bc2.phi()));
00219         float pz2(energy2 * cos(bc2.position().theta()));
00220 
00221         float ptpair(sqrt((px1 + px2) * (px1 + px2) + (py1 + py2) * (py1 + py2)));
00222         if(ptpair < 2.5) continue;
00223 
00224         float epair(energy1 + energy2);
00225         float pzpair(abs(pz1 + pz2));
00226         
00227         if(epair < pzpair + 1.e-10) continue;
00228         float eta(0.5 * log((epair + pzpair)/(epair - pzpair)));
00229         float phi(atan2(px1 + px2, py1 + py2));
00230 
00231         float iso(0.);
00232         for(reco::BasicClusterCollection::const_iterator bcItr(_bcs.begin()); bcItr != _bcs.end(); ++bcItr){
00233           float dEta(bcItr->eta() - eta);
00234           float dPhi(bcItr->phi() - phi);
00235           if(sqrt(dEta * dEta + dPhi * dPhi) < 0.2) iso += bcItr->energy() * sin(bcItr->position().theta());
00236         }
00237         if(iso > 0.5) continue;
00238 
00239         float mass(sqrt(epair * epair - pzpair * pzpair - ptpair * ptpair));
00240         MEs_[kPi0]->fill(mass);
00241         MEs_[kJPsi]->fill(mass);
00242       }
00243     }
00244   }
00245 
00246   void
00247   ClusterTask::runOnSuperClusters(const reco::SuperClusterCollection &_scs, Collections _collection)
00248   {
00249     using namespace std;
00250 
00251     const EcalRecHitCollection *hits(0);
00252     bool isBarrel;
00253     if(_collection == kEBSuperCluster){
00254       hits = ebHits_;
00255       isBarrel = true;
00256     }else{
00257       hits = eeHits_;
00258       isBarrel = false;
00259     }
00260 
00261     reco::SuperCluster const* leading(0);
00262     reco::SuperCluster const* subLeading(0);
00263 
00264     int nSC[] = {0, 0};
00265 
00266     for(reco::SuperClusterCollection::const_iterator scItr(_scs.begin()); scItr != _scs.end(); ++scItr){
00267       const math::XYZPoint &position(scItr->position());
00268 
00269       DetId id(scItr->seed()->seed());
00270       if(id.null()){
00271         GlobalPoint gp(position.x(), position.y(), position.z());
00272         const CaloSubdetectorGeometry* subgeom(isBarrel ? ebGeometry_ : eeGeometry_);
00273 
00274         id = subgeom->getClosestCell(gp);
00275       }
00276 
00277       if(id.null() || (id.subdetId() == EcalBarrel && !isBarrel) || (id.subdetId() == EcalEndcap && isBarrel)) continue;
00278 
00279       float energy(scItr->energy());
00280 
00281       MEs_[kSCE]->fill(id, energy);
00282       if(energy < lowEMax_) MEs_[kSCELow]->fill(id, energy);
00283 
00284       MEs_[kSCNBCs]->fill(id, scItr->clustersSize());
00285       MEs_[kSCNcrystals]->fill(id, scItr->size());
00286 
00287       if(!hits) continue;
00288       EcalRecHitCollection::const_iterator seedItr(hits->find(id));
00289       if(seedItr == hits->end()) continue;
00290 
00291       MEs_[kSCSeedEnergy]->fill(id, seedItr->energy());
00292       MEs_[kSCClusterVsSeed]->fill(id, seedItr->energy(), energy);
00293 
00294       MEs_[kSCSeedOccupancy]->fill(id);
00295 
00296       if(_scs.size() == 1)
00297         MEs_[kSingleCrystalCluster]->fill(id);
00298 
00299       float e3x3(EcalClusterTools::e3x3(*scItr->seed(), hits, topology_));
00300       MEs_[kSCR9]->fill(id, e3x3 / energy);
00301 
00302       int zside(position.z() > 0 ? 1 : 0);
00303       nSC[zside]++;
00304 
00305       if(ievt_ % massCalcPrescale_ != 0) continue;
00306 
00307       float et(energy * sin(scItr->position().theta()));
00308       if(!leading || et > leading->energy() * sin(leading->position().theta())){
00309         subLeading = leading;
00310         leading = &(*scItr);
00311       }
00312       else if(!subLeading || et > subLeading->energy() * sin(subLeading->position().theta())){
00313         subLeading = &(*scItr);
00314       }
00315     }
00316 
00317     if(_collection == kEBSuperCluster){
00318       MEs_[kSCNum]->fill((unsigned)BinService::kEB + 1, nSC[0] + nSC[1]);
00319     }else{
00320       MEs_[kSCNum]->fill((unsigned)BinService::kEEm + 1, nSC[0]);
00321       MEs_[kSCNum]->fill((unsigned)BinService::kEEp + 1, nSC[1]);
00322     }
00323 
00324     if(ievt_ % massCalcPrescale_ != 0) return;
00325 
00326     // implement isolation & cuts
00327     if(!leading || !subLeading) return;
00328     float e(leading->energy() + subLeading->energy());
00329     float px(leading->energy() * sin(leading->position().theta()) * cos(leading->phi()) + subLeading->energy() * sin(subLeading->position().theta()) * cos(subLeading->phi()));
00330     float py(leading->energy() * sin(leading->position().theta()) * sin(leading->phi()) + subLeading->energy() * sin(subLeading->position().theta()) * sin(subLeading->phi()));
00331     float pz(leading->energy() * cos(leading->position().theta()) + subLeading->energy() * cos(subLeading->position().theta()));
00332     float mass(sqrt(e * e - px * px - py * py - pz * pz));
00333     MEs_[kZ]->fill(mass);
00334     MEs_[kHighMass]->fill(mass);
00335 
00336   }
00337 
00338   /*static*/
00339   void
00340   ClusterTask::setMEData(std::vector<MEData>& _data)
00341   {
00342     BinService::AxisSpecs xaxis, yaxis, zaxis;
00343 
00344     zaxis.low = 0.;
00345     zaxis.high = 50.;
00346     _data[kBCEMap] = MEData("BCEMap", BinService::kEcal3P, BinService::kSuperCrystal, MonitorElement::DQM_KIND_TPROFILE2D, 0, 0, &zaxis);
00347     _data[kBCEMapProjEta] = MEData("BCEMap", BinService::kEcal3P, BinService::kProjEta, MonitorElement::DQM_KIND_TPROFILE);
00348     _data[kBCEMapProjPhi] = MEData("BCEMap", BinService::kEcal3P, BinService::kProjPhi, MonitorElement::DQM_KIND_TPROFILE);
00349 
00350     _data[kBCOccupancy] = MEData("BCOccupancy", BinService::kEcal3P, BinService::kSuperCrystal, MonitorElement::DQM_KIND_TH2F);
00351     _data[kBCOccupancyProjEta] = MEData("BCOccupancy", BinService::kEcal3P, BinService::kProjEta, MonitorElement::DQM_KIND_TH1F);
00352     _data[kBCOccupancyProjPhi] = MEData("BCOccupancy", BinService::kEcal3P, BinService::kProjPhi, MonitorElement::DQM_KIND_TH1F);
00353 
00354     zaxis.high = 30.;
00355     _data[kBCSizeMap] = MEData("BCSizeMap", BinService::kEcal3P, BinService::kSuperCrystal, MonitorElement::DQM_KIND_TPROFILE2D, 0, 0, &zaxis);
00356     _data[kBCSizeMapProjEta] = MEData("BCSizeMap", BinService::kEcal3P, BinService::kProjEta, MonitorElement::DQM_KIND_TPROFILE);
00357     _data[kBCSizeMapProjPhi] = MEData("BCSizeMap", BinService::kEcal3P, BinService::kProjPhi, MonitorElement::DQM_KIND_TPROFILE);
00358 
00359     xaxis.nbins = 50;
00360     xaxis.low = 0.;
00361     xaxis.high = 150.;
00362     _data[kBCE] = MEData("BCE", BinService::kEcal3P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00363 
00364     xaxis.nbins = 20;
00365     xaxis.low = 0.;
00366     xaxis.high = 100.;
00367     _data[kBCNum] = MEData("BCNum", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00368 
00369     xaxis.nbins = 50;
00370     xaxis.low = 0.;
00371     xaxis.high = 100.;
00372     _data[kBCSize] = MEData("BCSize", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00373 
00374     xaxis.nbins = 50;
00375     xaxis.low = 0.;
00376     xaxis.high = 150.;
00377     _data[kSCE] = MEData("SCE", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00378 
00379     xaxis.nbins = 50;
00380     xaxis.low = 0.;
00381     xaxis.high = 10.;
00382     _data[kSCELow] = MEData("SCELow", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00383 
00384     xaxis.nbins = 50;
00385     xaxis.low = 0.;
00386     xaxis.high = 150.;
00387     _data[kSCSeedEnergy] = MEData("SCSeedEnergy", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00388 
00389     yaxis.nbins = 50;
00390     yaxis.low = 0.;
00391     yaxis.high = 150.;
00392     _data[kSCClusterVsSeed] = MEData("SCClusterVsSeed", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH2F, &xaxis, &yaxis);
00393 
00394     _data[kSCSeedOccupancy] = MEData("SCSeedOccupancy", BinService::kEcal3P, BinService::kSuperCrystal, MonitorElement::DQM_KIND_TH2F);
00395     _data[kSingleCrystalCluster] = MEData("SCSingleCrystalCluster", BinService::kEcal3P, BinService::kSuperCrystal, MonitorElement::DQM_KIND_TH2F);
00396 
00397     xaxis.nbins = 20;
00398     xaxis.low = 0.;
00399     xaxis.high = 20.;
00400     _data[kSCNum] = MEData("SCNum", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00401 
00402     xaxis.nbins = 15;
00403     xaxis.low = 0.;
00404     xaxis.high = 15.;
00405     _data[kSCNBCs] = MEData("SCNBCs", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00406     xaxis.nbins = 50;
00407     xaxis.low = 0.;
00408     xaxis.high = 150.;
00409     _data[kSCNcrystals] = MEData("SCNcrystals", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00410 
00411     xaxis.nbins = 50;
00412     xaxis.low = 0.;
00413     xaxis.high = 1.2;
00414     _data[kSCR9] = MEData("SCR9", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00415 
00416 
00417     xaxis.nbins = 50;
00418     xaxis.low = 0.;
00419     xaxis.high = 0.5;
00420     _data[kPi0] = MEData("Pi0", BinService::nObjType, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00421     xaxis.low = 2.9;
00422     xaxis.high = 3.3;
00423     _data[kJPsi] = MEData("JPsi", BinService::nObjType, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00424     xaxis.low = 40.;
00425     xaxis.high = 110.;
00426     _data[kZ] = MEData("Z", BinService::nObjType, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00427     xaxis.low = 110.;
00428     xaxis.high = 3000.;
00429     _data[kHighMass] = MEData("HighMass", BinService::nObjType, BinService::kUser, MonitorElement::DQM_KIND_TH1F, &xaxis);
00430   }
00431 
00432   DEFINE_ECALDQM_WORKER(ClusterTask);
00433 }
00434
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/DQM/EcalBarrelMonitorTasks/src/ClusterTask.cc
