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