CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ClusterTask.cc
Go to the documentation of this file.
1 #include "../interface/ClusterTask.h"
2 
3 #include <cmath>
4 
16 
17 namespace ecaldqm {
18 
19  ClusterTask::ClusterTask(const edm::ParameterSet &_params, const::edm::ParameterSet& _paths) :
20  DQWorkerTask(_params, _paths, "ClusterTask"),
21  topology_(0),
22  ebGeometry_(0),
23  eeGeometry_(0),
24  ebHits_(0),
25  eeHits_(0),
26  ievt_(0),
27  lowEMax_(0.),
28  massCalcPrescale_(0)
29  {
31  (0x1 << kRun) |
32  (0x1 << kEBRecHit) |
33  (0x1 << kEERecHit) |
34  (0x1 << kEBBasicCluster) |
35  (0x1 << kEEBasicCluster) |
36  (0x1 << kEBSuperCluster) |
37  (0x1 << kEESuperCluster);
38 
39  dependencies_.push_back(std::pair<Collections, Collections>(kEBSuperCluster, kEBRecHit));
40  dependencies_.push_back(std::pair<Collections, Collections>(kEESuperCluster, kEERecHit));
41 
42  edm::ParameterSet const& taskParams(_params.getUntrackedParameterSet(name_));
43 
44  lowEMax_ = taskParams.getUntrackedParameter<double>("lowEMax");
45  massCalcPrescale_ = taskParams.getUntrackedParameter<int>("massCalcPrescale");
46  }
47 
49  {
50  }
51 
52  void
54  {
56  _es.get<CaloTopologyRecord>().get(topoHndl);
57  topology_ = topoHndl.product();
58  if(!topology_)
59  throw cms::Exception("EventSetup") << "CaloTopology missing" << std::endl;
60 
62  _es.get<CaloGeometryRecord>().get(geomHndl);
63  ebGeometry_ = geomHndl->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
64  eeGeometry_ = geomHndl->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
65  if(!ebGeometry_ || !eeGeometry_)
66  throw cms::Exception("EventSetup") << "CaloSubdetectorGeometry missing" << std::endl;
67 
68  ievt_ = 0;
69  }
70 
71  void
73  {
74  ebHits_ = 0;
75  eeHits_ = 0;
76 
77  ievt_++;
78  }
79 
80  void
82  {
84 
85  MEs_[kBCE]->setAxisTitle("energy (GeV)", 1);
86  MEs_[kBCNum]->setAxisTitle("number of clusters", 1);
87  MEs_[kBCSize]->setAxisTitle("number of clusters", 1);
88  MEs_[kSCE]->setAxisTitle("energy (GeV)", 1);
89  MEs_[kSCELow]->setAxisTitle("energy (GeV)", 1);
90  MEs_[kSCSeedEnergy]->setAxisTitle("energy (GeV)", 1);
91  MEs_[kSCClusterVsSeed]->setAxisTitle("seed crystal energy (GeV)", 1);
92  MEs_[kSCClusterVsSeed]->setAxisTitle("cluster energy (GeV)", 2);
93  MEs_[kSCNum]->setAxisTitle("number of clusters", 1);
94  MEs_[kSCNBCs]->setAxisTitle("cluster size", 1);
95  MEs_[kSCNcrystals]->setAxisTitle("cluster size", 1);
96  MEs_[kSCR9]->setAxisTitle("R9", 1);
97  MEs_[kPi0]->setAxisTitle("mass (GeV)", 1);
98  MEs_[kJPsi]->setAxisTitle("mass (GeV)", 1);
99  MEs_[kZ]->setAxisTitle("mass (GeV)", 1);
100  MEs_[kHighMass]->setAxisTitle("mass (GeV)", 1);
101 
102  for(int i(0); i < 2; i++)
103  MEs_[kSCELow]->getME(i)->getTH1()->GetXaxis()->SetLimits(0., lowEMax_);
104  }
105 
106  bool
107  ClusterTask::filterRunType(const std::vector<short>& _runType)
108  {
109  for(int iFED(0); iFED < 54; iFED++){
110  if ( _runType[iFED] == EcalDCCHeaderBlock::COSMIC ||
111  _runType[iFED] == EcalDCCHeaderBlock::MTCC ||
112  _runType[iFED] == EcalDCCHeaderBlock::COSMICS_GLOBAL ||
113  _runType[iFED] == EcalDCCHeaderBlock::PHYSICS_GLOBAL ||
114  _runType[iFED] == EcalDCCHeaderBlock::COSMICS_LOCAL ||
115  _runType[iFED] == EcalDCCHeaderBlock::PHYSICS_LOCAL ) return true;
116  }
117 
118  return false;
119  }
120 
121  void
123  {
124  switch(_collection){
125  case kEBRecHit:
126  ebHits_ = &_hits;
127  break;
128  case kEERecHit:
129  eeHits_ = &_hits;
130  break;
131  default:
132  break;
133  }
134  }
135 
136  void
138  {
139  using namespace std;
140 
141  int nBC[] = {0, 0};
142  bool isBarrel(_collection == kEBBasicCluster);
143 
144  vector<reco::BasicCluster const*> lowMassCands;
145 
146  for(reco::BasicClusterCollection::const_iterator bcItr(_bcs.begin()); bcItr != _bcs.end(); ++bcItr){
147  const math::XYZPoint &position(bcItr->position());
148 
149  DetId id(bcItr->seed());
150  if(id.null()){
151  GlobalPoint gp(position.x(), position.y(), position.z());
152  const CaloSubdetectorGeometry* subgeom(isBarrel ? ebGeometry_ : eeGeometry_);
153 
154  id = subgeom->getClosestCell(gp);
155  }
156 
157  if(id.null() || (id.subdetId() == EcalBarrel && !isBarrel) || (id.subdetId() == EcalEndcap && isBarrel)) continue;
158 
159  float energy(bcItr->energy());
160 
161  MEs_[kBCE]->fill(id, energy);
162 
163  MEs_[kBCEMap]->fill(id, energy);
164  MEs_[kBCEMapProjEta]->fill(id, energy);
165  MEs_[kBCEMapProjPhi]->fill(id, energy);
166 
167  MEs_[kBCOccupancy]->fill(id);
168  MEs_[kBCOccupancyProjEta]->fill(id);
169  MEs_[kBCOccupancyProjPhi]->fill(id);
170 
171  float size(bcItr->size());
172 
173  MEs_[kBCSize]->fill(id, size);
174 
175  MEs_[kBCSizeMap]->fill(id, size);
176  MEs_[kBCSizeMapProjEta]->fill(id, size);
177  MEs_[kBCSizeMapProjPhi]->fill(id, size);
178 
179  int zside(position.z() > 0 ? 1 : 0);
180  nBC[zside]++;
181 
182  if(ievt_ % massCalcPrescale_ != 0) continue;
183 
184  if(energy > 10.) continue;
185 
186  EcalRecHitCollection::const_iterator hitItr(isBarrel ? ebHits_->find(id) : eeHits_->find(id));
187  if(hitItr == (isBarrel ? ebHits_->end() : eeHits_->end())) continue;
188 
189  // cuts here must be parametrized
190  if(hitItr->energy() < 0.5) continue;
191 
192  if(hitItr->energy() / energy > 0.95) continue;
193 
194  lowMassCands.push_back(&(*bcItr));
195  }
196 
197  if(isBarrel){
198  MEs_[kBCNum]->fill((unsigned)BinService::kEB + 1, nBC[0] + nBC[1]);
199  }else{
200  MEs_[kBCNum]->fill((unsigned)BinService::kEEm + 1, nBC[0]);
201  MEs_[kBCNum]->fill((unsigned)BinService::kEEp + 1, nBC[1]);
202  }
203 
204  if(ievt_ % massCalcPrescale_ != 0) return;
205 
206  for(vector<reco::BasicCluster const*>::iterator bcItr1(lowMassCands.begin()); bcItr1 != lowMassCands.end(); ++bcItr1){
207  reco::BasicCluster const& bc1(**bcItr1);
208  float energy1(bc1.energy());
209  float px1(energy1 * sin(bc1.position().theta()) * cos(bc1.phi()));
210  float py1(energy1 * sin(bc1.position().theta()) * sin(bc1.phi()));
211  float pz1(energy1 * cos(bc1.position().theta()));
212 
213  for(vector<reco::BasicCluster const*>::iterator bcItr2(lowMassCands.begin()); bcItr2 != lowMassCands.end(); ++bcItr2){
214  if(*bcItr1 == *bcItr2) continue;
215  reco::BasicCluster const& bc2(**bcItr2);
216  float energy2(bc2.energy());
217  float px2(energy2 * sin(bc2.position().theta()) * cos(bc2.phi()));
218  float py2(energy2 * sin(bc2.position().theta()) * sin(bc2.phi()));
219  float pz2(energy2 * cos(bc2.position().theta()));
220 
221  float ptpair(sqrt((px1 + px2) * (px1 + px2) + (py1 + py2) * (py1 + py2)));
222  if(ptpair < 2.5) continue;
223 
224  float epair(energy1 + energy2);
225  float pzpair(abs(pz1 + pz2));
226 
227  if(epair < pzpair + 1.e-10) continue;
228  float eta(0.5 * log((epair + pzpair)/(epair - pzpair)));
229  float phi(atan2(px1 + px2, py1 + py2));
230 
231  float iso(0.);
232  for(reco::BasicClusterCollection::const_iterator bcItr(_bcs.begin()); bcItr != _bcs.end(); ++bcItr){
233  float dEta(bcItr->eta() - eta);
234  float dPhi(bcItr->phi() - phi);
235  if(sqrt(dEta * dEta + dPhi * dPhi) < 0.2) iso += bcItr->energy() * sin(bcItr->position().theta());
236  }
237  if(iso > 0.5) continue;
238 
239  float mass(sqrt(epair * epair - pzpair * pzpair - ptpair * ptpair));
240  MEs_[kPi0]->fill(mass);
241  MEs_[kJPsi]->fill(mass);
242  }
243  }
244  }
245 
246  void
248  {
249  using namespace std;
250 
251  const EcalRecHitCollection *hits(0);
252  bool isBarrel;
253  if(_collection == kEBSuperCluster){
254  hits = ebHits_;
255  isBarrel = true;
256  }else{
257  hits = eeHits_;
258  isBarrel = false;
259  }
260 
261  reco::SuperCluster const* leading(0);
262  reco::SuperCluster const* subLeading(0);
263 
264  int nSC[] = {0, 0};
265 
266  for(reco::SuperClusterCollection::const_iterator scItr(_scs.begin()); scItr != _scs.end(); ++scItr){
267  const math::XYZPoint &position(scItr->position());
268 
269  DetId id(scItr->seed()->seed());
270  if(id.null()){
271  GlobalPoint gp(position.x(), position.y(), position.z());
272  const CaloSubdetectorGeometry* subgeom(isBarrel ? ebGeometry_ : eeGeometry_);
273 
274  id = subgeom->getClosestCell(gp);
275  }
276 
277  if(id.null() || (id.subdetId() == EcalBarrel && !isBarrel) || (id.subdetId() == EcalEndcap && isBarrel)) continue;
278 
279  float energy(scItr->energy());
280 
281  MEs_[kSCE]->fill(id, energy);
282  if(energy < lowEMax_) MEs_[kSCELow]->fill(id, energy);
283 
284  MEs_[kSCNBCs]->fill(id, scItr->clustersSize());
285  MEs_[kSCNcrystals]->fill(id, scItr->size());
286 
287  if(!hits) continue;
288  EcalRecHitCollection::const_iterator seedItr(hits->find(id));
289  if(seedItr == hits->end()) continue;
290 
291  MEs_[kSCSeedEnergy]->fill(id, seedItr->energy());
292  MEs_[kSCClusterVsSeed]->fill(id, seedItr->energy(), energy);
293 
294  MEs_[kSCSeedOccupancy]->fill(id);
295 
296  if(_scs.size() == 1)
298 
299  float e3x3(EcalClusterTools::e3x3(*scItr->seed(), hits, topology_));
300  MEs_[kSCR9]->fill(id, e3x3 / energy);
301 
302  int zside(position.z() > 0 ? 1 : 0);
303  nSC[zside]++;
304 
305  if(ievt_ % massCalcPrescale_ != 0) continue;
306 
307  float et(energy * sin(scItr->position().theta()));
308  if(!leading || et > leading->energy() * sin(leading->position().theta())){
309  subLeading = leading;
310  leading = &(*scItr);
311  }
312  else if(!subLeading || et > subLeading->energy() * sin(subLeading->position().theta())){
313  subLeading = &(*scItr);
314  }
315  }
316 
317  if(_collection == kEBSuperCluster){
318  MEs_[kSCNum]->fill((unsigned)BinService::kEB + 1, nSC[0] + nSC[1]);
319  }else{
320  MEs_[kSCNum]->fill((unsigned)BinService::kEEm + 1, nSC[0]);
321  MEs_[kSCNum]->fill((unsigned)BinService::kEEp + 1, nSC[1]);
322  }
323 
324  if(ievt_ % massCalcPrescale_ != 0) return;
325 
326  // implement isolation & cuts
327  if(!leading || !subLeading) return;
328  float e(leading->energy() + subLeading->energy());
329  float px(leading->energy() * sin(leading->position().theta()) * cos(leading->phi()) + subLeading->energy() * sin(subLeading->position().theta()) * cos(subLeading->phi()));
330  float py(leading->energy() * sin(leading->position().theta()) * sin(leading->phi()) + subLeading->energy() * sin(subLeading->position().theta()) * sin(subLeading->phi()));
331  float pz(leading->energy() * cos(leading->position().theta()) + subLeading->energy() * cos(subLeading->position().theta()));
332  float mass(sqrt(e * e - px * px - py * py - pz * pz));
333  MEs_[kZ]->fill(mass);
334  MEs_[kHighMass]->fill(mass);
335 
336  }
337 
338  /*static*/
339  void
340  ClusterTask::setMEData(std::vector<MEData>& _data)
341  {
342  BinService::AxisSpecs xaxis, yaxis, zaxis;
343 
344  zaxis.low = 0.;
345  zaxis.high = 50.;
349 
353 
354  zaxis.high = 30.;
358 
359  xaxis.nbins = 50;
360  xaxis.low = 0.;
361  xaxis.high = 150.;
363 
364  xaxis.nbins = 20;
365  xaxis.low = 0.;
366  xaxis.high = 100.;
368 
369  xaxis.nbins = 50;
370  xaxis.low = 0.;
371  xaxis.high = 100.;
373 
374  xaxis.nbins = 50;
375  xaxis.low = 0.;
376  xaxis.high = 150.;
378 
379  xaxis.nbins = 50;
380  xaxis.low = 0.;
381  xaxis.high = 10.;
383 
384  xaxis.nbins = 50;
385  xaxis.low = 0.;
386  xaxis.high = 150.;
388 
389  yaxis.nbins = 50;
390  yaxis.low = 0.;
391  yaxis.high = 150.;
392  _data[kSCClusterVsSeed] = MEData("SCClusterVsSeed", BinService::kEcal2P, BinService::kUser, MonitorElement::DQM_KIND_TH2F, &xaxis, &yaxis);
393 
396 
397  xaxis.nbins = 20;
398  xaxis.low = 0.;
399  xaxis.high = 20.;
401 
402  xaxis.nbins = 15;
403  xaxis.low = 0.;
404  xaxis.high = 15.;
406  xaxis.nbins = 50;
407  xaxis.low = 0.;
408  xaxis.high = 150.;
410 
411  xaxis.nbins = 50;
412  xaxis.low = 0.;
413  xaxis.high = 1.2;
415 
416 
417  xaxis.nbins = 50;
418  xaxis.low = 0.;
419  xaxis.high = 0.5;
421  xaxis.low = 2.9;
422  xaxis.high = 3.3;
424  xaxis.low = 40.;
425  xaxis.high = 110.;
427  xaxis.low = 110.;
428  xaxis.high = 3000.;
430  }
431 
433 }
434 
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:123
int i
Definition: DBlmapReader.cc:9
std::vector< std::pair< Collections, Collections > > dependencies_
Definition: DQWorkerTask.h:31
string fill
Definition: lumiContext.py:319
const EcalRecHitCollection * eeHits_
Definition: ClusterTask.h:70
bool filterRunType(const std::vector< short > &)
Definition: ClusterTask.cc:107
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
std::vector< EcalRecHit >::const_iterator const_iterator
#define abs(x)
Definition: mlp_lapack.h:159
T eta() const
static int position[TOTALCHAMBERS][3]
Definition: ReadPGInfo.cc:509
const CaloSubdetectorGeometry * eeGeometry_
Definition: ClusterTask.h:69
static void setMEData(std::vector< MEData > &)
Definition: ClusterTask.cc:340
void beginEvent(const edm::Event &, const edm::EventSetup &)
Definition: ClusterTask.cc:72
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
const EcalRecHitCollection * ebHits_
Definition: ClusterTask.h:70
std::vector< SuperCluster > SuperClusterCollection
collection of SuperCluser objectr
ParameterSet const & getUntrackedParameterSet(std::string const &name, ParameterSet const &defaultValue) const
void beginRun(const edm::Run &, const edm::EventSetup &)
Definition: ClusterTask.cc:53
T sqrt(T t)
Definition: SSEVec.h:48
void runOnRecHits(const EcalRecHitCollection &, Collections)
Definition: ClusterTask.cc:122
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
static float e3x3(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology)
double energy() const
cluster energy
Definition: CaloCluster.h:120
void runOnSuperClusters(const reco::SuperClusterCollection &, Collections)
Definition: ClusterTask.cc:247
const_iterator end() const
Definition: DetId.h:20
virtual void bookMEs()
Definition: DQWorker.cc:48
std::vector< MESet * > MEs_
Definition: DQWorker.h:56
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:13
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
std::vector< BasicCluster > BasicClusterCollection
collection of BasicCluster objects
const CaloTopology * topology_
Definition: ClusterTask.h:67
iterator find(key_type k)
ClusterTask(const edm::ParameterSet &, const edm::ParameterSet &)
Definition: ClusterTask.cc:19
void runOnBasicClusters(const reco::BasicClusterCollection &, Collections)
Definition: ClusterTask.cc:137
DEFINE_ECALDQM_WORKER(CertificationClient)
double phi() const
azimuthal angle of cluster centroid
Definition: CaloCluster.h:163
tuple size
Write out results.
std::string name_
Definition: DQWorker.h:55
Definition: Run.h:36
Definition: DDAxes.h:10
const CaloSubdetectorGeometry * ebGeometry_
Definition: ClusterTask.h:68