CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CaloTowersCreationAlgo.cc
Go to the documentation of this file.
9 #include "Math/Interpolator.h"
10 #include<iostream>
11 
13  : theEBthreshold(-1000.),
14  theEEthreshold(-1000.),
15 
16  theUseEtEBTresholdFlag(false),
17  theUseEtEETresholdFlag(false),
18  theUseSymEBTresholdFlag(false),
19  theUseSymEETresholdFlag(false),
20 
21 
22  theHcalThreshold(-1000.),
23  theHBthreshold(-1000.),
24  theHESthreshold(-1000.),
25  theHEDthreshold(-1000.),
26  theHOthreshold0(-1000.),
27  theHOthresholdPlus1(-1000.),
28  theHOthresholdMinus1(-1000.),
29  theHOthresholdPlus2(-1000.),
30  theHOthresholdMinus2(-1000.),
31  theHF1threshold(-1000.),
32  theHF2threshold(-1000.),
33  theEBGrid(std::vector<double>(5,10.)),
34  theEBWeights(std::vector<double>(5,1.)),
35  theEEGrid(std::vector<double>(5,10.)),
36  theEEWeights(std::vector<double>(5,1.)),
37  theHBGrid(std::vector<double>(5,10.)),
38  theHBWeights(std::vector<double>(5,1.)),
39  theHESGrid(std::vector<double>(5,10.)),
40  theHESWeights(std::vector<double>(5,1.)),
41  theHEDGrid(std::vector<double>(5,10.)),
42  theHEDWeights(std::vector<double>(5,1.)),
43  theHOGrid(std::vector<double>(5,10.)),
44  theHOWeights(std::vector<double>(5,1.)),
45  theHF1Grid(std::vector<double>(5,10.)),
46  theHF1Weights(std::vector<double>(5,1.)),
47  theHF2Grid(std::vector<double>(5,10.)),
48  theHF2Weights(std::vector<double>(5,1.)),
49  theEBweight(1.),
50  theEEweight(1.),
51  theHBweight(1.),
52  theHESweight(1.),
53  theHEDweight(1.),
54  theHOweight(1.),
55  theHF1weight(1.),
56  theHF2weight(1.),
57  theEcutTower(-1000.),
58  theEBSumThreshold(-1000.),
59  theEESumThreshold(-1000.),
60  theHcalTopology(0),
61  theGeometry(0),
62  theTowerConstituentsMap(0),
63  theHOIsUsed(true),
64  // (for momentum reconstruction algorithm)
65  theMomConstrMethod(0),
66  theMomHBDepth(0.),
67  theMomHEDepth(0.),
68  theMomEBDepth(0.),
69  theMomEEDepth(0.)
70 {
71 }
72 
73 CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold, double EEthreshold,
74 
75  bool useEtEBTreshold,
76  bool useEtEETreshold,
77  bool useSymEBTreshold,
78  bool useSymEETreshold,
79 
80  double HcalThreshold,
81  double HBthreshold, double HESthreshold, double HEDthreshold,
82  double HOthreshold0, double HOthresholdPlus1, double HOthresholdMinus1,
83  double HOthresholdPlus2, double HOthresholdMinus2,
84  double HF1threshold, double HF2threshold,
85  double EBweight, double EEweight,
86  double HBweight, double HESweight, double HEDweight,
87  double HOweight, double HF1weight, double HF2weight,
88  double EcutTower, double EBSumThreshold, double EESumThreshold,
89  bool useHO,
90  // (momentum reconstruction algorithm)
91  int momConstrMethod,
92  double momHBDepth,
93  double momHEDepth,
94  double momEBDepth,
95  double momEEDepth)
96 
97  : theEBthreshold(EBthreshold),
98  theEEthreshold(EEthreshold),
99 
100  theUseEtEBTresholdFlag(useEtEBTreshold),
101  theUseEtEETresholdFlag(useEtEETreshold),
102  theUseSymEBTresholdFlag(useSymEBTreshold),
103  theUseSymEETresholdFlag(useSymEETreshold),
104 
105  theHcalThreshold(HcalThreshold),
106  theHBthreshold(HBthreshold),
107  theHESthreshold(HESthreshold),
108  theHEDthreshold(HEDthreshold),
109  theHOthreshold0(HOthreshold0),
110  theHOthresholdPlus1(HOthresholdPlus1),
111  theHOthresholdMinus1(HOthresholdMinus1),
112  theHOthresholdPlus2(HOthresholdPlus2),
113  theHOthresholdMinus2(HOthresholdMinus2),
114  theHF1threshold(HF1threshold),
115  theHF2threshold(HF2threshold),
116  theEBGrid(std::vector<double>(5,10.)),
117  theEBWeights(std::vector<double>(5,1.)),
118  theEEGrid(std::vector<double>(5,10.)),
119  theEEWeights(std::vector<double>(5,1.)),
120  theHBGrid(std::vector<double>(5,10.)),
121  theHBWeights(std::vector<double>(5,1.)),
122  theHESGrid(std::vector<double>(5,10.)),
123  theHESWeights(std::vector<double>(5,1.)),
124  theHEDGrid(std::vector<double>(5,10.)),
125  theHEDWeights(std::vector<double>(5,1.)),
126  theHOGrid(std::vector<double>(5,10.)),
127  theHOWeights(std::vector<double>(5,1.)),
128  theHF1Grid(std::vector<double>(5,10.)),
129  theHF1Weights(std::vector<double>(5,1.)),
130  theHF2Grid(std::vector<double>(5,10.)),
131  theHF2Weights(std::vector<double>(5,1.)),
132  theEBweight(EBweight),
133  theEEweight(EEweight),
134  theHBweight(HBweight),
135  theHESweight(HESweight),
136  theHEDweight(HEDweight),
137  theHOweight(HOweight),
138  theHF1weight(HF1weight),
139  theHF2weight(HF2weight),
140  theEcutTower(EcutTower),
141  theEBSumThreshold(EBSumThreshold),
142  theEESumThreshold(EESumThreshold),
143  theHOIsUsed(useHO),
144  // (momentum reconstruction algorithm)
145  theMomConstrMethod(momConstrMethod),
146  theMomHBDepth(momHBDepth),
147  theMomHEDepth(momHEDepth),
148  theMomEBDepth(momEBDepth),
149  theMomEEDepth(momEEDepth)
150 
151 {
152 }
153 
154 CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold, double EEthreshold,
155 
156  bool useEtEBTreshold,
157  bool useEtEETreshold,
158  bool useSymEBTreshold,
159  bool useSymEETreshold,
160 
161  double HcalThreshold,
162  double HBthreshold, double HESthreshold, double HEDthreshold,
163  double HOthreshold0, double HOthresholdPlus1, double HOthresholdMinus1,
164  double HOthresholdPlus2, double HOthresholdMinus2,
165  double HF1threshold, double HF2threshold,
166  std::vector<double> EBGrid, std::vector<double> EBWeights,
167  std::vector<double> EEGrid, std::vector<double> EEWeights,
168  std::vector<double> HBGrid, std::vector<double> HBWeights,
169  std::vector<double> HESGrid, std::vector<double> HESWeights,
170  std::vector<double> HEDGrid, std::vector<double> HEDWeights,
171  std::vector<double> HOGrid, std::vector<double> HOWeights,
172  std::vector<double> HF1Grid, std::vector<double> HF1Weights,
173  std::vector<double> HF2Grid, std::vector<double> HF2Weights,
174  double EBweight, double EEweight,
175  double HBweight, double HESweight, double HEDweight,
176  double HOweight, double HF1weight, double HF2weight,
177  double EcutTower, double EBSumThreshold, double EESumThreshold,
178  bool useHO,
179  // (for the momentum construction algorithm)
180  int momConstrMethod,
181  double momHBDepth,
182  double momHEDepth,
183  double momEBDepth,
184  double momEEDepth)
185 
186  : theEBthreshold(EBthreshold),
187  theEEthreshold(EEthreshold),
188 
189  theUseEtEBTresholdFlag(useEtEBTreshold),
190  theUseEtEETresholdFlag(useEtEETreshold),
191  theUseSymEBTresholdFlag(useSymEBTreshold),
192  theUseSymEETresholdFlag(useSymEETreshold),
193 
194  theHcalThreshold(HcalThreshold),
195  theHBthreshold(HBthreshold),
196  theHESthreshold(HESthreshold),
197  theHEDthreshold(HEDthreshold),
198  theHOthreshold0(HOthreshold0),
199  theHOthresholdPlus1(HOthresholdPlus1),
200  theHOthresholdMinus1(HOthresholdMinus1),
201  theHOthresholdPlus2(HOthresholdPlus2),
202  theHOthresholdMinus2(HOthresholdMinus2),
203  theHF1threshold(HF1threshold),
204  theHF2threshold(HF2threshold),
205  theEBGrid(EBGrid),
206  theEBWeights(EBWeights),
207  theEEGrid(EEGrid),
208  theEEWeights(EEWeights),
209  theHBGrid(HBGrid),
210  theHBWeights(HBWeights),
211  theHESGrid(HESGrid),
212  theHESWeights(HESWeights),
213  theHEDGrid(HEDGrid),
214  theHEDWeights(HEDWeights),
215  theHOGrid(HOGrid),
216  theHOWeights(HOWeights),
217  theHF1Grid(HF1Grid),
218  theHF1Weights(HF1Weights),
219  theHF2Grid(HF2Grid),
220  theHF2Weights(HF2Weights),
221  theEBweight(EBweight),
222  theEEweight(EEweight),
223  theHBweight(HBweight),
224  theHESweight(HESweight),
225  theHEDweight(HEDweight),
226  theHOweight(HOweight),
227  theHF1weight(HF1weight),
228  theHF2weight(HF2weight),
229  theEcutTower(EcutTower),
230  theEBSumThreshold(EBSumThreshold),
231  theEESumThreshold(EESumThreshold),
232  theHOIsUsed(useHO),
233  // (momentum reconstruction algorithm)
234  theMomConstrMethod(momConstrMethod),
235  theMomHBDepth(momHBDepth),
236  theMomHEDepth(momHEDepth),
237  theMomEBDepth(momEBDepth),
238  theMomEEDepth(momEEDepth)
239 
240 {
241 }
242 
243 
246  theHcalTopology = topo;
247  theGeometry = geo;
249 }
250 
252  theTowerMap.clear();
253  hcalDropChMap.clear();
254 }
255 
257  for(HBHERecHitCollection::const_iterator hbheItr = hbhe.begin();
258  hbheItr != hbhe.end(); ++hbheItr)
259  assignHit(&(*hbheItr));
260 }
261 
263  for(HORecHitCollection::const_iterator hoItr = ho.begin();
264  hoItr != ho.end(); ++hoItr)
265  assignHit(&(*hoItr));
266 }
267 
269  for(HFRecHitCollection::const_iterator hfItr = hf.begin();
270  hfItr != hf.end(); ++hfItr)
271  assignHit(&(*hfItr));
272 }
273 
276  ecItr != ec.end(); ++ecItr)
277  assignHit(&(*ecItr));
278 }
279 
280 // this method should not be used any more as the towers in the changed format
281 // can not be properly rescaled with the "rescale" method.
282 // "rescale was replaced by "rescaleTowers"
283 //
285  for(CaloTowerCollection::const_iterator ctcItr = ctc.begin();
286  ctcItr != ctc.end(); ++ctcItr) {
287  rescale(&(*ctcItr));
288  }
289 }
290 
291 
292 
294  // now copy this map into the final collection
295  for(MetaTowerMap::const_iterator mapItr = theTowerMap.begin();
296  mapItr != theTowerMap.end(); ++mapItr) {
297 
298  // Convert only if there is at least one constituent in the metatower.
299  // The check of constituents size in the coverted tower is still needed!
300  if ( (mapItr->second).metaConstituents.size()<1) continue;
301 
302  CaloTower ct=convert(mapItr->first,mapItr->second);
303  if (ct.constituentsSize()>0 && ct.energy()>theEcutTower) {
304  result.push_back(ct);
305  }
306  }
307  theTowerMap.clear(); // save the memory
308 }
309 
310 
312 
313  for (CaloTowerCollection::const_iterator ctcItr = ctc.begin();
314  ctcItr != ctc.end(); ++ctcItr) {
315 
316  CaloTowerDetId twrId = ctcItr->id();
317  double newE_em = ctcItr->emEnergy();
318  double newE_had = ctcItr->hadEnergy();
319  double newE_outer = ctcItr->outerEnergy();
320 
321  double threshold = 0.0; // not used: we do not change thresholds
322  double weight = 1.0;
323 
324  // HF
325  if (ctcItr->ietaAbs()>=30) {
326  double E_short = 0.5 * newE_had; // from the definitions for HF
327  double E_long = newE_em + 0.5 * newE_had; //
328  // scale
329  E_long *= theHF1weight;
330  E_short *= theHF2weight;
331  // convert
332  newE_em = E_long - E_short;
333  newE_had = 2.0 * E_short;
334  }
335 
336  else { // barrel/endcap
337 
338  // find if its in EB, or EE; determine from first ecal constituent found
339  for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
340  DetId constId = ctcItr->constituent(iConst);
341  if (constId.det()!=DetId::Ecal) continue;
342  getThresholdAndWeight(constId, threshold, weight);
343  newE_em *= weight;
344  break;
345  }
346  // HO
347  for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
348  DetId constId = ctcItr->constituent(iConst);
349  if (constId.det()!=DetId::Hcal) continue;
350  if (HcalDetId(constId).subdet()!=HcalOuter) continue;
351  getThresholdAndWeight(constId, threshold, weight);
352  newE_outer *= weight;
353  break;
354  }
355  // HB/HE
356  for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
357  DetId constId = ctcItr->constituent(iConst);
358  if (constId.det()!=DetId::Hcal) continue;
359  if (HcalDetId(constId).subdet()==HcalOuter) continue;
360  getThresholdAndWeight(constId, threshold, weight);
361  newE_had *= weight;
362  if (ctcItr->ietaAbs()>16) newE_outer *= weight;
363  break;
364  }
365 
366  } // barrel/endcap region
367 
368  // now make the new tower
369 
370  double newE_hadTot = (theHOIsUsed && twrId.ietaAbs()<16)? newE_had+newE_outer : newE_had;
371 
372  GlobalPoint emPoint = ctcItr->emPosition();
373  GlobalPoint hadPoint = ctcItr->emPosition();
374 
375  double f_em = 1.0/cosh(emPoint.eta());
376  double f_had = 1.0/cosh(hadPoint.eta());
377 
379 
380  if (ctcItr->ietaAbs()<30) {
381  if (newE_em>0) towerP4 += CaloTower::PolarLorentzVector(newE_em*f_em, emPoint.eta(), emPoint.phi(), 0);
382  if (newE_hadTot>0) towerP4 += CaloTower::PolarLorentzVector(newE_hadTot*f_had, hadPoint.eta(), hadPoint.phi(), 0);
383  }
384  else {
385  double newE_tot = newE_em + newE_had;
386  // for HF we use common point for ecal, hcal shower positions regardless of the method
387  if (newE_tot>0) towerP4 += CaloTower::PolarLorentzVector(newE_tot*f_had, hadPoint.eta(), hadPoint.phi(), 0);
388  }
389 
390 
391 
392  CaloTower rescaledTower(twrId, newE_em, newE_had, newE_outer, -1, -1, towerP4, emPoint, hadPoint);
393  // copy the timings, have to convert back to int, 1 unit = 0.01 ns
394  rescaledTower.setEcalTime( int(ctcItr->ecalTime()*100.0 + 0.5) );
395  rescaledTower.setHcalTime( int(ctcItr->hcalTime()*100.0 + 0.5) );
396 
397  std::vector<DetId> contains;
398  for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
399  contains.push_back(ctcItr->constituent(iConst));
400  }
401  rescaledTower.addConstituents(contains);
402 
403  rescaledTower.setCaloTowerStatus(ctcItr->towerStatusWord());
404 
405  ctcResult.push_back(rescaledTower);
406 
407  } // end of loop over towers
408 
409 
410 }
411 
412 
413 
415  DetId detId = recHit->detid();
416 
417  unsigned int chStatusForCT = (detId.det()==DetId::Hcal)?
420 
421  // this is for skipping channls: mostly needed for the creation of
422  // bad towers from hits i the bad channel collections.
423  if (chStatusForCT==CaloTowersCreationAlgo::IgnoredChan) return;
424 
425  double threshold, weight;
426  getThresholdAndWeight(detId, threshold, weight);
427 
428  double energy = recHit->energy(); // original RecHit energy is used to apply thresholds
429  double e = energy * weight; // energies scaled by user weight: used in energy assignments
430 
431 
432  // SPECIAL handling of tower 28/depth 3 --> half into tower 28 and half into tower 29
433  if (detId.det()==DetId::Hcal &&
434  HcalDetId(detId).subdet()==HcalEndcap &&
435  HcalDetId(detId).depth()==3 &&
436  HcalDetId(detId).ietaAbs()==28) {
437 
438  CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
439  if (towerDetId.null()) return;
440  MetaTower & tower28 = find(towerDetId);
441  CaloTowerDetId towerDetId29 = CaloTowerDetId(towerDetId.ieta()+
442  towerDetId.zside(),
443  towerDetId.iphi());
444  MetaTower & tower29 = find(towerDetId29);
445 
447 
448  // bad channels are counted regardless of energy threshold
449 
450  if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
451  tower28.numBadHcalCells += 1;
452  tower29.numBadHcalCells += 1;
453  }
454 
455  else if (0.5*energy >= threshold) { // not bad channel: use energy if above threshold
456 
457  if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
458  tower28.numRecHcalCells += 1;
459  tower29.numRecHcalCells += 1;
460  }
461  else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
462  tower28.numProbHcalCells += 1;
463  tower29.numProbHcalCells += 1;
464  }
465 
466  // NOTE DIVIDE BY 2!!!
467  double e28 = 0.5 * e;
468  double e29 = 0.5 * e;
469 
470  tower28.E_had += e28;
471  tower28.E += e28;
472  std::pair<DetId,double> mc(detId,e28);
473  tower28.metaConstituents.push_back(mc);
474 
475  tower29.E_had += e29;
476  tower29.E += e29;
477  tower29.metaConstituents.push_back(mc);
478 
479  // time info: do not use in averaging if timing error is found: need
480  // full set of status info to implement: use only "good" channels for now
481 
482  if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
483  tower28.hadSumTimeTimesE += ( e28 * recHit->time() );
484  tower28.hadSumEForTime += e28;
485  tower29.hadSumTimeTimesE += ( e29 * recHit->time() );
486  tower29.hadSumEForTime += e29;
487  }
488 
489  // store the energy in layer 3 also in E_outer
490  tower28.E_outer += e28;
491  tower29.E_outer += e29;
492  } // not a "bad" hit
493 
494  } // end of special case
495 
496  else {
497 
498  CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
499 
500  if (towerDetId.null()) return;
501  MetaTower & tower = find(towerDetId);
502 
503  DetId::Detector det = detId.det();
504 
505  if (det == DetId::Ecal) {
506 
508 
509  // For ECAL we count all bad channels after the metatower is complete
510 
511  // Include options for symmetric thresholds and cut on Et
512  // for ECAL RecHits
513 
514  bool passEmThreshold = false;
515 
516  if (detId.subdetId() == EcalBarrel) {
517  if (theUseEtEBTresholdFlag) energy /= cosh( (theGeometry->getGeometry(detId)->getPosition()).eta() ) ;
518  if (theUseSymEBTresholdFlag) passEmThreshold = (fabs(energy) >= threshold);
519  else passEmThreshold = (energy >= threshold);
520 
521  }
522  else if (detId.subdetId() == EcalEndcap) {
523  if (theUseEtEETresholdFlag) energy /= cosh( (theGeometry->getGeometry(detId)->getPosition()).eta() ) ;
524  if (theUseSymEETresholdFlag) passEmThreshold = (fabs(energy) >= threshold);
525  else passEmThreshold = (energy >= threshold);
526  }
527 
528 
529  // if (chStatusForCT != CaloTowersCreationAlgo::BadChan && energy >= threshold) {
530  if (chStatusForCT != CaloTowersCreationAlgo::BadChan && passEmThreshold) {
531  tower.E_em += e;
532  tower.E += e;
533 
534  if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
535  tower.numRecEcalCells += 1;
536  }
537  else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
538  tower.numProbEcalCells += 1;
539  }
540 
541  // change when full status info is available
542  // for now use only good channels
543 
544  // add e>0 check (new options allow e<0)
545  if (chStatusForCT == CaloTowersCreationAlgo::GoodChan && e>0 ) {
546  tower.emSumTimeTimesE += ( e * recHit->time() );
547  tower.emSumEForTime += e; // see above
548  }
549 
550  std::pair<DetId,double> mc(detId,e);
551  tower.metaConstituents.push_back(mc);
552  }
553 
554  } // end of ECAL
555 
556  // HCAL
557  else {
558  HcalDetId hcalDetId(detId);
559 
561 
562  if(hcalDetId.subdet() == HcalOuter) {
563 
564  if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
565  if (theHOIsUsed) tower.numBadHcalCells += 1;
566  }
567 
568  else if (energy >= threshold) {
569  tower.E_outer += e; // store HO energy even if HO is not used
570  // add energy of the tower and/or flag if theHOIsUsed
571  if(theHOIsUsed) {
572  tower.E += e;
573 
574  if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
575  tower.numRecHcalCells += 1;
576  }
577  else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
578  tower.numProbHcalCells += 1;
579  }
580  } // HO is used
581 
582 
583  // add HO to constituents even if it is not used: JetMET wants to keep these towers
584  std::pair<DetId,double> mc(detId,e);
585  tower.metaConstituents.push_back(mc);
586 
587  } // not a bad channel, energy above threshold
588 
589  } // HO hit
590 
591  // HF calculates EM fraction differently
592  else if(hcalDetId.subdet() == HcalForward) {
593 
594  if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
595  tower.numBadHcalCells += 1;
596  }
597 
598  else if (energy >= threshold) {
599  if (hcalDetId.depth() == 1) {
600  // long fiber, so E_EM = E(Long) - E(Short)
601  tower.E_em += e;
602  }
603  else {
604  // short fiber, EHAD = 2 * E(Short)
605  tower.E_em -= e;
606  tower.E_had += 2. * e;
607  }
608  tower.E += e;
609  if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
610  tower.numRecHcalCells += 1;
611  }
612  else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
613  tower.numProbHcalCells += 1;
614  }
615 
616  // put the timing in HCAL -> have to check timing errors when available
617  // for now use only good channels
618  if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
619  tower.hadSumTimeTimesE += ( e * recHit->time() );
620  tower.hadSumEForTime += e;
621  }
622 
623  std::pair<DetId,double> mc(detId,e);
624  tower.metaConstituents.push_back(mc);
625 
626  } // not a bad HF channel, energy above threshold
627 
628  } // HF hit
629 
630  else {
631  // HCAL situation normal in HB/HE
632 
633  if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
634  tower.numBadHcalCells += 1;
635  }
636  else if (energy >= threshold) {
637  tower.E_had += e;
638  tower.E += e;
639  if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
640  tower.numRecHcalCells += 1;
641  }
642  else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
643  tower.numProbHcalCells += 1;
644  }
645 
646  // Timing information: need specific accessors
647  // for now use only good channels
648  if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
649  tower.hadSumTimeTimesE += ( e * recHit->time() );
650  tower.hadSumEForTime += e;
651  }
652  // store energy in highest depth for towers 18-27 (for electron,photon ID in endcap)
653  // also, store energy in HE part of tower 16 (for JetMET cleanup)
654  if (HcalDetId(detId).subdet()==HcalEndcap) {
655  if ( (HcalDetId(detId).depth()==2 && HcalDetId(detId).ietaAbs()>=18 && HcalDetId(detId).ietaAbs()<27) ||
656  (HcalDetId(detId).depth()==3 && HcalDetId(detId).ietaAbs()==27) ||
657  (HcalDetId(detId).depth()==3 && HcalDetId(detId).ietaAbs()==16) ) {
658  tower.E_outer += e;
659  }
660  }
661 
662  std::pair<DetId,double> mc(detId,e);
663  tower.metaConstituents.push_back(mc);
664 
665  } // not a "bad" channel, energy above threshold
666 
667  } // channel in HBHE (excluding twrs 28,29)
668 
669  }
670 
671  } // recHit normal case (not in HE towers 28,29)
672 
673 } // end of assignHit method
674 
675 
676 
677 
678 // This method is not flexible enough for the new CaloTower format.
679 // For now make a quick compatibility "fix" : WILL NOT WORK CORRECTLY with anything
680 // except the default simple p4 assignment!!!
681 // Must be rewritten for full functionality.
683  double threshold, weight;
684 
685  CaloTowerDetId towerDetId = ct->id();
686 // if (towerDetId.null()) return;
687  MetaTower & tower = find(towerDetId);
688 
689  tower.E_em = 0.;
690  tower.E_had = 0.;
691  tower.E_outer = 0.;
692  for (unsigned int i=0; i<ct->constituentsSize(); i++) {
693  DetId detId = ct->constituent(i);
694  getThresholdAndWeight(detId, threshold, weight);
695  DetId::Detector det = detId.det();
696  if(det == DetId::Ecal) {
697  tower.E_em = ct->emEnergy()*weight;
698  }
699  else {
700  HcalDetId hcalDetId(detId);
701  if(hcalDetId.subdet() == HcalForward) {
702  if (hcalDetId.depth()==1) tower.E_em = ct->emEnergy()*weight;
703  if (hcalDetId.depth()==2) tower.E_had = ct->hadEnergy()*weight;
704  }
705  else if(hcalDetId.subdet() == HcalOuter) {
706  tower.E_outer = ct->outerEnergy()*weight;
707  }
708  else {
709  tower.E_had = ct->hadEnergy()*weight;
710  }
711  }
712  tower.E = tower.E_had+tower.E_em+tower.E_outer;
713 
714  // this is to be compliant with the new MetaTower setup
715  // used only for the default simple vector assignment
716  std::pair<DetId, double> mc(detId, 0);
717  tower.metaConstituents.push_back(mc);
718  }
719 
720  // preserve time inforamtion
721  tower.emSumTimeTimesE = ct->ecalTime();
722  tower.hadSumTimeTimesE = ct->hcalTime();
723  tower.emSumEForTime = 1.0;
724  tower.hadSumEForTime = 1.0;
725 }
726 
728  E(0),E_em(0),E_had(0),E_outer(0), emSumTimeTimesE(0), hadSumTimeTimesE(0), emSumEForTime(0), hadSumEForTime(0),
729  numBadEcalCells(0), numRecEcalCells(0), numProbEcalCells(0), numBadHcalCells(0), numRecHcalCells(0), numProbHcalCells(0) { }
730 
731 
733  MetaTowerMap::iterator itr = theTowerMap.find(detId);
734  if(itr == theTowerMap.end()) {
735  // need to build a new tower
736  MetaTower t;
737 
738  // store it in the map
739  theTowerMap.insert(std::pair<CaloTowerDetId, CaloTowersCreationAlgo::MetaTower>(detId, t));
740  itr = theTowerMap.find(detId);
741 
742 
743  }
744  return itr->second;
745 }
746 
748 
749  double ecalThres=(id.ietaAbs()<=17)?(theEBSumThreshold):(theEESumThreshold);
750  double E=mt.E;
751  double E_em=mt.E_em;
752  double E_had=mt.E_had;
753  double E_outer=mt.E_outer;
754 
755  // conditional assignment of depths for barrel/endcap
756  // Some additional tuning may be required in the transitional region
757  // 14<|iEta|<19
758  if (id.ietaAbs()<=17) {
761  }
762  else {
765  }
766 
767 
768  // Note: E_outer is used to save HO energy OR energy in the outermost depths in endcap region
769  // In the methods with separate treatment of EM and HAD components:
770  // - HO is not used to determine direction, however HO energy is added to get "total had energy"
771  // => Check if the tower is within HO coverage before adding E_outer to the "total had" energy
772  // else the energy will be double counted
773  // When summing up the energy of the tower these checks are performed in the loops over RecHits
774 
775 
776  float ecalTime = (mt.emSumEForTime>0)? mt.emSumTimeTimesE/mt.emSumEForTime : -9999;
777  float hcalTime = (mt.hadSumEForTime>0)? mt.hadSumTimeTimesE/mt.hadSumEForTime : -9999;
778 
779  std::vector<std::pair<DetId,double> > metaContains=mt.metaConstituents;
780 
781  if (id.ietaAbs()<=29 && E_em<ecalThres) { // ignore EM threshold in HF
782  E-=E_em;
783  E_em=0;
784  std::vector<std::pair<DetId,double> > metaContains_noecal;
785 
786  for (std::vector<std::pair<DetId,double> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i)
787  if (i->first.det()!=DetId::Ecal) metaContains_noecal.push_back(*i);
788  metaContains.swap(metaContains_noecal);
789  }
790 
791  if (id.ietaAbs()<=29 && E_had<theHcalThreshold) {
792  E-=E_had;
793 
794  if (theHOIsUsed && id.ietaAbs()<16) E-=E_outer; // not subtracted before, think it should be done
795 
796  E_had=0;
797  E_outer=0;
798  std::vector<std::pair<DetId,double> > metaContains_nohcal;
799 
800  for (std::vector<std::pair<DetId,double> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i)
801  if (i->first.det()!=DetId::Hcal) metaContains_nohcal.push_back(*i);
802  metaContains.swap(metaContains_nohcal);
803  }
804 
805 
806  double E_had_tot = (theHOIsUsed && id.ietaAbs()<16)? E_had+E_outer : E_had;
807 
808 
809  // create CaloTower using the selected algorithm
810 
811  GlobalPoint emPoint, hadPoint;
812 
814 
815 
816  switch (theMomConstrMethod) {
817 
818  case 0 :
819  { // Simple 4-momentum assignment
821 
822  double pf=1.0/cosh(p.eta());
823  if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);
824 
825  emPoint = p;
826  hadPoint = p;
827  } // end case 0
828  break;
829 
830  case 1 :
831  { // separate 4-vectors for ECAL, HCAL, add to get the 4-vector of the tower (=>tower has mass!)
832  if (id.ietaAbs()<=29) {
833  if (E_em>0) {
834  emPoint = emShwrPos(metaContains, theMomEmDepth, E_em);
835  double emPf = 1.0/cosh(emPoint.eta());
836  towerP4 += CaloTower::PolarLorentzVector(E_em*emPf, emPoint.eta(), emPoint.phi(), 0);
837  }
838  if ( (E_had + E_outer) >0) {
839  hadPoint = hadShwrPos(id, theMomHadDepth);
840  double hadPf = 1.0/cosh(hadPoint.eta());
841 
842  if (E_had_tot>0) {
843  towerP4 += CaloTower::PolarLorentzVector(E_had_tot*hadPf, hadPoint.eta(), hadPoint.phi(), 0);
844  }
845  }
846  }
847  else { // forward detector: use the CaloTower position
849  double pf=1.0/cosh(p.eta());
850  if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0); // simple momentum assignment, same position
851  emPoint = p;
852  hadPoint = p;
853  }
854  } // end case 1
855  break;
856 
857  case 2:
858  { // use ECAL position for the tower (when E_cal>0), else default CaloTower position (massless tower)
859  if (id.ietaAbs()<=29) {
860  if (E_em>0) emPoint = emShwrLogWeightPos(metaContains, theMomEmDepth, E_em);
861  else emPoint = theTowerGeometry->getGeometry(id)->getPosition();
862 
863  double sumPf = 1.0/cosh(emPoint.eta());
864  if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*sumPf, emPoint.eta(), emPoint.phi(), 0);
865 
866  hadPoint = emPoint;
867  }
868  else { // forward detector: use the CaloTower position
870  double pf=1.0/cosh(p.eta());
871  if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0); // simple momentum assignment, same position
872  emPoint = p;
873  hadPoint = p;
874  }
875  } // end case 2
876  break;
877 
878  } // end of decision on p4 reconstruction method
879 
880 
881  CaloTower retval(id, E_em, E_had, E_outer, -1, -1, towerP4, emPoint, hadPoint);
882 
883  // set the timings
884  retval.setEcalTime(compactTime(ecalTime));
885  retval.setHcalTime(compactTime(hcalTime));
886 
887  // set the CaloTower status word =====================================
888  // Channels must be counter exclusively in the defined cathegories
889  // "Bad" channels (not used in energy assignment) can be flagged during
890  // CaloTower creation only if specified in the configuration file
891 
892  unsigned int numBadHcalChan = mt.numBadHcalCells;
893  // unsigned int numBadEcalChan = mt.numBadEcalCells;
894  unsigned int numBadEcalChan = 0; //
895 
896  unsigned int numRecHcalChan = mt.numRecHcalCells;
897  unsigned int numRecEcalChan = mt.numRecEcalCells;
898  unsigned int numProbHcalChan = mt.numProbHcalCells;
899  unsigned int numProbEcalChan = mt.numProbEcalCells;
900 
901  // now add dead/off/... channels not used in RecHit reconstruction for HCAL
902  if (hcalDropChMap.find(id) != hcalDropChMap.end()) numBadHcalChan += hcalDropChMap[id];
903 
904 
905  // for ECAL the number of all bad channels is obtained here -----------------------
906 
907  // get all possible constituents of the tower
908  std::vector<DetId> allConstituents = theTowerConstituentsMap->constituentsOf(id);
909 
910  for (std::vector<DetId>::iterator ac_it=allConstituents.begin();
911  ac_it!=allConstituents.end(); ++ac_it) {
912 
913  if (ac_it->det()!=DetId::Ecal) continue;
914 
915  int thisEcalSevLvl = -999;
916 
917  if (ac_it->subdetId() == EcalBarrel && theEbHandle.isValid()) {
918  thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEbHandle, *theEcalChStatus);
919  }
920  else if (ac_it->subdetId() == EcalEndcap && theEeHandle.isValid()) {
921  thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEeHandle, *theEcalChStatus);
922  }
923 
924  // "bad" include "recovered" if the flag not to use them was set
925  // this is consistent with the treatment of what channels are accepted
926  if (thisEcalSevLvl>int(theEcalAcceptSeverityLevel) ||
928 
929  ++numBadEcalChan;
930  }
931 
932  }
933  //--------------------------------------------------------------------------------------
934 
935  retval.setCaloTowerStatus(numBadHcalChan, numBadEcalChan,
936  numRecHcalChan, numRecEcalChan,
937  numProbHcalChan, numProbEcalChan);
938 
939  double maxCellE = -999.0; // for storing the hottest cell E in the calotower
940 
941  std::vector<DetId> contains;
942  for (std::vector<std::pair<DetId,double> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i) {
943 
944  contains.push_back(i->first);
945 
946  if (maxCellE < i->second) {
947  // need an extra check because of the funny towers that are empty except for the presence of an HO
948  // hit in the constituents (JetMET wanted them saved)
949  // This constituent is only used for storing the tower, but should not be concidered as a hot cell canditate for
950  // configurations with useHO = false
951 
952 
953  if (i->first.det()==DetId::Ecal) { // ECAL
954  maxCellE = i->second;
955  }
956  else { // HCAL
957  if (HcalDetId(i->first).subdet() != HcalOuter)
958  maxCellE = i->second;
959  else if (theHOIsUsed) maxCellE = i->second;
960  }
961 
962  } // found higher E cell
963 
964  } // loop over matacontains
965 
966  retval.addConstituents(contains);
967  retval.setHottestCellE(maxCellE);
968 
969  return retval;
970 }
971 
972 
973 
974 void CaloTowersCreationAlgo::getThresholdAndWeight(const DetId & detId, double & threshold, double & weight) const {
975  DetId::Detector det = detId.det();
976  weight=0; // in case the hit is not identified
977 
978  if(det == DetId::Ecal) {
979  // may or may not be EB. We'll find out.
980 
981  EcalSubdetector subdet = (EcalSubdetector)(detId.subdetId());
982  if(subdet == EcalBarrel) {
983  threshold = theEBthreshold;
984  weight = theEBweight;
985  if (weight <= 0.) {
986  ROOT::Math::Interpolator my(theEBGrid,theEBWeights,ROOT::Math::Interpolation::kAKIMA);
987  weight = my.Eval(theEBEScale);
988  }
989  }
990  else if(subdet == EcalEndcap) {
991  threshold = theEEthreshold;
992  weight = theEEweight;
993  if (weight <= 0.) {
994  ROOT::Math::Interpolator my(theEEGrid,theEEWeights,ROOT::Math::Interpolation::kAKIMA);
995  weight = my.Eval(theEEEScale);
996  }
997  }
998  }
999  else if(det == DetId::Hcal) {
1000  HcalDetId hcalDetId(detId);
1001  HcalSubdetector subdet = hcalDetId.subdet();
1002 
1003  if(subdet == HcalBarrel) {
1004  threshold = theHBthreshold;
1005  weight = theHBweight;
1006  if (weight <= 0.) {
1007  ROOT::Math::Interpolator my(theHBGrid,theHBWeights,ROOT::Math::Interpolation::kAKIMA);
1008  weight = my.Eval(theHBEScale);
1009  }
1010  }
1011 
1012  else if(subdet == HcalEndcap) {
1013  // check if it's single or double tower
1014  if(hcalDetId.ietaAbs() < theHcalTopology->firstHEDoublePhiRing()) {
1015  threshold = theHESthreshold;
1016  weight = theHESweight;
1017  if (weight <= 0.) {
1018  ROOT::Math::Interpolator my(theHESGrid,theHESWeights,ROOT::Math::Interpolation::kAKIMA);
1019  weight = my.Eval(theHESEScale);
1020  }
1021  }
1022  else {
1023  threshold = theHEDthreshold;
1024  weight = theHEDweight;
1025  if (weight <= 0.) {
1026  ROOT::Math::Interpolator my(theHEDGrid,theHEDWeights,ROOT::Math::Interpolation::kAKIMA);
1027  weight = my.Eval(theHEDEScale);
1028  }
1029  }
1030  }
1031 
1032  else if(subdet == HcalOuter) {
1033  //check if it's ring 0 or +1 or +2 or -1 or -2
1034  if(hcalDetId.ietaAbs() <= 4) threshold = theHOthreshold0;
1035  else if(hcalDetId.ieta() < 0) {
1036  // set threshold for ring -1 or -2
1037  threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdMinus1 : theHOthresholdMinus2;
1038  } else {
1039  // set threshold for ring +1 or +2
1040  threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdPlus1 : theHOthresholdPlus2;
1041  }
1042  weight = theHOweight;
1043  if (weight <= 0.) {
1044  ROOT::Math::Interpolator my(theHOGrid,theHOWeights,ROOT::Math::Interpolation::kAKIMA);
1045  weight = my.Eval(theHOEScale);
1046  }
1047  }
1048 
1049  else if(subdet == HcalForward) {
1050  if(hcalDetId.depth() == 1) {
1051  threshold = theHF1threshold;
1052  weight = theHF1weight;
1053  if (weight <= 0.) {
1054  ROOT::Math::Interpolator my(theHF1Grid,theHF1Weights,ROOT::Math::Interpolation::kAKIMA);
1055  weight = my.Eval(theHF1EScale);
1056  }
1057  } else {
1058  threshold = theHF2threshold;
1059  weight = theHF2weight;
1060  if (weight <= 0.) {
1061  ROOT::Math::Interpolator my(theHF2Grid,theHF2Weights,ROOT::Math::Interpolation::kAKIMA);
1062  weight = my.Eval(theHF2EScale);
1063  }
1064  }
1065  }
1066  }
1067  else {
1068  std::cout << "BAD CELL det " << det << std::endl;
1069  }
1070 }
1071 
1073  if (scale>0.00001) *&theEBEScale = scale;
1074  else *&theEBEScale = 50.;
1075 }
1076 
1078  if (scale>0.00001) *&theEEEScale = scale;
1079  else *&theEEEScale = 50.;
1080 }
1081 
1083  if (scale>0.00001) *&theHBEScale = scale;
1084  else *&theHBEScale = 50.;
1085 }
1086 
1088  if (scale>0.00001) *&theHESEScale = scale;
1089  else *&theHESEScale = 50.;
1090 }
1091 
1093  if (scale>0.00001) *&theHEDEScale = scale;
1094  else *&theHEDEScale = 50.;
1095 }
1096 
1098  if (scale>0.00001) *&theHOEScale = scale;
1099  else *&theHOEScale = 50.;
1100 }
1101 
1103  if (scale>0.00001) *&theHF1EScale = scale;
1104  else *&theHF1EScale = 50.;
1105 }
1106 
1108  if (scale>0.00001) *&theHF2EScale = scale;
1109  else *&theHF2EScale = 50.;
1110 }
1111 
1112 
1114  const CaloCellGeometry* cellGeometry = theGeometry->getGeometry(detId);
1115  GlobalPoint point = cellGeometry->getPosition(); // face of the cell
1116 
1117  if (fracDepth<0) fracDepth=0;
1118  else if (fracDepth>1) fracDepth=1;
1119 
1120  if (fracDepth>0.0) {
1121  CaloCellGeometry::CornersVec cv = cellGeometry->getCorners();
1122  GlobalPoint backPoint = GlobalPoint( 0.25*( cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x() ),
1123  0.25*( cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y() ),
1124  0.25*( cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z() ) );
1125  point += fracDepth * (backPoint-point);
1126  }
1127 
1128  return point;
1129 }
1130 
1132  const CaloCellGeometry* cellGeometry = theGeometry->getGeometry(detId);
1133  GlobalPoint point = cellGeometry->getPosition(); // face of the cell
1134 
1135  if (fracDepth<0) fracDepth=0;
1136  else if (fracDepth>1) fracDepth=1;
1137 
1138  if (fracDepth>0.0) {
1139  CaloCellGeometry::CornersVec cv = cellGeometry->getCorners();
1140  GlobalPoint backPoint = GlobalPoint( 0.25*( cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x() ),
1141  0.25*( cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y() ),
1142  0.25*( cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z() ) );
1143  point += fracDepth * (backPoint-point);
1144  }
1145 
1146  return point;
1147 }
1148 
1149 
1150 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(std::vector<std::pair<DetId,double> >& metaContains,
1151  float fracDepth, double hadE) {
1152 
1153  // this is based on available RecHits, can lead to different actual depths if
1154  // hits in multi-depth towers are not all there
1155  if (hadE<=0) return GlobalPoint(0,0,0);
1156 
1157  double hadX = 0.0;
1158  double hadY = 0.0;
1159  double hadZ = 0.0;
1160 
1161  int nConst = 0;
1162 
1163  std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
1164  for (; mc_it!=metaContains.end(); ++mc_it) {
1165  if (mc_it->first.det() != DetId::Hcal) continue;
1166  // do not use HO for deirection calculations for now
1167  if (HcalDetId(mc_it->first).subdet() == HcalOuter) continue;
1168  ++nConst;
1169 
1170  GlobalPoint p = hadSegmentShwrPos(mc_it->first, fracDepth);
1171 
1172  // longitudinal segmentation: do not weight by energy,
1173  // get the geometrical position
1174  hadX += p.x();
1175  hadY += p.y();
1176  hadZ += p.z();
1177  }
1178 
1179  return GlobalPoint(hadX/nConst, hadY/nConst, hadZ/nConst);
1180 }
1181 
1182 
1184 
1185  // set depth using geometry of cells that are associated with the
1186  // tower (regardless if they have non-zero energies)
1187 
1188 // if (hadE <= 0) return GlobalPoint(0, 0, 0);
1189 
1190  if (fracDepth < 0) fracDepth = 0;
1191  else if (fracDepth > 1) fracDepth = 1;
1192 
1193  GlobalPoint point(0,0,0);
1194 
1195  int iEta = towerId.ieta();
1196  int iPhi = towerId.iphi();
1197 
1198  HcalDetId frontCellId, backCellId;
1199 
1200  if (towerId.ietaAbs() <= 14) {
1201  // barrel, one depth only
1202  frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
1203  backCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
1204  }
1205  else if (towerId.ietaAbs() == 15) {
1206  // barrel, two depths
1207  frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
1208  backCellId = HcalDetId(HcalBarrel, iEta, iPhi, 2);
1209  }
1210  else if (towerId.ietaAbs() == 16) {
1211  // barrel and endcap: two depths HB, one depth HE
1212  frontCellId = HcalDetId(HcalBarrel, iEta, iPhi, 1);
1213  backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 3); // this cell is in endcap!
1214  }
1215  else if (towerId.ietaAbs() == 17) {
1216  // endcap, one depth only
1217  frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
1218  backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
1219  }
1220  else if (towerId.ietaAbs() >= 18 && towerId.ietaAbs() <= 26) {
1221  // endcap: two depths
1222  frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
1223  backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 2);
1224  }
1225  else if (towerId.ietaAbs() <= 29) {
1226  // endcap: three depths
1227  frontCellId = HcalDetId(HcalEndcap, iEta, iPhi, 1);
1228  // there is no iEta=29 for depth 3
1229  if (iEta == 29) iEta = 28;
1230  if (iEta == -29) iEta = -28;
1231  backCellId = HcalDetId(HcalEndcap, iEta, iPhi, 3);
1232  }
1233  else if (towerId.ietaAbs() >= 30) {
1234  // forward, take the goemetry for long fibers
1235  frontCellId = HcalDetId(HcalForward, iEta, iPhi, 1);
1236  backCellId = HcalDetId(HcalForward, iEta, iPhi, 1);
1237  }
1238  else {
1239  // should not get here
1240  return point;
1241  }
1242 
1243  point = hadShwPosFromCells(DetId(frontCellId), DetId(backCellId), fracDepth);
1244 
1245  return point;
1246 }
1247 
1248 GlobalPoint CaloTowersCreationAlgo::hadShwPosFromCells(DetId frontCellId, DetId backCellId, float fracDepth) {
1249 
1250  // uses the "front" and "back" cells
1251  // to determine the axis. point set by the predefined depth.
1252 
1253  const CaloCellGeometry* frontCellGeometry = theGeometry->getGeometry(DetId(frontCellId));
1254  const CaloCellGeometry* backCellGeometry = theGeometry->getGeometry(DetId(backCellId));
1255 
1256  GlobalPoint point = frontCellGeometry->getPosition();
1257 
1258  CaloCellGeometry::CornersVec cv = backCellGeometry->getCorners();
1259 
1260  GlobalPoint backPoint = GlobalPoint(0.25 * (cv[4].x() + cv[5].x() + cv[6].x() + cv[7].x()),
1261  0.25 * (cv[4].y() + cv[5].y() + cv[6].y() + cv[7].y()),
1262  0.25 * (cv[4].z() + cv[5].z() + cv[6].z() + cv[7].z()));
1263 
1264  point += fracDepth * (backPoint - point);
1265 
1266  return point;
1267 }
1268 
1269 
1270 GlobalPoint CaloTowersCreationAlgo::emShwrPos(std::vector<std::pair<DetId,double> >& metaContains,
1271  float fracDepth, double emE) {
1272 
1273  if (emE<=0) return GlobalPoint(0,0,0);
1274 
1275  double emX = 0.0;
1276  double emY = 0.0;
1277  double emZ = 0.0;
1278 
1279  double eSum = 0;
1280 
1281  std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
1282  for (; mc_it!=metaContains.end(); ++mc_it) {
1283  if (mc_it->first.det() != DetId::Ecal) continue;
1284  GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
1285  double e = mc_it->second;
1286 
1287  if (e>0) {
1288  emX += p.x() * e;
1289  emY += p.y() * e;
1290  emZ += p.z() * e;
1291  eSum += e;
1292  }
1293 
1294  }
1295 
1296  return GlobalPoint(emX/eSum, emY/eSum, emZ/eSum);
1297 }
1298 
1299 
1300 GlobalPoint CaloTowersCreationAlgo::emShwrLogWeightPos(std::vector<std::pair<DetId,double> >& metaContains,
1301  float fracDepth, double emE) {
1302 
1303  double emX = 0.0;
1304  double emY = 0.0;
1305  double emZ = 0.0;
1306 
1307  double weight = 0;
1308  double sumWeights = 0;
1309  double sumEmE = 0; // add crystals with E/E_EM > 1.5%
1310  double crystalThresh = 0.015 * emE;
1311 
1312  std::vector<std::pair<DetId,double> >::iterator mc_it = metaContains.begin();
1313  for (; mc_it!=metaContains.end(); ++mc_it) {
1314  if (mc_it->second < 0) continue;
1315  if (mc_it->first.det() == DetId::Ecal && mc_it->second > crystalThresh) sumEmE += mc_it->second;
1316  }
1317 
1318  for (mc_it = metaContains.begin(); mc_it!=metaContains.end(); ++mc_it) {
1319 
1320  if (mc_it->first.det() != DetId::Ecal || mc_it->second < crystalThresh) continue;
1321 
1322  GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
1323 
1324  weight = 4.2 + log(mc_it->second/sumEmE);
1325  sumWeights += weight;
1326 
1327  emX += p.x() * weight;
1328  emY += p.y() * weight;
1329  emZ += p.z() * weight;
1330  }
1331 
1332  return GlobalPoint(emX/sumWeights, emY/sumWeights, emZ/sumWeights);
1333 }
1334 
1335 
1336 
1337 
1338 
1340 
1341  const float timeUnit = 0.01; // discretization (ns)
1342 
1343  if (time> 300.0) return 30000;
1344  if (time< -300.0) return -30000;
1345 
1346  return int(time/timeUnit + 0.5);
1347 
1348 }
1349 
1350 
1351 
1352 //========================================================
1353 //
1354 // Bad/anomolous cell handling
1355 
1356 
1357 
1358 
1360 
1361  // This method fills the map of number of dead channels for the calotower,
1362  // The key of the map is CaloTowerDetId.
1363  // By definition these channels are not going to be in the RecHit collections.
1364 
1365  std::vector<DetId> allChanInStatusCont = theHcalChStatus->getAllChannels();
1366 
1367  for (std::vector<DetId>::iterator it = allChanInStatusCont.begin(); it!=allChanInStatusCont.end(); ++it) {
1368 
1369  const uint32_t dbStatusFlag = theHcalChStatus->getValues(*it)->getValue();
1370 
1371  if (theHcalSevLvlComputer->dropChannel(dbStatusFlag)) {
1372 
1374 
1375  hcalDropChMap[twrId] +=1;
1376 
1377  // special case for tower 29: if HCAL hit is in depth 3 add to twr 29 as well
1378  if (HcalDetId(*it).subdet()==HcalEndcap &&
1379  HcalDetId(*it).depth()==3 &&
1380  HcalDetId(*it).ietaAbs()==28) {
1381 
1382  CaloTowerDetId twrId29 = CaloTowerDetId(twrId.ieta()+twrId.zside(), twrId.iphi());
1383  hcalDropChMap[twrId29] +=1;
1384  }
1385 
1386  }
1387 
1388  }
1389 
1390  return;
1391 }
1392 
1393 
1394 
1396 
1398 
1399  const DetId id = hit->detid();
1400 
1401  const uint32_t recHitFlag = hit->flags();
1402  const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1403 
1404  int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1405  bool isRecovered = theHcalSevLvlComputer->recoveredRecHit(id, recHitFlag);
1406 
1407 
1408  // For use with hits rejected in the default reconstruction
1409  if (useRejectedHitsOnly) {
1410 
1411  if (!isRecovered) {
1412  if (severityLevel <= int(theHcalAcceptSeverityLevel) ||
1414  // this hit was either already accepted or is worse than
1415  }
1416  else {
1417 
1419  // skip recovered hits either because they were already used or because there was an explicit instruction
1421  }
1422  else if (useRejectedRecoveredHcalHits) {
1424  }
1425 
1426  } // recovered channels
1427 
1428  // clasify channels as problematic: no good hits are supposed to be present in the
1429  // extra rechit collections
1431 
1432  } // treatment of rejected hits
1433 
1434 
1435 
1436 
1437  // this is for the regular reconstruction sequence
1438 
1439  if (severityLevel == 0) return CaloTowersCreationAlgo::GoodChan;
1440 
1441  if (isRecovered) {
1442  return (theRecoveredHcalHitsAreUsed) ?
1444  }
1445  else {
1446  if (severityLevel > int(theHcalAcceptSeverityLevel)) {
1448  }
1449  else {
1451  }
1452  }
1453 
1454 }
1455 
1456 
1457 
1459 
1460  const DetId id = hit->detid();
1461 
1462  // uint16_t dbStatus = theEcalChStatus->find(id)->getStatusCode();
1463  // uint32_t rhFlags = hit->flags();
1464  // int severityLevel = theEcalSevLvlAlgo->severityLevel(rhFlags, dbStatus);
1465  // The methods above will become private and cannot be usef for flagging ecal spikes.
1466  // Use the recommended interface - we leave the parameters for spilke removal to be specified by ECAL.
1467 
1468 
1469  int severityLevel = 999;
1470 
1471  if (id.subdetId() == EcalBarrel) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEbHandle, *theEcalChStatus);
1472  else if (id.subdetId() == EcalEndcap) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEeHandle, *theEcalChStatus);
1473 
1474  // there should be no other ECAL types used in this reconstruction
1475 
1476  // The definition of ECAL severity levels uses categories that
1477  // are similar to the defined for CaloTower. (However, the categorization
1478  // for CaloTowers depends on the specified maximum acceptabel severity and therefore cannnot
1479  // be exact correspondence between the two. ECAL has additional categories describing modes of failure.)
1480  // This approach is different from the initial idea and from
1481  // the implementation for HCAL. Still make the logic similar to HCAL so that one has the ability to
1482  // exclude problematic channels as defined by ECAL.
1483  // For definitions of ECAL severity levels see RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h
1484 
1485 
1486  bool isRecovered = (severityLevel == EcalSeverityLevelAlgo::kRecovered);
1487 
1488 
1489  // for use with rejected by the regular reconstruction hits
1490  if (useRejectedHitsOnly) {
1491 
1492  if (!isRecovered) {
1493  if (severityLevel <= int(theEcalAcceptSeverityLevel) ||
1495  // this hit was either accepted already, or counted as bad under the conditions
1496  }
1497  else {
1498 
1500  // skip recovered hits either because they were already used or because there was an explicit instruction
1502  }
1503  else if (useRejectedRecoveredEcalHits) {
1505  }
1506 
1507  } // recovered channels
1508 
1509  // clasify channels as problematic
1511 
1512  } // treatment of rejected hits
1513 
1514 
1515 
1516 
1517 
1518  // for normal reconstruction
1520 
1521  if (isRecovered) {
1522  return (theRecoveredEcalHitsAreUsed) ?
1524  }
1525  else {
1526  if (severityLevel > int(theEcalAcceptSeverityLevel)) {
1528  }
1529  else {
1531  }
1532  }
1533 
1534 
1535 }
1536 
std::vector< double > theHBGrid
const EcalChannelStatus * theEcalChStatus
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:43
int i
Definition: DBlmapReader.cc:9
const HcalChannelQuality * theHcalChStatus
size_t constituentsSize() const
Definition: CaloTower.h:74
unsigned int theEcalAcceptSeverityLevelForRejectedHit
bool contains(EventRange const &lh, EventID const &rh)
Definition: EventRange.cc:38
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:32
int ietaAbs() const
get the absolute value of the tower ieta
const DetId & detid() const
Definition: CaloRecHit.h:21
DetId constituent(size_t i) const
Definition: CaloTower.h:75
void setCaloTowerStatus(unsigned int numBadHcalChan, unsigned int numBadEcalChan, unsigned int numRecHcalChan, unsigned int numRecEcalChan, unsigned int numProbHcalChan, unsigned int numProbEcalChan)
Definition: CaloTower.cc:208
std::vector< double > theHESGrid
std::vector< double > theHEDGrid
GlobalPoint emCrystalShwrPos(DetId detId, float fracDepth)
float ecalTime() const
Definition: CaloTower.h:150
math::PtEtaPhiMLorentzVector PolarLorentzVector
Lorentz vector.
Definition: LeafCandidate.h:27
Geom::Phi< T > phi() const
Definition: PV3DBase.h:63
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
std::vector< T >::const_iterator const_iterator
T y() const
Definition: PV3DBase.h:57
void push_back(T const &t)
void finish(CaloTowerCollection &destCollection)
GlobalPoint emShwrLogWeightPos(std::vector< std::pair< DetId, double > > &metaContains, float fracDepth, double totEmE)
edm::Handle< EcalRecHitCollection > theEeHandle
std::vector< double > theHOWeights
MetaTower & find(const CaloTowerDetId &id)
looks for a given tower in the internal cache. If it can&#39;t find it, it makes it.
std::vector< double > theEEGrid
void rescale(const CaloTower *ct)
float time() const
Definition: CaloRecHit.h:20
CaloTower convert(const CaloTowerDetId &id, const MetaTower &mt)
unsigned int ecalChanStatusForCaloTower(const CaloRecHit *hit)
T eta() const
void setHottestCellE(double e)
Definition: CaloTower.h:66
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
std::vector< DetId > constituentsOf(const CaloTowerDetId &id) const
Get the constituent detids for this tower id ( not yet implemented )
std::vector< double > theEEWeights
GlobalPoint hadShwrPos(std::vector< std::pair< DetId, double > > &metaContains, float fracDepth, double hadE)
std::vector< double > theHESWeights
U second(std::pair< T, U > const &p)
const CaloSubdetectorGeometry * theTowerGeometry
virtual double energy() const
energy
int depth() const
get the tower depth
Definition: HcalDetId.h:42
void addConstituents(const std::vector< DetId > &ids)
Definition: CaloTower.cc:186
unsigned int hcalChanStatusForCaloTower(const CaloRecHit *hit)
CaloTowerDetId towerOf(const DetId &id) const
Get the tower id for this det id (or null if not known)
Definition: DDAxes.h:10
double emEnergy() const
Definition: CaloTower.h:79
std::vector< double > theHF2Grid
std::vector< double > theHEDWeights
void rescaleTowers(const CaloTowerCollection &ctInput, CaloTowerCollection &ctResult)
static const int SubdetId
float threshold
Definition: crabWrap.py:319
float energy() const
Definition: CaloRecHit.h:19
std::vector< DetId > getAllChannels() const
T z() const
Definition: PV3DBase.h:58
bool recoveredRecHit(const DetId &myid, const uint32_t &myflag) const
tuple result
Definition: query.py:137
const CaloGeometry * theGeometry
uint32_t flags() const
Definition: CaloRecHit.h:22
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
int iphi() const
get the tower iphi
HcalSubdetector
Definition: HcalAssistant.h:32
std::vector< double > theHOGrid
bool dropChannel(const uint32_t &mystatus) const
const CaloTowerConstituentsMap * theTowerConstituentsMap
std::vector< double > theHF1Weights
int firstHEDoublePhiRing() const
Definition: HcalTopology.h:70
GlobalPoint hadSegmentShwrPos(DetId detId, float fracDepth)
std::vector< std::pair< DetId, double > > metaConstituents
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
bool theHOIsUsed
only affects energy and ET calculation. HO is still recorded in the tower
void process(const HBHERecHitCollection &hbhe)
int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:36
double hadEnergy() const
Definition: CaloTower.h:80
const_iterator end() const
GlobalPoint emShwrPos(std::vector< std::pair< DetId, double > > &metaContains, float fracDepth, double totEmE)
std::vector< double > theHBWeights
const HcalSeverityLevelComputer * theHcalSevLvlComputer
Definition: DetId.h:20
CaloTowerDetId id() const
Definition: CaloTower.h:72
float hcalTime() const
Definition: CaloTower.h:151
Log< T >::type log(const T &t)
Definition: Log.h:22
Detector
Definition: DetId.h:26
GlobalPoint hadShwPosFromCells(DetId frontCell, DetId backCell, float fracDepth)
bool null() const
is this a null id ?
Definition: DetId.h:47
int zside() const
get the z-side of the tower (1/-1)
void getThresholdAndWeight(const DetId &detId, double &threshold, double &weight) const
helper method to look up the appropriate threshold &amp; weight
int getSeverityLevel(const DetId &myid, const uint32_t &myflag, const uint32_t &mystatus) const
const HcalTopology * theHcalTopology
T eta() const
Definition: PV3DBase.h:70
void setEcalTime(int t)
Definition: CaloTower.h:53
std::map< CaloTowerDetId, int > hcalDropChMap
const EcalSeverityLevelAlgo * theEcalSevLvlAlgo
const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id.
Definition: CaloGeometry.cc:76
std::vector< double > theHF1Grid
perl if(1 lt scalar(@::datatypes))
Definition: edlooper.cc:31
void assignHit(const CaloRecHit *recHit)
adds a single hit to the tower
std::vector< double > theEBGrid
tuple cout
Definition: gather_cfg.py:41
static int severityLevel(const DetId, const EcalRecHitCollection &, const EcalChannelStatus &, float recHitEtThreshold=5., SpikeId spId=kSwissCross, float spIdThreshold=0.95, float recHitEnergyThresholdForTiming=2., float recHitEnergyThresholdForEE=15, float spIdThresholdIEta85=0.999)
std::vector< double > theHF2Weights
int ieta() const
get the tower ieta
uint32_t getValue() const
unsigned int theHcalAcceptSeverityLevelForRejectedHit
const Item * getValues(DetId fId) const
Detector det() const
get the detector field from this detid
Definition: DetId.h:37
void setHcalTime(int t)
Definition: CaloTower.h:54
void setGeometry(const CaloTowerConstituentsMap *cttopo, const HcalTopology *htopo, const CaloGeometry *geo)
EcalSubdetector
const GlobalPoint & getPosition() const
T x() const
Definition: PV3DBase.h:56
std::vector< double > theEBWeights
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
double outerEnergy() const
Definition: CaloTower.h:81
const_iterator begin() const
virtual const CornersVec & getCorners() const =0
edm::Handle< EcalRecHitCollection > theEbHandle
math::PtEtaPhiMLorentzVector PolarLorentzVector
Lorentz vector.
Definition: Candidate.h:41