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
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
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
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
1.7.3