File indexing completed on 2024-09-13 22:52:41
0001 #include "CaloTowersCreationAlgo.h"
0002
0003 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0004 #include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
0005 #include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
0006 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0007 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "Math/Interpolator.h"
0011
0012 #include <cmath>
0013
0014
0015
0016 CaloTowersCreationAlgo::CaloTowersCreationAlgo()
0017 : theEBthreshold(-1000.),
0018 theEEthreshold(-1000.),
0019
0020 theUseEtEBTresholdFlag(false),
0021 theUseEtEETresholdFlag(false),
0022 theUseSymEBTresholdFlag(false),
0023 theUseSymEETresholdFlag(false),
0024
0025 theHcalThreshold(-1000.),
0026 theHBthreshold(-1000.),
0027 theHBthreshold1(-1000.),
0028 theHBthreshold2(-1000.),
0029 theHESthreshold(-1000.),
0030 theHESthreshold1(-1000.),
0031 theHEDthreshold(-1000.),
0032 theHEDthreshold1(-1000.),
0033 theHOthreshold0(-1000.),
0034 theHOthresholdPlus1(-1000.),
0035 theHOthresholdMinus1(-1000.),
0036 theHOthresholdPlus2(-1000.),
0037 theHOthresholdMinus2(-1000.),
0038 theHF1threshold(-1000.),
0039 theHF2threshold(-1000.),
0040 theEBGrid(std::vector<double>(5, 10.)),
0041 theEBWeights(std::vector<double>(5, 1.)),
0042 theEEGrid(std::vector<double>(5, 10.)),
0043 theEEWeights(std::vector<double>(5, 1.)),
0044 theHBGrid(std::vector<double>(5, 10.)),
0045 theHBWeights(std::vector<double>(5, 1.)),
0046 theHESGrid(std::vector<double>(5, 10.)),
0047 theHESWeights(std::vector<double>(5, 1.)),
0048 theHEDGrid(std::vector<double>(5, 10.)),
0049 theHEDWeights(std::vector<double>(5, 1.)),
0050 theHOGrid(std::vector<double>(5, 10.)),
0051 theHOWeights(std::vector<double>(5, 1.)),
0052 theHF1Grid(std::vector<double>(5, 10.)),
0053 theHF1Weights(std::vector<double>(5, 1.)),
0054 theHF2Grid(std::vector<double>(5, 10.)),
0055 theHF2Weights(std::vector<double>(5, 1.)),
0056 theEBweight(1.),
0057 theEEweight(1.),
0058 theHBweight(1.),
0059 theHESweight(1.),
0060 theHEDweight(1.),
0061 theHOweight(1.),
0062 theHF1weight(1.),
0063 theHF2weight(1.),
0064 theEcutTower(-1000.),
0065 theEBSumThreshold(-1000.),
0066 theEESumThreshold(-1000.),
0067 theEBEScale(50.),
0068 theEEEScale(50.),
0069 theHBEScale(50.),
0070 theHESEScale(50.),
0071 theHEDEScale(50.),
0072 theHOEScale(50.),
0073 theHF1EScale(50.),
0074 theHF2EScale(50.),
0075 theHcalTopology(nullptr),
0076 theGeometry(nullptr),
0077 theTowerConstituentsMap(nullptr),
0078 theHcalAcceptSeverityLevel(0),
0079 theRecoveredHcalHitsAreUsed(false),
0080 theRecoveredEcalHitsAreUsed(false),
0081 useRejectedHitsOnly(false),
0082 theHcalAcceptSeverityLevelForRejectedHit(0),
0083 useRejectedRecoveredHcalHits(0),
0084 useRejectedRecoveredEcalHits(0),
0085 missingHcalRescaleFactorForEcal(0.),
0086 theHOIsUsed(true),
0087
0088 theMomConstrMethod(0),
0089 theMomHBDepth(0.),
0090 theMomHEDepth(0.),
0091 theMomEBDepth(0.),
0092 theMomEEDepth(0.),
0093 theHcalPhase(0) {}
0094
0095 CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold,
0096 double EEthreshold,
0097
0098 bool useEtEBTreshold,
0099 bool useEtEETreshold,
0100 bool useSymEBTreshold,
0101 bool useSymEETreshold,
0102
0103 double HcalThreshold,
0104 double HBthreshold,
0105 double HBthreshold1,
0106 double HBthreshold2,
0107 double HESthreshold,
0108 double HESthreshold1,
0109 double HEDthreshold,
0110 double HEDthreshold1,
0111 double HOthreshold0,
0112 double HOthresholdPlus1,
0113 double HOthresholdMinus1,
0114 double HOthresholdPlus2,
0115 double HOthresholdMinus2,
0116 double HF1threshold,
0117 double HF2threshold,
0118 double EBweight,
0119 double EEweight,
0120 double HBweight,
0121 double HESweight,
0122 double HEDweight,
0123 double HOweight,
0124 double HF1weight,
0125 double HF2weight,
0126 double EcutTower,
0127 double EBSumThreshold,
0128 double EESumThreshold,
0129 bool useHO,
0130
0131 int momConstrMethod,
0132 double momHBDepth,
0133 double momHEDepth,
0134 double momEBDepth,
0135 double momEEDepth,
0136 int hcalPhase)
0137
0138 : theEBthreshold(EBthreshold),
0139 theEEthreshold(EEthreshold),
0140
0141 theUseEtEBTresholdFlag(useEtEBTreshold),
0142 theUseEtEETresholdFlag(useEtEETreshold),
0143 theUseSymEBTresholdFlag(useSymEBTreshold),
0144 theUseSymEETresholdFlag(useSymEETreshold),
0145
0146 theHcalThreshold(HcalThreshold),
0147 theHBthreshold(HBthreshold),
0148 theHBthreshold1(HBthreshold1),
0149 theHBthreshold2(HBthreshold2),
0150 theHESthreshold(HESthreshold),
0151 theHESthreshold1(HESthreshold1),
0152 theHEDthreshold(HEDthreshold),
0153 theHEDthreshold1(HEDthreshold1),
0154 theHOthreshold0(HOthreshold0),
0155 theHOthresholdPlus1(HOthresholdPlus1),
0156 theHOthresholdMinus1(HOthresholdMinus1),
0157 theHOthresholdPlus2(HOthresholdPlus2),
0158 theHOthresholdMinus2(HOthresholdMinus2),
0159 theHF1threshold(HF1threshold),
0160 theHF2threshold(HF2threshold),
0161 theEBGrid(std::vector<double>(5, 10.)),
0162 theEBWeights(std::vector<double>(5, 1.)),
0163 theEEGrid(std::vector<double>(5, 10.)),
0164 theEEWeights(std::vector<double>(5, 1.)),
0165 theHBGrid(std::vector<double>(5, 10.)),
0166 theHBWeights(std::vector<double>(5, 1.)),
0167 theHESGrid(std::vector<double>(5, 10.)),
0168 theHESWeights(std::vector<double>(5, 1.)),
0169 theHEDGrid(std::vector<double>(5, 10.)),
0170 theHEDWeights(std::vector<double>(5, 1.)),
0171 theHOGrid(std::vector<double>(5, 10.)),
0172 theHOWeights(std::vector<double>(5, 1.)),
0173 theHF1Grid(std::vector<double>(5, 10.)),
0174 theHF1Weights(std::vector<double>(5, 1.)),
0175 theHF2Grid(std::vector<double>(5, 10.)),
0176 theHF2Weights(std::vector<double>(5, 1.)),
0177 theEBweight(EBweight),
0178 theEEweight(EEweight),
0179 theHBweight(HBweight),
0180 theHESweight(HESweight),
0181 theHEDweight(HEDweight),
0182 theHOweight(HOweight),
0183 theHF1weight(HF1weight),
0184 theHF2weight(HF2weight),
0185 theEcutTower(EcutTower),
0186 theEBSumThreshold(EBSumThreshold),
0187 theEESumThreshold(EESumThreshold),
0188 theEBEScale(50.),
0189 theEEEScale(50.),
0190 theHBEScale(50.),
0191 theHESEScale(50.),
0192 theHEDEScale(50.),
0193 theHOEScale(50.),
0194 theHF1EScale(50.),
0195 theHF2EScale(50.),
0196 theHcalTopology(nullptr),
0197 theGeometry(nullptr),
0198 theTowerConstituentsMap(nullptr),
0199 theHcalAcceptSeverityLevel(0),
0200 theRecoveredHcalHitsAreUsed(false),
0201 theRecoveredEcalHitsAreUsed(false),
0202 useRejectedHitsOnly(false),
0203 theHcalAcceptSeverityLevelForRejectedHit(0),
0204 useRejectedRecoveredHcalHits(0),
0205 useRejectedRecoveredEcalHits(0),
0206 missingHcalRescaleFactorForEcal(0.),
0207 theHOIsUsed(useHO),
0208
0209 theMomConstrMethod(momConstrMethod),
0210 theMomHBDepth(momHBDepth),
0211 theMomHEDepth(momHEDepth),
0212 theMomEBDepth(momEBDepth),
0213 theMomEEDepth(momEEDepth),
0214 theHcalPhase(hcalPhase) {}
0215
0216 CaloTowersCreationAlgo::CaloTowersCreationAlgo(double EBthreshold,
0217 double EEthreshold,
0218 bool useEtEBTreshold,
0219 bool useEtEETreshold,
0220 bool useSymEBTreshold,
0221 bool useSymEETreshold,
0222
0223 double HcalThreshold,
0224 double HBthreshold,
0225 double HBthreshold1,
0226 double HBthreshold2,
0227 double HESthreshold,
0228 double HESthreshold1,
0229 double HEDthreshold,
0230 double HEDthreshold1,
0231 double HOthreshold0,
0232 double HOthresholdPlus1,
0233 double HOthresholdMinus1,
0234 double HOthresholdPlus2,
0235 double HOthresholdMinus2,
0236 double HF1threshold,
0237 double HF2threshold,
0238 const std::vector<double>& EBGrid,
0239 const std::vector<double>& EBWeights,
0240 const std::vector<double>& EEGrid,
0241 const std::vector<double>& EEWeights,
0242 const std::vector<double>& HBGrid,
0243 const std::vector<double>& HBWeights,
0244 const std::vector<double>& HESGrid,
0245 const std::vector<double>& HESWeights,
0246 const std::vector<double>& HEDGrid,
0247 const std::vector<double>& HEDWeights,
0248 const std::vector<double>& HOGrid,
0249 const std::vector<double>& HOWeights,
0250 const std::vector<double>& HF1Grid,
0251 const std::vector<double>& HF1Weights,
0252 const std::vector<double>& HF2Grid,
0253 const std::vector<double>& HF2Weights,
0254 double EBweight,
0255 double EEweight,
0256 double HBweight,
0257 double HESweight,
0258 double HEDweight,
0259 double HOweight,
0260 double HF1weight,
0261 double HF2weight,
0262 double EcutTower,
0263 double EBSumThreshold,
0264 double EESumThreshold,
0265 bool useHO,
0266
0267 int momConstrMethod,
0268 double momHBDepth,
0269 double momHEDepth,
0270 double momEBDepth,
0271 double momEEDepth,
0272 int hcalPhase)
0273
0274 : theEBthreshold(EBthreshold),
0275 theEEthreshold(EEthreshold),
0276
0277 theUseEtEBTresholdFlag(useEtEBTreshold),
0278 theUseEtEETresholdFlag(useEtEETreshold),
0279 theUseSymEBTresholdFlag(useSymEBTreshold),
0280 theUseSymEETresholdFlag(useSymEETreshold),
0281
0282 theHcalThreshold(HcalThreshold),
0283 theHBthreshold(HBthreshold),
0284 theHBthreshold1(HBthreshold1),
0285 theHBthreshold2(HBthreshold2),
0286 theHESthreshold(HESthreshold),
0287 theHESthreshold1(HESthreshold1),
0288 theHEDthreshold(HEDthreshold),
0289 theHEDthreshold1(HEDthreshold1),
0290 theHOthreshold0(HOthreshold0),
0291 theHOthresholdPlus1(HOthresholdPlus1),
0292 theHOthresholdMinus1(HOthresholdMinus1),
0293 theHOthresholdPlus2(HOthresholdPlus2),
0294 theHOthresholdMinus2(HOthresholdMinus2),
0295 theHF1threshold(HF1threshold),
0296 theHF2threshold(HF2threshold),
0297 theEBGrid(EBGrid),
0298 theEBWeights(EBWeights),
0299 theEEGrid(EEGrid),
0300 theEEWeights(EEWeights),
0301 theHBGrid(HBGrid),
0302 theHBWeights(HBWeights),
0303 theHESGrid(HESGrid),
0304 theHESWeights(HESWeights),
0305 theHEDGrid(HEDGrid),
0306 theHEDWeights(HEDWeights),
0307 theHOGrid(HOGrid),
0308 theHOWeights(HOWeights),
0309 theHF1Grid(HF1Grid),
0310 theHF1Weights(HF1Weights),
0311 theHF2Grid(HF2Grid),
0312 theHF2Weights(HF2Weights),
0313 theEBweight(EBweight),
0314 theEEweight(EEweight),
0315 theHBweight(HBweight),
0316 theHESweight(HESweight),
0317 theHEDweight(HEDweight),
0318 theHOweight(HOweight),
0319 theHF1weight(HF1weight),
0320 theHF2weight(HF2weight),
0321 theEcutTower(EcutTower),
0322 theEBSumThreshold(EBSumThreshold),
0323 theEESumThreshold(EESumThreshold),
0324 theEBEScale(50.),
0325 theEEEScale(50.),
0326 theHBEScale(50.),
0327 theHESEScale(50.),
0328 theHEDEScale(50.),
0329 theHOEScale(50.),
0330 theHF1EScale(50.),
0331 theHF2EScale(50.),
0332 theHcalTopology(nullptr),
0333 theGeometry(nullptr),
0334 theTowerConstituentsMap(nullptr),
0335 theHcalAcceptSeverityLevel(0),
0336 theRecoveredHcalHitsAreUsed(false),
0337 theRecoveredEcalHitsAreUsed(false),
0338 useRejectedHitsOnly(false),
0339 theHcalAcceptSeverityLevelForRejectedHit(0),
0340 useRejectedRecoveredHcalHits(0),
0341 useRejectedRecoveredEcalHits(0),
0342 missingHcalRescaleFactorForEcal(0.),
0343 theHOIsUsed(useHO),
0344
0345 theMomConstrMethod(momConstrMethod),
0346 theMomHBDepth(momHBDepth),
0347 theMomHEDepth(momHEDepth),
0348 theMomEBDepth(momEBDepth),
0349 theMomEEDepth(momEEDepth),
0350 theHcalPhase(hcalPhase) {
0351
0352
0353
0354 }
0355
0356 void CaloTowersCreationAlgo::setThresFromDB(const EcalPFRecHitThresholds* EcalCuts, const HcalPFCuts* HcalCuts) {
0357 ecalCuts = EcalCuts;
0358 hcalCuts = HcalCuts;
0359 }
0360
0361 void CaloTowersCreationAlgo::setGeometry(const CaloTowerTopology* cttopo,
0362 const CaloTowerConstituentsMap* ctmap,
0363 const HcalTopology* htopo,
0364 const CaloGeometry* geo) {
0365 theTowerTopology = cttopo;
0366 theTowerConstituentsMap = ctmap;
0367 theHcalTopology = htopo;
0368 theGeometry = geo;
0369 theTowerGeometry = geo->getSubdetectorGeometry(DetId::Calo, CaloTowerDetId::SubdetId);
0370
0371
0372 ecalBadChs.resize(theTowerTopology->sizeForDenseIndexing(), 0);
0373 }
0374
0375 void CaloTowersCreationAlgo::begin() {
0376 theTowerMap.clear();
0377 theTowerMapSize = 0;
0378
0379 }
0380
0381 void CaloTowersCreationAlgo::process(const HBHERecHitCollection& hbhe) {
0382 for (HBHERecHitCollection::const_iterator hbheItr = hbhe.begin(); hbheItr != hbhe.end(); ++hbheItr)
0383 assignHitHcal(&(*hbheItr));
0384 }
0385
0386 void CaloTowersCreationAlgo::process(const HORecHitCollection& ho) {
0387 for (HORecHitCollection::const_iterator hoItr = ho.begin(); hoItr != ho.end(); ++hoItr)
0388 assignHitHcal(&(*hoItr));
0389 }
0390
0391 void CaloTowersCreationAlgo::process(const HFRecHitCollection& hf) {
0392 for (HFRecHitCollection::const_iterator hfItr = hf.begin(); hfItr != hf.end(); ++hfItr)
0393 assignHitHcal(&(*hfItr));
0394 }
0395
0396 void CaloTowersCreationAlgo::process(const EcalRecHitCollection& ec) {
0397 for (EcalRecHitCollection::const_iterator ecItr = ec.begin(); ecItr != ec.end(); ++ecItr)
0398 assignHitEcal(&(*ecItr));
0399 }
0400
0401
0402
0403
0404
0405 void CaloTowersCreationAlgo::process(const CaloTowerCollection& ctc) {
0406 for (CaloTowerCollection::const_iterator ctcItr = ctc.begin(); ctcItr != ctc.end(); ++ctcItr) {
0407 rescale(&(*ctcItr));
0408 }
0409 }
0410
0411 void CaloTowersCreationAlgo::finish(CaloTowerCollection& result) {
0412
0413 result.reserve(theTowerMapSize);
0414
0415
0416
0417
0418 for (auto const& mt : theTowerMap) {
0419
0420
0421 if (!mt.empty()) {
0422 convert(mt.id, mt, result);
0423 }
0424 }
0425
0426
0427
0428
0429 theTowerMap.clear();
0430 theTowerMapSize = 0;
0431 }
0432
0433 void CaloTowersCreationAlgo::rescaleTowers(const CaloTowerCollection& ctc, CaloTowerCollection& ctcResult) {
0434 for (CaloTowerCollection::const_iterator ctcItr = ctc.begin(); ctcItr != ctc.end(); ++ctcItr) {
0435 CaloTowerDetId twrId = ctcItr->id();
0436 double newE_em = ctcItr->emEnergy();
0437 double newE_had = ctcItr->hadEnergy();
0438 double newE_outer = ctcItr->outerEnergy();
0439
0440 double threshold = 0.0;
0441 double weight = 1.0;
0442
0443
0444 if (ctcItr->ietaAbs() >= theTowerTopology->firstHFRing()) {
0445 double E_short = 0.5 * newE_had;
0446 double E_long = newE_em + 0.5 * newE_had;
0447
0448 E_long *= theHF1weight;
0449 E_short *= theHF2weight;
0450
0451 newE_em = E_long - E_short;
0452 newE_had = 2.0 * E_short;
0453 }
0454
0455 else {
0456
0457
0458 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0459 DetId constId = ctcItr->constituent(iConst);
0460 if (constId.det() != DetId::Ecal)
0461 continue;
0462 getThresholdAndWeight(constId, threshold, weight);
0463 newE_em *= weight;
0464 break;
0465 }
0466
0467 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0468 DetId constId = ctcItr->constituent(iConst);
0469 if (constId.det() != DetId::Hcal)
0470 continue;
0471 if (HcalDetId(constId).subdet() != HcalOuter)
0472 continue;
0473 getThresholdAndWeight(constId, threshold, weight);
0474 newE_outer *= weight;
0475 break;
0476 }
0477
0478 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0479 DetId constId = ctcItr->constituent(iConst);
0480 if (constId.det() != DetId::Hcal)
0481 continue;
0482 if (HcalDetId(constId).subdet() == HcalOuter)
0483 continue;
0484 getThresholdAndWeight(constId, threshold, weight);
0485 newE_had *= weight;
0486 if (ctcItr->ietaAbs() > theTowerTopology->firstHERing())
0487 newE_outer *= weight;
0488 break;
0489 }
0490
0491 }
0492
0493
0494
0495 double newE_hadTot =
0496 (theHOIsUsed && twrId.ietaAbs() <= theTowerTopology->lastHORing()) ? newE_had + newE_outer : newE_had;
0497
0498 GlobalPoint emPoint = ctcItr->emPosition();
0499 GlobalPoint hadPoint = ctcItr->emPosition();
0500
0501 double f_em = 1.0 / cosh(emPoint.eta());
0502 double f_had = 1.0 / cosh(hadPoint.eta());
0503
0504 CaloTower::PolarLorentzVector towerP4;
0505
0506 if (ctcItr->ietaAbs() < theTowerTopology->firstHFRing()) {
0507 if (newE_em > 0)
0508 towerP4 += CaloTower::PolarLorentzVector(newE_em * f_em, emPoint.eta(), emPoint.phi(), 0);
0509 if (newE_hadTot > 0)
0510 towerP4 += CaloTower::PolarLorentzVector(newE_hadTot * f_had, hadPoint.eta(), hadPoint.phi(), 0);
0511 } else {
0512 double newE_tot = newE_em + newE_had;
0513
0514 if (newE_tot > 0)
0515 towerP4 += CaloTower::PolarLorentzVector(newE_tot * f_had, hadPoint.eta(), hadPoint.phi(), 0);
0516 }
0517
0518 CaloTower rescaledTower(twrId, newE_em, newE_had, newE_outer, -1, -1, towerP4, emPoint, hadPoint);
0519
0520 rescaledTower.setEcalTime(int(ctcItr->ecalTime() * 100.0 + 0.5));
0521 rescaledTower.setHcalTime(int(ctcItr->hcalTime() * 100.0 + 0.5));
0522
0523 rescaledTower.setHcalSubdet(theTowerTopology->lastHBRing(),
0524 theTowerTopology->lastHERing(),
0525 theTowerTopology->lastHFRing(),
0526 theTowerTopology->lastHORing());
0527
0528 std::vector<DetId> contains;
0529 for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0530 contains.push_back(ctcItr->constituent(iConst));
0531 }
0532 rescaledTower.addConstituents(contains);
0533
0534 rescaledTower.setCaloTowerStatus(ctcItr->towerStatusWord());
0535
0536 ctcResult.push_back(rescaledTower);
0537
0538 }
0539 }
0540
0541 void CaloTowersCreationAlgo::assignHitHcal(const CaloRecHit* recHit) {
0542 DetId detId = recHit->detid();
0543 DetId detIdF(detId);
0544 if (detId.det() == DetId::Hcal && theHcalTopology->getMergePositionFlag()) {
0545 detIdF = theHcalTopology->idFront(HcalDetId(detId));
0546 #ifdef EDM_ML_DEBUG
0547 std::cout << "AssignHitHcal DetId " << HcalDetId(detId) << " Front " << HcalDetId(detIdF) << std::endl;
0548 #endif
0549 }
0550
0551 unsigned int chStatusForCT = hcalChanStatusForCaloTower(recHit);
0552
0553
0554
0555 if (chStatusForCT == CaloTowersCreationAlgo::IgnoredChan)
0556 return;
0557
0558 double threshold, weight;
0559 getThresholdAndWeight(detId, threshold, weight);
0560
0561 double energy = recHit->energy();
0562 double e = energy * weight;
0563
0564
0565 bool merge(false);
0566 if (detIdF.det() == DetId::Hcal && HcalDetId(detIdF).subdet() == HcalEndcap &&
0567 (theHcalPhase == 0 || theHcalPhase == 1) &&
0568
0569 HcalDetId(detIdF).ietaAbs() == theHcalTopology->lastHERing() - 1) {
0570 merge = theHcalTopology->mergedDepth29(HcalDetId(detIdF));
0571 #ifdef EDM_ML_DEBUG
0572 std::cout << "Merge " << HcalDetId(detIdF) << ":" << merge << std::endl;
0573 #endif
0574 }
0575 if (merge) {
0576
0577
0578
0579
0580 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0581 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0582 if (towerDetId.null())
0583 return;
0584 MetaTower& tower28 = find(towerDetId);
0585 CaloTowerDetId towerDetId29(towerDetId.ieta() + towerDetId.zside(), towerDetId.iphi());
0586 MetaTower& tower29 = find(towerDetId29);
0587 tower28.numBadHcalCells += 1;
0588 tower29.numBadHcalCells += 1;
0589 }
0590
0591 else if (0.5 * energy >= threshold) {
0592
0593 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0594 if (towerDetId.null())
0595 return;
0596 MetaTower& tower28 = find(towerDetId);
0597 CaloTowerDetId towerDetId29(towerDetId.ieta() + towerDetId.zside(), towerDetId.iphi());
0598 MetaTower& tower29 = find(towerDetId29);
0599
0600 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0601 tower28.numRecHcalCells += 1;
0602 tower29.numRecHcalCells += 1;
0603 } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0604 tower28.numProbHcalCells += 1;
0605 tower29.numProbHcalCells += 1;
0606 }
0607
0608
0609 double e28 = 0.5 * e;
0610 double e29 = 0.5 * e;
0611
0612 tower28.E_had += e28;
0613 tower28.E += e28;
0614 std::pair<DetId, float> mc(detId, e28);
0615 tower28.metaConstituents.push_back(mc);
0616
0617 tower29.E_had += e29;
0618 tower29.E += e29;
0619 tower29.metaConstituents.push_back(mc);
0620
0621
0622
0623
0624 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0625 tower28.hadSumTimeTimesE += (e28 * recHit->time());
0626 tower28.hadSumEForTime += e28;
0627 tower29.hadSumTimeTimesE += (e29 * recHit->time());
0628 tower29.hadSumEForTime += e29;
0629 }
0630
0631
0632 tower28.E_outer += e28;
0633 tower29.E_outer += e29;
0634 }
0635 }
0636
0637 else {
0638 HcalDetId hcalDetId(detId);
0639
0640
0641
0642 if (hcalDetId.subdet() == HcalOuter) {
0643 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0644 if (towerDetId.null())
0645 return;
0646 MetaTower& tower = find(towerDetId);
0647
0648 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0649 if (theHOIsUsed)
0650 tower.numBadHcalCells += 1;
0651 }
0652
0653 else if (energy >= threshold) {
0654 tower.E_outer += e;
0655
0656 if (theHOIsUsed) {
0657 tower.E += e;
0658
0659 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0660 tower.numRecHcalCells += 1;
0661 } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0662 tower.numProbHcalCells += 1;
0663 }
0664 }
0665
0666
0667 std::pair<DetId, float> mc(detId, e);
0668 tower.metaConstituents.push_back(mc);
0669
0670 }
0671
0672 }
0673
0674
0675 else if (hcalDetId.subdet() == HcalForward) {
0676 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0677 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0678 if (towerDetId.null())
0679 return;
0680 MetaTower& tower = find(towerDetId);
0681 tower.numBadHcalCells += 1;
0682 }
0683
0684 else if (energy >= threshold) {
0685 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0686 if (towerDetId.null())
0687 return;
0688 MetaTower& tower = find(towerDetId);
0689
0690 if (hcalDetId.depth() == 1) {
0691
0692 tower.E_em += e;
0693 } else {
0694
0695 tower.E_em -= e;
0696 tower.E_had += 2. * e;
0697 }
0698 tower.E += e;
0699 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0700 tower.numRecHcalCells += 1;
0701 } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0702 tower.numProbHcalCells += 1;
0703 }
0704
0705
0706
0707 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0708 tower.hadSumTimeTimesE += (e * recHit->time());
0709 tower.hadSumEForTime += e;
0710 }
0711
0712 std::pair<DetId, float> mc(detId, e);
0713 tower.metaConstituents.push_back(mc);
0714
0715 }
0716
0717 }
0718
0719 else {
0720
0721 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0722 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0723 if (towerDetId.null())
0724 return;
0725 MetaTower& tower = find(towerDetId);
0726 tower.numBadHcalCells += 1;
0727 } else if (energy >= threshold) {
0728 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0729 if (towerDetId.null())
0730 return;
0731 MetaTower& tower = find(towerDetId);
0732 tower.E_had += e;
0733 tower.E += e;
0734 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0735 tower.numRecHcalCells += 1;
0736 } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0737 tower.numProbHcalCells += 1;
0738 }
0739
0740
0741
0742 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0743 tower.hadSumTimeTimesE += (e * recHit->time());
0744 tower.hadSumEForTime += e;
0745 }
0746
0747
0748 HcalDetId hcalDetId(detId);
0749 if (hcalDetId.subdet() == HcalEndcap) {
0750 if (theHcalPhase == 0) {
0751 if ((hcalDetId.depth() == 2 && hcalDetId.ietaAbs() >= 18 && hcalDetId.ietaAbs() < 27) ||
0752 (hcalDetId.depth() == 3 && hcalDetId.ietaAbs() == 27) ||
0753 (hcalDetId.depth() == 3 && hcalDetId.ietaAbs() == 16)) {
0754 tower.E_outer += e;
0755 }
0756 }
0757
0758 else if (theHcalPhase == 1) {
0759 if ((hcalDetId.depth() >= 3 && hcalDetId.ietaAbs() >= 18 && hcalDetId.ietaAbs() < 26) ||
0760 (hcalDetId.depth() >= 4 && (hcalDetId.ietaAbs() == 26 || hcalDetId.ietaAbs() == 27)) ||
0761 (hcalDetId.depth() == 3 && hcalDetId.ietaAbs() == 17) ||
0762 (hcalDetId.depth() == 4 && hcalDetId.ietaAbs() == 16)) {
0763 tower.E_outer += e;
0764 }
0765 }
0766 }
0767
0768 std::pair<DetId, float> mc(detId, e);
0769 tower.metaConstituents.push_back(mc);
0770
0771 }
0772
0773 }
0774
0775 }
0776
0777 }
0778
0779 void CaloTowersCreationAlgo::assignHitEcal(const EcalRecHit* recHit) {
0780 DetId detId = recHit->detid();
0781
0782 unsigned int chStatusForCT;
0783 bool ecalIsBad = false;
0784 std::tie(chStatusForCT, ecalIsBad) = ecalChanStatusForCaloTower(recHit);
0785
0786
0787
0788 if (chStatusForCT == CaloTowersCreationAlgo::IgnoredChan)
0789 return;
0790
0791 double threshold, weight;
0792 getThresholdAndWeight(detId, threshold, weight);
0793
0794 double energy = recHit->energy();
0795 double e = energy * weight;
0796
0797
0798
0799
0800
0801
0802
0803
0804 bool passEmThreshold = false;
0805
0806 if (detId.subdetId() == EcalBarrel) {
0807 if (theUseEtEBTresholdFlag)
0808 energy /= cosh((theGeometry->getGeometry(detId)->getPosition()).eta());
0809 if (theUseSymEBTresholdFlag)
0810 passEmThreshold = (fabs(energy) >= threshold);
0811 else
0812 passEmThreshold = (energy >= threshold);
0813
0814 } else if (detId.subdetId() == EcalEndcap) {
0815 if (theUseEtEETresholdFlag)
0816 energy /= cosh((theGeometry->getGeometry(detId)->getPosition()).eta());
0817 if (theUseSymEETresholdFlag)
0818 passEmThreshold = (fabs(energy) >= threshold);
0819 else
0820 passEmThreshold = (energy >= threshold);
0821 }
0822
0823 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0824 if (towerDetId.null())
0825 return;
0826 MetaTower& tower = find(towerDetId);
0827
0828
0829
0830
0831
0832
0833
0834 if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0835 auto thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel(detId);
0836
0837 auto sevit = std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), thisEcalSevLvl);
0838 if (sevit == theEcalSeveritiesToBeExcluded.end())
0839 ++tower.numBadEcalCells;
0840 }
0841
0842
0843 if (chStatusForCT != CaloTowersCreationAlgo::BadChan && passEmThreshold) {
0844 tower.E_em += e;
0845 tower.E += e;
0846
0847 if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0848 tower.numRecEcalCells += 1;
0849 } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0850 tower.numProbEcalCells += 1;
0851 }
0852
0853
0854
0855
0856
0857 if (chStatusForCT == CaloTowersCreationAlgo::GoodChan && e > 0) {
0858 tower.emSumTimeTimesE += (e * recHit->time());
0859 tower.emSumEForTime += e;
0860 }
0861
0862 std::pair<DetId, float> mc(detId, e);
0863 tower.metaConstituents.push_back(mc);
0864 }
0865 }
0866
0867
0868
0869
0870
0871 void CaloTowersCreationAlgo::rescale(const CaloTower* ct) {
0872 double threshold, weight;
0873 CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(ct->id());
0874 if (towerDetId.null())
0875 return;
0876 MetaTower& tower = find(towerDetId);
0877
0878 tower.E_em = 0.;
0879 tower.E_had = 0.;
0880 tower.E_outer = 0.;
0881 for (unsigned int i = 0; i < ct->constituentsSize(); i++) {
0882 DetId detId = ct->constituent(i);
0883 getThresholdAndWeight(detId, threshold, weight);
0884 DetId::Detector det = detId.det();
0885 if (det == DetId::Ecal) {
0886 tower.E_em = ct->emEnergy() * weight;
0887 } else {
0888 HcalDetId hcalDetId(detId);
0889 if (hcalDetId.subdet() == HcalForward) {
0890 if (hcalDetId.depth() == 1)
0891 tower.E_em = ct->emEnergy() * weight;
0892 if (hcalDetId.depth() == 2)
0893 tower.E_had = ct->hadEnergy() * weight;
0894 } else if (hcalDetId.subdet() == HcalOuter) {
0895 tower.E_outer = ct->outerEnergy() * weight;
0896 } else {
0897 tower.E_had = ct->hadEnergy() * weight;
0898 }
0899 }
0900 tower.E = tower.E_had + tower.E_em + tower.E_outer;
0901
0902
0903
0904 std::pair<DetId, float> mc(detId, 0);
0905 tower.metaConstituents.push_back(mc);
0906 }
0907
0908
0909 tower.emSumTimeTimesE = ct->ecalTime();
0910 tower.hadSumTimeTimesE = ct->hcalTime();
0911 tower.emSumEForTime = 1.0;
0912 tower.hadSumEForTime = 1.0;
0913 }
0914
0915 CaloTowersCreationAlgo::MetaTower& CaloTowersCreationAlgo::find(const CaloTowerDetId& detId) {
0916 if (theTowerMap.empty()) {
0917 theTowerMap.resize(theTowerTopology->sizeForDenseIndexing());
0918 }
0919
0920 auto& mt = theTowerMap[theTowerTopology->denseIndex(detId)];
0921
0922 if (mt.empty()) {
0923 mt.id = detId;
0924 mt.metaConstituents.reserve(detId.ietaAbs() < theTowerTopology->firstHFRing() ? 12 : 2);
0925 ++theTowerMapSize;
0926 }
0927
0928 return mt;
0929 }
0930
0931 void CaloTowersCreationAlgo::convert(const CaloTowerDetId& id, const MetaTower& mt, CaloTowerCollection& collection) {
0932 assert(id.rawId() != 0);
0933
0934 double ecalThres = (id.ietaAbs() <= 17) ? (theEBSumThreshold) : (theEESumThreshold);
0935 double E = mt.E;
0936 double E_em = mt.E_em;
0937 double E_had = mt.E_had;
0938 double E_outer = mt.E_outer;
0939
0940
0941
0942
0943
0944
0945
0946
0947 std::vector<std::pair<DetId, float> > metaContains = mt.metaConstituents;
0948 if (id.ietaAbs() < theTowerTopology->firstHFRing() && E_em < ecalThres) {
0949 E -= E_em;
0950 E_em = 0;
0951 std::vector<std::pair<DetId, float> > metaContains_noecal;
0952
0953 for (std::vector<std::pair<DetId, float> >::iterator i = metaContains.begin(); i != metaContains.end(); ++i)
0954 if (i->first.det() != DetId::Ecal)
0955 metaContains_noecal.push_back(*i);
0956 metaContains.swap(metaContains_noecal);
0957 }
0958 if (id.ietaAbs() < theTowerTopology->firstHFRing() && E_had < theHcalThreshold) {
0959 E -= E_had;
0960
0961 if (theHOIsUsed && id.ietaAbs() <= theTowerTopology->lastHORing())
0962 E -= E_outer;
0963
0964 E_had = 0;
0965 E_outer = 0;
0966 std::vector<std::pair<DetId, float> > metaContains_nohcal;
0967
0968 for (std::vector<std::pair<DetId, float> >::iterator i = metaContains.begin(); i != metaContains.end(); ++i)
0969 if (i->first.det() != DetId::Hcal)
0970 metaContains_nohcal.push_back(*i);
0971 metaContains.swap(metaContains_nohcal);
0972 }
0973
0974 if (metaContains.empty())
0975 return;
0976
0977 if (missingHcalRescaleFactorForEcal > 0 && E_had == 0 && E_em > 0) {
0978 auto match = hcalDropChMap.find(id);
0979 if (match != hcalDropChMap.end() && match->second.second) {
0980 E_had = missingHcalRescaleFactorForEcal * E_em;
0981 E += E_had;
0982 }
0983 }
0984
0985 double E_had_tot = (theHOIsUsed && id.ietaAbs() <= theTowerTopology->lastHORing()) ? E_had + E_outer : E_had;
0986
0987
0988
0989 GlobalPoint emPoint, hadPoint;
0990
0991
0992 Basic3DVectorF towerP4;
0993 bool massless = true;
0994
0995 float mass2 = 0;
0996
0997
0998
0999
1000 double momEmDepth = 0.;
1001 double momHadDepth = 0.;
1002 if (id.ietaAbs() <= 17) {
1003 momHadDepth = theMomHBDepth;
1004 momEmDepth = theMomEBDepth;
1005 } else {
1006 momHadDepth = theMomHEDepth;
1007 momEmDepth = theMomEEDepth;
1008 }
1009
1010 switch (theMomConstrMethod) {
1011
1012 case 0: {
1013 GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1014 towerP4 = p.basicVector().unit();
1015 towerP4.mathVector()[3] = 1.f;
1016 towerP4 *= E;
1017
1018
1019
1020
1021 emPoint = p;
1022 hadPoint = p;
1023 }
1024 break;
1025
1026 case 1: {
1027 if (id.ietaAbs() < theTowerTopology->firstHFRing()) {
1028 Basic3DVectorF emP4;
1029 if (E_em > 0) {
1030 emPoint = emShwrPos(metaContains, momEmDepth, E_em);
1031 emP4 = emPoint.basicVector().unit();
1032 emP4.mathVector()[3] = 1.f;
1033 towerP4 = emP4 * E_em;
1034
1035
1036
1037 }
1038 if ((E_had + E_outer) > 0) {
1039 massless = (E_em <= 0);
1040 hadPoint = hadShwrPos(id, momHadDepth);
1041 auto lP4 = hadPoint.basicVector().unit();
1042 lP4.mathVector()[3] = 1.f;
1043 if (!massless) {
1044 auto diff = lP4 - emP4;
1045 mass2 = std::sqrt(E_em * E_had_tot * diff.mag2());
1046 }
1047 lP4 *= E_had_tot;
1048 towerP4 += lP4;
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059 }
1060 } else {
1061 GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1062 towerP4 = p.basicVector().unit();
1063 towerP4.mathVector()[3] = 1.f;
1064 towerP4 *= E;
1065
1066
1067 emPoint = p;
1068 hadPoint = p;
1069 }
1070 }
1071 break;
1072
1073 case 2: {
1074 if (id.ietaAbs() < theTowerTopology->firstHFRing()) {
1075 if (E_em > 0)
1076 emPoint = emShwrLogWeightPos(metaContains, momEmDepth, E_em);
1077 else
1078 emPoint = theTowerGeometry->getGeometry(id)->getPosition();
1079 towerP4 = emPoint.basicVector().unit();
1080 towerP4.mathVector()[3] = 1.f;
1081 towerP4 *= E;
1082
1083
1084
1085
1086 hadPoint = emPoint;
1087 } else {
1088 GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1089 towerP4 = p.basicVector().unit();
1090 towerP4.mathVector()[3] = 1.f;
1091 towerP4 *= E;
1092
1093
1094
1095 emPoint = p;
1096 hadPoint = p;
1097 }
1098 }
1099 break;
1100
1101 }
1102
1103
1104 if UNLIKELY ((towerP4[3] == 0) & (E_outer > 0)) {
1105 float val = theHOIsUsed ? 0 : 1E-9;
1106 collection.emplace_back(id,
1107 E_em,
1108 E_had,
1109 E_outer,
1110 -1,
1111 -1,
1112 CaloTower::PolarLorentzVector(val, hadPoint.eta(), hadPoint.phi(), 0),
1113 emPoint,
1114 hadPoint);
1115 } else {
1116 collection.emplace_back(
1117 id, E_em, E_had, E_outer, -1, -1, GlobalVector(towerP4), towerP4[3], mass2, emPoint, hadPoint);
1118 }
1119 auto& caloTower = collection.back();
1120
1121
1122
1123
1124
1125 if (caloTower.energy() < theEcutTower) {
1126 collection.pop_back();
1127 return;
1128 }
1129
1130
1131 float ecalTime = (mt.emSumEForTime > 0) ? mt.emSumTimeTimesE / mt.emSumEForTime : -9999;
1132 float hcalTime = (mt.hadSumEForTime > 0) ? mt.hadSumTimeTimesE / mt.hadSumEForTime : -9999;
1133 caloTower.setEcalTime(compactTime(ecalTime));
1134 caloTower.setHcalTime(compactTime(hcalTime));
1135
1136 caloTower.setHcalSubdet(theTowerTopology->lastHBRing(),
1137 theTowerTopology->lastHERing(),
1138 theTowerTopology->lastHFRing(),
1139 theTowerTopology->lastHORing());
1140
1141
1142
1143
1144
1145
1146 unsigned int numBadHcalChan = mt.numBadHcalCells;
1147
1148 unsigned int numBadEcalChan = 0;
1149
1150 unsigned int numRecHcalChan = mt.numRecHcalCells;
1151 unsigned int numRecEcalChan = mt.numRecEcalCells;
1152 unsigned int numProbHcalChan = mt.numProbHcalCells;
1153 unsigned int numProbEcalChan = mt.numProbEcalCells;
1154
1155
1156 HcalDropChMap::iterator dropChItr = hcalDropChMap.find(id);
1157 if (dropChItr != hcalDropChMap.end())
1158 numBadHcalChan += dropChItr->second.first;
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210 numBadEcalChan = ecalBadChs[theTowerTopology->denseIndex(id)] + mt.numBadEcalCells;
1211
1212
1213
1214 caloTower.setCaloTowerStatus(
1215 numBadHcalChan, numBadEcalChan, numRecHcalChan, numRecEcalChan, numProbHcalChan, numProbEcalChan);
1216
1217 double maxCellE = -999.0;
1218
1219 std::vector<DetId> contains;
1220 contains.reserve(metaContains.size());
1221 for (std::vector<std::pair<DetId, float> >::iterator i = metaContains.begin(); i != metaContains.end(); ++i) {
1222 contains.push_back(i->first);
1223
1224 if (maxCellE < i->second) {
1225
1226
1227
1228
1229
1230 if (i->first.det() == DetId::Ecal) {
1231 maxCellE = i->second;
1232 } else {
1233 if (HcalDetId(i->first).subdet() != HcalOuter)
1234 maxCellE = i->second;
1235 else if (theHOIsUsed)
1236 maxCellE = i->second;
1237 }
1238
1239 }
1240
1241 }
1242
1243 caloTower.setConstituents(std::move(contains));
1244 caloTower.setHottestCellE(maxCellE);
1245
1246
1247
1248 }
1249
1250 void CaloTowersCreationAlgo::getThresholdAndWeight(const DetId& detId, double& threshold, double& weight) const {
1251 DetId::Detector det = detId.det();
1252 weight = 0;
1253
1254 if (det == DetId::Ecal) {
1255
1256
1257 EcalSubdetector subdet = (EcalSubdetector)(detId.subdetId());
1258 if (subdet == EcalBarrel) {
1259 if (ecalCuts == nullptr) {
1260 threshold = theEBthreshold;
1261 } else {
1262 threshold = (*ecalCuts)[detId];
1263 }
1264 weight = theEBweight;
1265 if (weight <= 0.) {
1266 ROOT::Math::Interpolator my(theEBGrid, theEBWeights, ROOT::Math::Interpolation::kAKIMA);
1267 weight = my.Eval(theEBEScale);
1268 }
1269 } else if (subdet == EcalEndcap) {
1270 if (ecalCuts == nullptr) {
1271 threshold = theEEthreshold;
1272 } else {
1273 threshold = (*ecalCuts)[detId];
1274 }
1275 weight = theEEweight;
1276 if (weight <= 0.) {
1277 ROOT::Math::Interpolator my(theEEGrid, theEEWeights, ROOT::Math::Interpolation::kAKIMA);
1278 weight = my.Eval(theEEEScale);
1279 }
1280 }
1281 } else if (det == DetId::Hcal) {
1282 HcalDetId hcalDetId(detId);
1283 HcalSubdetector subdet = hcalDetId.subdet();
1284 int depth = hcalDetId.depth();
1285
1286 if (subdet == HcalBarrel) {
1287 if (hcalCuts == nullptr) {
1288 threshold = (depth == 1) ? theHBthreshold1 : (depth == 2) ? theHBthreshold2 : theHBthreshold;
1289 } else {
1290 const HcalPFCut* item = hcalCuts->getValues(hcalDetId.rawId());
1291 threshold = item->noiseThreshold();
1292 }
1293 weight = theHBweight;
1294 if (weight <= 0.) {
1295 ROOT::Math::Interpolator my(theHBGrid, theHBWeights, ROOT::Math::Interpolation::kAKIMA);
1296 weight = my.Eval(theHBEScale);
1297 }
1298 }
1299
1300 else if (subdet == HcalEndcap) {
1301
1302 if (hcalDetId.ietaAbs() < theHcalTopology->firstHEDoublePhiRing()) {
1303 if (hcalCuts == nullptr) {
1304 threshold = (depth == 1) ? theHESthreshold1 : theHESthreshold;
1305 } else {
1306 const HcalPFCut* item = hcalCuts->getValues(hcalDetId.rawId());
1307 threshold = item->noiseThreshold();
1308 }
1309 weight = theHESweight;
1310 if (weight <= 0.) {
1311 ROOT::Math::Interpolator my(theHESGrid, theHESWeights, ROOT::Math::Interpolation::kAKIMA);
1312 weight = my.Eval(theHESEScale);
1313 }
1314 } else {
1315 if (hcalCuts == nullptr) {
1316 threshold = (depth == 1) ? theHEDthreshold1 : theHEDthreshold;
1317 } else {
1318 const HcalPFCut* item = hcalCuts->getValues(hcalDetId.rawId());
1319 threshold = item->noiseThreshold();
1320 }
1321 weight = theHEDweight;
1322 if (weight <= 0.) {
1323 ROOT::Math::Interpolator my(theHEDGrid, theHEDWeights, ROOT::Math::Interpolation::kAKIMA);
1324 weight = my.Eval(theHEDEScale);
1325 }
1326 }
1327 }
1328
1329 else if (subdet == HcalOuter) {
1330
1331 if (hcalDetId.ietaAbs() <= 4)
1332 threshold = theHOthreshold0;
1333 else if (hcalDetId.ieta() < 0) {
1334
1335 threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdMinus1 : theHOthresholdMinus2;
1336 } else {
1337
1338 threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdPlus1 : theHOthresholdPlus2;
1339 }
1340 weight = theHOweight;
1341 if (weight <= 0.) {
1342 ROOT::Math::Interpolator my(theHOGrid, theHOWeights, ROOT::Math::Interpolation::kAKIMA);
1343 weight = my.Eval(theHOEScale);
1344 }
1345 }
1346
1347 else if (subdet == HcalForward) {
1348 if (hcalDetId.depth() == 1) {
1349 threshold = theHF1threshold;
1350 weight = theHF1weight;
1351 if (weight <= 0.) {
1352 ROOT::Math::Interpolator my(theHF1Grid, theHF1Weights, ROOT::Math::Interpolation::kAKIMA);
1353 weight = my.Eval(theHF1EScale);
1354 }
1355 } else {
1356 threshold = theHF2threshold;
1357 weight = theHF2weight;
1358 if (weight <= 0.) {
1359 ROOT::Math::Interpolator my(theHF2Grid, theHF2Weights, ROOT::Math::Interpolation::kAKIMA);
1360 weight = my.Eval(theHF2EScale);
1361 }
1362 }
1363 }
1364 } else {
1365 edm::LogError("CaloTowersCreationAlgo") << "Bad cell: " << det << std::endl;
1366 }
1367 }
1368
1369 void CaloTowersCreationAlgo::setEBEScale(double scale) {
1370 if (scale > 0.00001)
1371 *&theEBEScale = scale;
1372 else
1373 *&theEBEScale = 50.;
1374 }
1375
1376 void CaloTowersCreationAlgo::setEEEScale(double scale) {
1377 if (scale > 0.00001)
1378 *&theEEEScale = scale;
1379 else
1380 *&theEEEScale = 50.;
1381 }
1382
1383 void CaloTowersCreationAlgo::setHBEScale(double scale) {
1384 if (scale > 0.00001)
1385 *&theHBEScale = scale;
1386 else
1387 *&theHBEScale = 50.;
1388 }
1389
1390 void CaloTowersCreationAlgo::setHESEScale(double scale) {
1391 if (scale > 0.00001)
1392 *&theHESEScale = scale;
1393 else
1394 *&theHESEScale = 50.;
1395 }
1396
1397 void CaloTowersCreationAlgo::setHEDEScale(double scale) {
1398 if (scale > 0.00001)
1399 *&theHEDEScale = scale;
1400 else
1401 *&theHEDEScale = 50.;
1402 }
1403
1404 void CaloTowersCreationAlgo::setHOEScale(double scale) {
1405 if (scale > 0.00001)
1406 *&theHOEScale = scale;
1407 else
1408 *&theHOEScale = 50.;
1409 }
1410
1411 void CaloTowersCreationAlgo::setHF1EScale(double scale) {
1412 if (scale > 0.00001)
1413 *&theHF1EScale = scale;
1414 else
1415 *&theHF1EScale = 50.;
1416 }
1417
1418 void CaloTowersCreationAlgo::setHF2EScale(double scale) {
1419 if (scale > 0.00001)
1420 *&theHF2EScale = scale;
1421 else
1422 *&theHF2EScale = 50.;
1423 }
1424
1425 GlobalPoint CaloTowersCreationAlgo::emCrystalShwrPos(DetId detId, float fracDepth) {
1426 auto cellGeometry = theGeometry->getGeometry(detId);
1427 GlobalPoint point = cellGeometry->getPosition();
1428
1429 if (fracDepth <= 0)
1430 return point;
1431 if (fracDepth > 1)
1432 fracDepth = 1;
1433
1434 const GlobalPoint& backPoint = cellGeometry->getBackPoint();
1435 point += fracDepth * (backPoint - point);
1436
1437 return point;
1438 }
1439
1440 GlobalPoint CaloTowersCreationAlgo::hadSegmentShwrPos(DetId detId, float fracDepth) {
1441
1442 return emCrystalShwrPos(detId, fracDepth);
1443 }
1444
1445 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(const std::vector<std::pair<DetId, float> >& metaContains,
1446 float fracDepth,
1447 double hadE) {
1448
1449
1450 #ifdef EDM_ML_DEBUG
1451 std::cout << "hadShwrPos called with " << metaContains.size() << " elements and energy " << hadE << ":" << fracDepth
1452 << std::endl;
1453 #endif
1454 if (hadE <= 0)
1455 return GlobalPoint(0, 0, 0);
1456
1457 double hadX = 0.0;
1458 double hadY = 0.0;
1459 double hadZ = 0.0;
1460
1461 int nConst = 0;
1462
1463 std::vector<std::pair<DetId, float> >::const_iterator mc_it = metaContains.begin();
1464 for (; mc_it != metaContains.end(); ++mc_it) {
1465 if (mc_it->first.det() != DetId::Hcal)
1466 continue;
1467
1468 if (HcalDetId(mc_it->first).subdet() == HcalOuter)
1469 continue;
1470 ++nConst;
1471
1472 GlobalPoint p = hadSegmentShwrPos(mc_it->first, fracDepth);
1473
1474
1475
1476 hadX += p.x();
1477 hadY += p.y();
1478 hadZ += p.z();
1479 }
1480
1481 return GlobalPoint(hadX / nConst, hadY / nConst, hadZ / nConst);
1482 }
1483
1484 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(CaloTowerDetId towerId, float fracDepth) {
1485
1486
1487
1488
1489 #ifdef EDM_ML_DEBUG
1490 std::cout << "hadShwrPos " << towerId << " frac " << fracDepth << std::endl;
1491 #endif
1492 if (fracDepth < 0)
1493 fracDepth = 0;
1494 else if (fracDepth > 1)
1495 fracDepth = 1;
1496
1497 GlobalPoint point(0, 0, 0);
1498
1499 int iEta = towerId.ieta();
1500 int iPhi = towerId.iphi();
1501
1502 HcalDetId frontCellId, backCellId;
1503
1504 if (towerId.ietaAbs() >= theTowerTopology->firstHFRing()) {
1505
1506 frontCellId = HcalDetId(HcalForward, towerId.zside() * theTowerTopology->convertCTtoHcal(abs(iEta)), iPhi, 1);
1507 backCellId = HcalDetId(HcalForward, towerId.zside() * theTowerTopology->convertCTtoHcal(abs(iEta)), iPhi, 1);
1508 } else {
1509
1510 std::vector<DetId> items = theTowerConstituentsMap->constituentsOf(towerId);
1511 int frontDepth = 1000;
1512 int backDepth = -1000;
1513 for (unsigned i = 0; i < items.size(); i++) {
1514 if (items[i].det() != DetId::Hcal)
1515 continue;
1516 HcalDetId hid(items[i]);
1517 if (hid.subdet() == HcalOuter)
1518 continue;
1519 if (!theHcalTopology->validHcal(hid, 2))
1520 continue;
1521
1522 if (theHcalTopology->idFront(hid).depth() < frontDepth) {
1523 frontCellId = hid;
1524 frontDepth = theHcalTopology->idFront(hid).depth();
1525 }
1526 if (theHcalTopology->idBack(hid).depth() > backDepth) {
1527 backCellId = hid;
1528 backDepth = theHcalTopology->idBack(hid).depth();
1529 }
1530 }
1531 #ifdef EDM_ML_DEBUG
1532 std::cout << "Front " << frontCellId << " Back " << backCellId << " Depths " << frontDepth << ":" << backDepth
1533 << std::endl;
1534 #endif
1535
1536 if (towerId.ietaAbs() == theTowerTopology->lastHERing() && (theHcalPhase == 0 || theHcalPhase == 1)) {
1537 CaloTowerDetId towerId28(towerId.ieta() - towerId.zside(), towerId.iphi());
1538 std::vector<DetId> items28 = theTowerConstituentsMap->constituentsOf(towerId28);
1539 #ifdef EDM_ML_DEBUG
1540 std::cout << towerId28 << " with " << items28.size() << " constituents:";
1541 for (unsigned k = 0; k < items28.size(); ++k)
1542 if (items28[k].det() == DetId::Hcal)
1543 std::cout << " " << HcalDetId(items28[k]);
1544 std::cout << std::endl;
1545 #endif
1546 for (unsigned i = 0; i < items28.size(); i++) {
1547 if (items28[i].det() != DetId::Hcal)
1548 continue;
1549 HcalDetId hid(items28[i]);
1550 if (hid.subdet() == HcalOuter)
1551 continue;
1552
1553 if (theHcalTopology->idBack(hid).depth() > backDepth) {
1554 backCellId = hid;
1555 backDepth = theHcalTopology->idBack(hid).depth();
1556 }
1557 }
1558 }
1559 #ifdef EDM_ML_DEBUG
1560 std::cout << "Back " << backDepth << " ID " << backCellId << std::endl;
1561 #endif
1562 }
1563 point = hadShwPosFromCells(DetId(frontCellId), DetId(backCellId), fracDepth);
1564
1565 return point;
1566 }
1567
1568 GlobalPoint CaloTowersCreationAlgo::hadShwPosFromCells(DetId frontCellId, DetId backCellId, float fracDepth) {
1569
1570
1571
1572 HcalDetId hid1(frontCellId), hid2(backCellId);
1573 if (theHcalTopology->getMergePositionFlag()) {
1574 hid1 = theHcalTopology->idFront(frontCellId);
1575 #ifdef EDM_ML_DEBUG
1576 std::cout << "Front " << HcalDetId(frontCellId) << " " << hid1 << "\n";
1577 #endif
1578 hid2 = theHcalTopology->idBack(backCellId);
1579 #ifdef EDM_ML_DEBUG
1580 std::cout << "Back " << HcalDetId(backCellId) << " " << hid2 << "\n";
1581 #endif
1582 }
1583
1584 auto frontCellGeometry = theGeometry->getGeometry(DetId(hid1));
1585 auto backCellGeometry = theGeometry->getGeometry(DetId(hid2));
1586
1587 GlobalPoint point = frontCellGeometry->getPosition();
1588 const GlobalPoint& backPoint = backCellGeometry->getBackPoint();
1589
1590 point += fracDepth * (backPoint - point);
1591
1592 return point;
1593 }
1594
1595 GlobalPoint CaloTowersCreationAlgo::emShwrPos(const std::vector<std::pair<DetId, float> >& metaContains,
1596 float fracDepth,
1597 double emE) {
1598 if (emE <= 0)
1599 return GlobalPoint(0, 0, 0);
1600
1601 double emX = 0.0;
1602 double emY = 0.0;
1603 double emZ = 0.0;
1604
1605 double eSum = 0;
1606
1607 std::vector<std::pair<DetId, float> >::const_iterator mc_it = metaContains.begin();
1608 for (; mc_it != metaContains.end(); ++mc_it) {
1609 if (mc_it->first.det() != DetId::Ecal)
1610 continue;
1611 GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
1612 double e = mc_it->second;
1613
1614 if (e > 0) {
1615 emX += p.x() * e;
1616 emY += p.y() * e;
1617 emZ += p.z() * e;
1618 eSum += e;
1619 }
1620 }
1621
1622 return GlobalPoint(emX / eSum, emY / eSum, emZ / eSum);
1623 }
1624
1625 GlobalPoint CaloTowersCreationAlgo::emShwrLogWeightPos(const std::vector<std::pair<DetId, float> >& metaContains,
1626 float fracDepth,
1627 double emE) {
1628 double emX = 0.0;
1629 double emY = 0.0;
1630 double emZ = 0.0;
1631
1632 double weight = 0;
1633 double sumWeights = 0;
1634 double sumEmE = 0;
1635 double crystalThresh = 0.015 * emE;
1636
1637 std::vector<std::pair<DetId, float> >::const_iterator mc_it = metaContains.begin();
1638 for (; mc_it != metaContains.end(); ++mc_it) {
1639 if (mc_it->second < 0)
1640 continue;
1641 if (mc_it->first.det() == DetId::Ecal && mc_it->second > crystalThresh)
1642 sumEmE += mc_it->second;
1643 }
1644
1645 for (mc_it = metaContains.begin(); mc_it != metaContains.end(); ++mc_it) {
1646 if (mc_it->first.det() != DetId::Ecal || mc_it->second < crystalThresh)
1647 continue;
1648
1649 GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
1650
1651 weight = 4.2 + log(mc_it->second / sumEmE);
1652 sumWeights += weight;
1653
1654 emX += p.x() * weight;
1655 emY += p.y() * weight;
1656 emZ += p.z() * weight;
1657 }
1658
1659 return GlobalPoint(emX / sumWeights, emY / sumWeights, emZ / sumWeights);
1660 }
1661
1662 int CaloTowersCreationAlgo::compactTime(float time) {
1663 const float timeUnit = 0.01;
1664
1665 if (time > 300.0)
1666 return 30000;
1667 if (time < -300.0)
1668 return -30000;
1669
1670 return int(time / timeUnit + 0.5);
1671 }
1672
1673
1674
1675
1676
1677 void CaloTowersCreationAlgo::makeHcalDropChMap() {
1678
1679
1680
1681 hcalDropChMap.clear();
1682 std::vector<DetId> allChanInStatusCont = theHcalChStatus->getAllChannels();
1683
1684 #ifdef EDM_ML_DEBUG
1685 std::cout << "DropChMap with " << allChanInStatusCont.size() << " channels" << std::endl;
1686 #endif
1687 for (std::vector<DetId>::iterator it = allChanInStatusCont.begin(); it != allChanInStatusCont.end(); ++it) {
1688 const uint32_t dbStatusFlag = theHcalChStatus->getValues(*it)->getValue();
1689 if (theHcalSevLvlComputer->dropChannel(dbStatusFlag)) {
1690 DetId id = theHcalTopology->mergedDepthDetId(HcalDetId(*it));
1691
1692 CaloTowerDetId twrId = theTowerConstituentsMap->towerOf(id);
1693
1694 hcalDropChMap[twrId].first += 1;
1695
1696 HcalDetId hid(*it);
1697
1698
1699 if (hid.subdet() == HcalEndcap && (theHcalPhase == 0 || theHcalPhase == 1) &&
1700 hid.ietaAbs() == theHcalTopology->lastHERing() - 1) {
1701 bool merge = theHcalTopology->mergedDepth29(hid);
1702 if (merge) {
1703 CaloTowerDetId twrId29(twrId.ieta() + twrId.zside(), twrId.iphi());
1704 hcalDropChMap[twrId29].first += 1;
1705 }
1706 }
1707 }
1708 }
1709
1710 if (missingHcalRescaleFactorForEcal > 0) {
1711 for (auto& pair : hcalDropChMap) {
1712 if (pair.second.first == 0)
1713 continue;
1714 int ngood = 0, nbad = 0;
1715 for (DetId id : theTowerConstituentsMap->constituentsOf(pair.first)) {
1716 if (id.det() != DetId::Hcal)
1717 continue;
1718 HcalDetId hid(id);
1719 if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap)
1720 continue;
1721 const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1722 if (dbStatusFlag == 0 || !theHcalSevLvlComputer->dropChannel(dbStatusFlag)) {
1723 ngood += 1;
1724 } else {
1725 nbad += 1;
1726 }
1727 }
1728 if (nbad > 0 && nbad >= ngood) {
1729
1730
1731
1732 pair.second.second = true;
1733 }
1734 }
1735 }
1736 }
1737
1738 void CaloTowersCreationAlgo::makeEcalBadChs() {
1739
1740
1741
1742
1743 for (auto ind = 0U; ind < theTowerTopology->sizeForDenseIndexing(); ++ind) {
1744 auto& numBadEcalChan = ecalBadChs[ind];
1745 numBadEcalChan = 0;
1746 auto id = theTowerTopology->detIdFromDenseIndex(ind);
1747
1748
1749
1750
1751 std::vector<DetId> allConstituents = theTowerConstituentsMap->constituentsOf(id);
1752
1753 for (std::vector<DetId>::iterator ac_it = allConstituents.begin(); ac_it != allConstituents.end(); ++ac_it) {
1754 if (ac_it->det() != DetId::Ecal)
1755 continue;
1756
1757 auto thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel(*ac_it);
1758
1759
1760 std::vector<int>::const_iterator sevit =
1761 std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), thisEcalSevLvl);
1762 if (sevit != theEcalSeveritiesToBeExcluded.end()) {
1763 ++numBadEcalChan;
1764 }
1765 }
1766
1767
1768 }
1769
1770
1771
1772
1773
1774
1775
1776
1777 }
1778
1779
1780
1781 unsigned int CaloTowersCreationAlgo::hcalChanStatusForCaloTower(const CaloRecHit* hit) {
1782 HcalDetId hid(hit->detid());
1783 DetId id = theHcalTopology->idFront(hid);
1784 #ifdef EDM_ML_DEBUG
1785 std::cout << "ChanStatusForCaloTower for " << hid << " to " << HcalDetId(id) << std::endl;
1786 #endif
1787 const uint32_t recHitFlag = hit->flags();
1788 const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1789
1790 int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1791 bool isRecovered = theHcalSevLvlComputer->recoveredRecHit(id, recHitFlag);
1792
1793
1794 if (useRejectedHitsOnly) {
1795 if (!isRecovered) {
1796 if (severityLevel <= int(theHcalAcceptSeverityLevel) ||
1797 severityLevel > int(theHcalAcceptSeverityLevelForRejectedHit))
1798 return CaloTowersCreationAlgo::IgnoredChan;
1799
1800 } else {
1801 if (theRecoveredHcalHitsAreUsed || !useRejectedRecoveredHcalHits) {
1802
1803 return CaloTowersCreationAlgo::IgnoredChan;
1804 } else if (useRejectedRecoveredHcalHits) {
1805 return CaloTowersCreationAlgo::RecoveredChan;
1806 }
1807
1808 }
1809
1810
1811
1812 return CaloTowersCreationAlgo::ProblematicChan;
1813
1814 }
1815
1816
1817
1818 if (severityLevel == 0)
1819 return CaloTowersCreationAlgo::GoodChan;
1820
1821 if (isRecovered) {
1822 return (theRecoveredHcalHitsAreUsed) ? CaloTowersCreationAlgo::RecoveredChan : CaloTowersCreationAlgo::BadChan;
1823 } else {
1824 if (severityLevel > int(theHcalAcceptSeverityLevel)) {
1825 return CaloTowersCreationAlgo::BadChan;
1826 } else {
1827 return CaloTowersCreationAlgo::ProblematicChan;
1828 }
1829 }
1830 }
1831
1832 std::tuple<unsigned int, bool> CaloTowersCreationAlgo::ecalChanStatusForCaloTower(const EcalRecHit* hit) {
1833
1834
1835
1836
1837
1838
1839
1840
1841
1842
1843 EcalRecHit const& rh = *reinterpret_cast<EcalRecHit const*>(hit);
1844 int severityLevel = theEcalSevLvlAlgo->severityLevel(rh);
1845
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860 bool isBad = (severityLevel == EcalSeverityLevel::kBad);
1861
1862 bool isRecovered = (severityLevel == EcalSeverityLevel::kRecovered);
1863
1864
1865
1866 std::vector<int>::const_iterator sevit =
1867 std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), severityLevel);
1868 bool accepted = (sevit == theEcalSeveritiesToBeExcluded.end());
1869
1870
1871
1872
1873
1874
1875 if (useRejectedHitsOnly) {
1876 if (!isRecovered) {
1877 if (accepted || std::find(theEcalSeveritiesToBeUsedInBadTowers.begin(),
1878 theEcalSeveritiesToBeUsedInBadTowers.end(),
1879 severityLevel) == theEcalSeveritiesToBeUsedInBadTowers.end())
1880 return std::make_tuple(CaloTowersCreationAlgo::IgnoredChan, isBad);
1881
1882 } else {
1883 if (theRecoveredEcalHitsAreUsed || !useRejectedRecoveredEcalHits) {
1884
1885 return std::make_tuple(CaloTowersCreationAlgo::IgnoredChan, isBad);
1886 ;
1887 } else if (useRejectedRecoveredEcalHits) {
1888 return std::make_tuple(CaloTowersCreationAlgo::RecoveredChan, isBad);
1889 }
1890
1891 }
1892
1893
1894 return std::make_tuple(CaloTowersCreationAlgo::ProblematicChan, isBad);
1895
1896 }
1897
1898
1899 if (severityLevel == EcalSeverityLevel::kGood)
1900 return std::make_tuple(CaloTowersCreationAlgo::GoodChan, false);
1901
1902 if (isRecovered) {
1903 return std::make_tuple(
1904 (theRecoveredEcalHitsAreUsed) ? CaloTowersCreationAlgo::RecoveredChan : CaloTowersCreationAlgo::BadChan, true);
1905 } else {
1906 return std::make_tuple(accepted ? CaloTowersCreationAlgo::ProblematicChan : CaloTowersCreationAlgo::BadChan, isBad);
1907 }
1908 }