Back to home page

Project CMSSW displayed by LXR

 
 

    


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 //#define EDM_ML_DEBUG
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       // (for momentum reconstruction algorithm)
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                                                // (momentum reconstruction algorithm)
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       // (momentum reconstruction algorithm)
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                                                // (for the momentum construction algorithm)
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       // (momentum reconstruction algorithm)
0345       theMomConstrMethod(momConstrMethod),
0346       theMomHBDepth(momHBDepth),
0347       theMomHEDepth(momHEDepth),
0348       theMomEBDepth(momEBDepth),
0349       theMomEEDepth(momEEDepth),
0350       theHcalPhase(hcalPhase) {
0351   // static int N = 0;
0352   // std::cout << "VI Algo " << ++N << std::endl;
0353   // nalgo=N;
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   //initialize ecal bad channel map
0372   ecalBadChs.resize(theTowerTopology->sizeForDenseIndexing(), 0);
0373 }
0374 
0375 void CaloTowersCreationAlgo::begin() {
0376   theTowerMap.clear();
0377   theTowerMapSize = 0;
0378   //hcalDropChMap.clear();
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 // this method should not be used any more as the towers in the changed format
0402 // can not be properly rescaled with the "rescale" method.
0403 // "rescale was replaced by "rescaleTowers"
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   // now copy this map into the final collection
0413   result.reserve(theTowerMapSize);
0414   // auto k=0U;
0415   //  if (!theEbHandle.isValid()) std::cout << "VI ebHandle not valid" << std::endl;
0416   // if (!theEeHandle.isValid()) std::cout << "VI eeHandle not valid" << std::endl;
0417 
0418   for (auto const& mt : theTowerMap) {
0419     // Convert only if there is at least one constituent in the metatower.
0420     // The check of constituents size in the coverted tower is still needed!
0421     if (!mt.empty()) {
0422       convert(mt.id, mt, result);
0423     }  // ++k;}
0424   }
0425 
0426   // assert(k==theTowerMapSize);
0427   // std::cout << "VI TowerMap " << theTowerMapSize << " " << k << std::endl;
0428 
0429   theTowerMap.clear();  // save the memory
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;  // not used: we do not change thresholds
0441     double weight = 1.0;
0442 
0443     // HF
0444     if (ctcItr->ietaAbs() >= theTowerTopology->firstHFRing()) {
0445       double E_short = 0.5 * newE_had;           // from the definitions for HF
0446       double E_long = newE_em + 0.5 * newE_had;  //
0447       // scale
0448       E_long *= theHF1weight;
0449       E_short *= theHF2weight;
0450       // convert
0451       newE_em = E_long - E_short;
0452       newE_had = 2.0 * E_short;
0453     }
0454 
0455     else {  // barrel/endcap
0456 
0457       // find if its in EB, or EE; determine from first ecal constituent found
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       // HO
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       // HB/HE
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     }  // barrel/endcap region
0492 
0493     // now make the new tower
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       // for HF we use common point for ecal, hcal shower positions regardless of the method
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     // copy the timings, have to convert back to int, 1 unit = 0.01 ns
0520     rescaledTower.setEcalTime(int(ctcItr->ecalTime() * 100.0 + 0.5));
0521     rescaledTower.setHcalTime(int(ctcItr->hcalTime() * 100.0 + 0.5));
0522     //add topology info
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   }  // end of loop over towers
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   // this is for skipping channls: mostly needed for the creation of
0554   // bad towers from hits i the bad channel collections.
0555   if (chStatusForCT == CaloTowersCreationAlgo::IgnoredChan)
0556     return;
0557 
0558   double threshold, weight;
0559   getThresholdAndWeight(detId, threshold, weight);
0560 
0561   double energy = recHit->energy();  // original RecHit energy is used to apply thresholds
0562   double e = energy * weight;        // energies scaled by user weight: used in energy assignments
0563 
0564   // SPECIAL handling of tower 28 merged depths --> half into tower 28 and half into tower 29
0565   bool merge(false);
0566   if (detIdF.det() == DetId::Hcal && HcalDetId(detIdF).subdet() == HcalEndcap &&
0567       (theHcalPhase == 0 || theHcalPhase == 1) &&
0568       //HcalDetId(detId).depth()==3 &&
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     //////////////////////////////    unsigned int chStatusForCT = hcalChanStatusForCaloTower(recHit);
0577 
0578     // bad channels are counted regardless of energy threshold
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) {  // not bad channel: use energy if above 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       // NOTE DIVIDE BY 2!!!
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       // time info: do not use in averaging if timing error is found: need
0622       // full set of status info to implement: use only "good" channels for now
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       // store the energy in layer 3 also in E_outer
0632       tower28.E_outer += e28;
0633       tower29.E_outer += e29;
0634     }  // not a "bad" hit
0635   }  // end of special case
0636 
0637   else {
0638     HcalDetId hcalDetId(detId);
0639 
0640     ///////////////////////      unsigned int chStatusForCT = hcalChanStatusForCaloTower(recHit);
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;  // store HO energy even if HO is not used
0655         // add energy of the tower and/or flag if theHOIsUsed
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         }  // HO is used
0665 
0666         // add HO to constituents even if it is not used: JetMET wants to keep these towers
0667         std::pair<DetId, float> mc(detId, e);
0668         tower.metaConstituents.push_back(mc);
0669 
0670       }  // not a bad channel, energy above threshold
0671 
0672     }  // HO hit
0673 
0674     // HF calculates EM fraction differently
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           // long fiber, so E_EM = E(Long) - E(Short)
0692           tower.E_em += e;
0693         } else {
0694           // short fiber, EHAD = 2 * E(Short)
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         // put the timing in HCAL -> have to check timing errors when available
0706         // for now use only good channels
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       }  // not a bad HF channel, energy above threshold
0716 
0717     }  // HF hit
0718 
0719     else {
0720       // HCAL situation normal in HB/HE
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         // Timing information: need specific accessors
0741         // for now use only good channels
0742         if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0743           tower.hadSumTimeTimesE += (e * recHit->time());
0744           tower.hadSumEForTime += e;
0745         }
0746         // store energy in highest depth for towers 18-27 (for electron,photon ID in endcap)
0747         // also, store energy in HE part of tower 16 (for JetMET cleanup)
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           //combine depths in phase0-like way
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       }  // not a "bad" channel, energy above threshold
0772 
0773     }  // channel in HBHE (excluding twrs 28,29)
0774 
0775   }  // recHit normal case (not in HE towers 28,29)
0776 
0777 }  // end of assignHitHcal method
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   // this is for skipping channls: mostly needed for the creation of
0787   // bad towers from hits i the bad channel collections.
0788   if (chStatusForCT == CaloTowersCreationAlgo::IgnoredChan)
0789     return;
0790 
0791   double threshold, weight;
0792   getThresholdAndWeight(detId, threshold, weight);
0793 
0794   double energy = recHit->energy();  // original RecHit energy is used to apply thresholds
0795   double e = energy * weight;        // energies scaled by user weight: used in energy assignments
0796 
0797   /////////////////////////////      unsigned int chStatusForCT = ecalChanStatusForCaloTower(recHit);
0798 
0799   // For ECAL we count all bad channels after the metatower is complete
0800 
0801   // Include options for symmetric thresholds and cut on Et
0802   // for ECAL RecHits
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   // count bad cells and avoid double counting with those from DB (Recovered are counted bad)
0829 
0830   // somehow misses some
0831   // if ( (chStatusForCT == CaloTowersCreationAlgo::BadChan) & (!ecalIsBad) ) ++tower.numBadEcalCells;
0832 
0833   // a bit slower...
0834   if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0835     auto thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel(detId);
0836     // check if the Ecal severity is ok to keep
0837     auto sevit = std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), thisEcalSevLvl);
0838     if (sevit == theEcalSeveritiesToBeExcluded.end())
0839       ++tower.numBadEcalCells;  // notinDB
0840   }
0841 
0842   //      if (chStatusForCT != CaloTowersCreationAlgo::BadChan && energy >= threshold) {
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     // change when full status info is available
0854     // for now use only good channels
0855 
0856     // add e>0 check (new options allow e<0)
0857     if (chStatusForCT == CaloTowersCreationAlgo::GoodChan && e > 0) {
0858       tower.emSumTimeTimesE += (e * recHit->time());
0859       tower.emSumEForTime += e;  // see above
0860     }
0861 
0862     std::pair<DetId, float> mc(detId, e);
0863     tower.metaConstituents.push_back(mc);
0864   }
0865 }  // end of assignHitEcal method
0866 
0867 // This method is not flexible enough for the new CaloTower format.
0868 // For now make a quick compatibility "fix" : WILL NOT WORK CORRECTLY with anything
0869 // except the default simple p4 assignment!!!
0870 // Must be rewritten for full functionality.
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     // this is to be compliant with the new MetaTower setup
0903     // used only for the default simple vector assignment
0904     std::pair<DetId, float> mc(detId, 0);
0905     tower.metaConstituents.push_back(mc);
0906   }
0907 
0908   // preserve time inforamtion
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   // Note: E_outer is used to save HO energy OR energy in the outermost depths in endcap region
0941   // In the methods with separate treatment of EM and HAD components:
0942   //  - HO is not used to determine direction, however HO energy is added to get "total had energy"
0943   //  => Check if the tower is within HO coverage before adding E_outer to the "total had" energy
0944   //     else the energy will be double counted
0945   // When summing up the energy of the tower these checks are performed in the loops over RecHits
0946 
0947   std::vector<std::pair<DetId, float> > metaContains = mt.metaConstituents;
0948   if (id.ietaAbs() < theTowerTopology->firstHFRing() && E_em < ecalThres) {  // ignore EM threshold in HF
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;  // not subtracted before, think it should be done
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   // create CaloTower using the selected algorithm
0988 
0989   GlobalPoint emPoint, hadPoint;
0990 
0991   // this is actually a 4D vector
0992   Basic3DVectorF towerP4;
0993   bool massless = true;
0994   // float mass1=0;
0995   float mass2 = 0;
0996 
0997   // conditional assignment of depths for barrel/endcap
0998   // Some additional tuning may be required in the transitional region
0999   // 14<|iEta|<19
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       // FIXME  : move to simple cartesian algebra
1012     case 0: {  // Simple 4-momentum assignment
1013       GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1014       towerP4 = p.basicVector().unit();
1015       towerP4.mathVector()[3] = 1.f;  // energy
1016       towerP4 *= E;
1017 
1018       // double pf=1.0/cosh(p.eta());
1019       // if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);
1020 
1021       emPoint = p;
1022       hadPoint = p;
1023     }  // end case 0
1024     break;
1025 
1026     case 1: {  // separate 4-vectors for ECAL, HCAL, add to get the 4-vector of the tower (=>tower has mass!)
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;  // energy
1033           towerP4 = emP4 * E_em;
1034 
1035           // double emPf = 1.0/cosh(emPoint.eta());
1036           // towerP4 += CaloTower::PolarLorentzVector(E_em*emPf, emPoint.eta(), emPoint.phi(), 0);
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;  // energy
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             if (!massless) {
1051                auto p = towerP4;
1052                double m2 = double(p[3]*p[3]) - double(p[0]*p[0])+double(p[1]*p[1])+double(p[2]*p[2]); mass1 = m2>0 ? std::sqrt(m2) : 0;
1053             }       
1054             */
1055           // double hadPf = 1.0/cosh(hadPoint.eta());
1056           // if (E_had_tot>0) {
1057           //  towerP4 += CaloTower::PolarLorentzVector(E_had_tot*hadPf, hadPoint.eta(), hadPoint.phi(), 0);
1058           // }
1059         }
1060       } else {  // forward detector: use the CaloTower position
1061         GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1062         towerP4 = p.basicVector().unit();
1063         towerP4.mathVector()[3] = 1.f;  // energy
1064         towerP4 *= E;
1065         // double pf=1.0/cosh(p.eta());
1066         // if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);  // simple momentum assignment, same position
1067         emPoint = p;
1068         hadPoint = p;
1069       }
1070     }  // end case 1
1071     break;
1072 
1073     case 2: {  // use ECAL position for the tower (when E_cal>0), else default CaloTower position (massless tower)
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;  // energy
1081         towerP4 *= E;
1082 
1083         // double sumPf = 1.0/cosh(emPoint.eta());
1084         /// if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*sumPf, emPoint.eta(), emPoint.phi(), 0);
1085 
1086         hadPoint = emPoint;
1087       } else {  // forward detector: use the CaloTower position
1088         GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1089         towerP4 = p.basicVector().unit();
1090         towerP4.mathVector()[3] = 1.f;  // energy
1091         towerP4 *= E;
1092 
1093         // double pf=1.0/cosh(p.eta());
1094         // if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);  // simple momentum assignment, same position
1095         emPoint = p;
1096         hadPoint = p;
1097       }
1098     }  // end case 2
1099     break;
1100 
1101   }  // end of decision on p4 reconstruction method
1102 
1103   // insert in collection (remove and return if below threshold)
1104   if UNLIKELY ((towerP4[3] == 0) & (E_outer > 0)) {
1105     float val = theHOIsUsed ? 0 : 1E-9;  // to keep backwards compatibility for theHOIsUsed == true
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   // if (!massless) std::cout << "massive " << id <<' ' << mass1 <<' ' << mass2 <<' ' <<  caloTower.mass() << std::endl;
1122   // std::cout << "CaloTowerVI " <<theMomConstrMethod <<' ' << id <<' '<< E_em <<' '<< E_had <<' '<< E_outer <<' '<< GlobalVector(towerP4) <<' '<< towerP4[3] <<' '<< emPoint <<' '<< hadPoint << std::endl;
1123   //if (towerP4[3]==0) std::cout << "CaloTowerVIzero " << theEcutTower << ' ' << collection.back().eta() <<' '<< collection.back().phi() << std::endl;
1124 
1125   if (caloTower.energy() < theEcutTower) {
1126     collection.pop_back();
1127     return;
1128   }
1129 
1130   // set the timings
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   //add topology info
1136   caloTower.setHcalSubdet(theTowerTopology->lastHBRing(),
1137                           theTowerTopology->lastHERing(),
1138                           theTowerTopology->lastHFRing(),
1139                           theTowerTopology->lastHORing());
1140 
1141   // set the CaloTower status word =====================================
1142   // Channels must be counter exclusively in the defined cathegories
1143   // "Bad" channels (not used in energy assignment) can be flagged during
1144   // CaloTower creation only if specified in the configuration file
1145 
1146   unsigned int numBadHcalChan = mt.numBadHcalCells;
1147   //    unsigned int numBadEcalChan  = mt.numBadEcalCells;
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   // now add dead/off/... channels not used in RecHit reconstruction for HCAL
1156   HcalDropChMap::iterator dropChItr = hcalDropChMap.find(id);
1157   if (dropChItr != hcalDropChMap.end())
1158     numBadHcalChan += dropChItr->second.first;
1159 
1160   // for ECAL the number of all bad channels is obtained here -----------------------
1161 
1162   /*
1163     // old hyper slow algorithm    
1164     // get all possible constituents of the tower
1165     std::vector<DetId> allConstituents = theTowerConstituentsMap->constituentsOf(id);
1166 
1167     for (std::vector<DetId>::iterator ac_it=allConstituents.begin(); 
1168      ac_it!=allConstituents.end(); ++ac_it) {
1169 
1170       if (ac_it->det()!=DetId::Ecal) continue;
1171  
1172       int thisEcalSevLvl = -999;
1173      
1174       if (ac_it->subdetId() == EcalBarrel && theEbHandle.isValid()) {
1175     thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEbHandle);//, *theEcalChStatus);
1176       }
1177       else if (ac_it->subdetId() == EcalEndcap && theEeHandle.isValid()) {
1178     thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEeHandle);//, *theEcalChStatus);
1179       }
1180  
1181       // check if the Ecal severity is ok to keep
1182       std::vector<int>::const_iterator sevit = std::find(theEcalSeveritiesToBeExcluded.begin(),
1183                              theEcalSeveritiesToBeExcluded.end(),
1184                              thisEcalSevLvl);
1185       if (sevit!=theEcalSeveritiesToBeExcluded.end()) {
1186     ++numBadEcalChan;
1187       }
1188 
1189      }
1190      
1191      // compare with fast version
1192  
1193    //  hcal: 
1194     int inEcals[2] = {0,0};
1195     for (std::vector<std::pair<DetId,float> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i) {
1196       DetId detId = i->first;
1197       if(detId.det() == DetId::Ecal){
1198     if( detId.subdetId()==EcalBarrel ) inEcals[0] =1;
1199     else if( detId.subdetId()==EcalEndcap ) inEcals[1] =1;
1200       }
1201     }
1202 
1203      auto  numBadEcalChanNew = ecalBadChs[theTowerTopology->denseIndex(id)]+mt.numBadEcalCells; // - mt.numRecEcalCells
1204      if (int(numBadEcalChanNew)!=int(numBadEcalChan)) { 
1205        std::cout << "VI wrong " << ((inEcals[1]==1) ? "EE" : "" ) << id << " " << numBadEcalChanNew << " " << numBadEcalChan 
1206          << " " << mt.numBadEcalCells << " " << mt.numRecEcalCells << std::endl;
1207      }
1208      */
1209 
1210   numBadEcalChan = ecalBadChs[theTowerTopology->denseIndex(id)] + mt.numBadEcalCells;  // - mt.numRecEcalCells
1211 
1212   //--------------------------------------------------------------------------------------
1213 
1214   caloTower.setCaloTowerStatus(
1215       numBadHcalChan, numBadEcalChan, numRecHcalChan, numRecEcalChan, numProbHcalChan, numProbEcalChan);
1216 
1217   double maxCellE = -999.0;  // for storing the hottest cell E in the calotower
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       // need an extra check because of the funny towers that are empty except for the presence of an HO
1226       // hit in the constituents (JetMET wanted them saved)
1227       // This constituent is only used for storing the tower, but should not be concidered as a hot cell canditate for
1228       // configurations with useHO = false
1229 
1230       if (i->first.det() == DetId::Ecal) {  // ECAL
1231         maxCellE = i->second;
1232       } else {  // HCAL
1233         if (HcalDetId(i->first).subdet() != HcalOuter)
1234           maxCellE = i->second;
1235         else if (theHOIsUsed)
1236           maxCellE = i->second;
1237       }
1238 
1239     }  // found higher E cell
1240 
1241   }  // loop over matacontains
1242 
1243   caloTower.setConstituents(std::move(contains));
1244   caloTower.setHottestCellE(maxCellE);
1245 
1246   // std::cout << "CaloTowerVI " << nalgo << ' ' << caloTower.id() << ((inEcals[1]==1) ? "EE " : " " )  << caloTower.pt() << ' ' << caloTower.et() << ' ' << caloTower.mass() << ' '
1247   //           << caloTower.constituentsSize() <<' '<< caloTower.towerStatusWord() << std::endl;
1248 }
1249 
1250 void CaloTowersCreationAlgo::getThresholdAndWeight(const DetId& detId, double& threshold, double& weight) const {
1251   DetId::Detector det = detId.det();
1252   weight = 0;  // in case the hit is not identified
1253 
1254   if (det == DetId::Ecal) {
1255     // may or may not be EB.  We'll find out.
1256 
1257     EcalSubdetector subdet = (EcalSubdetector)(detId.subdetId());
1258     if (subdet == EcalBarrel) {
1259       if (ecalCuts == nullptr) {  // this means that ecalRecHitThresh_ is false
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) {  // this means cutsFromDB is false
1288         threshold = (depth == 1) ? theHBthreshold1 : (depth == 2) ? theHBthreshold2 : theHBthreshold;
1289       } else {  // hcalCuts is not nullptr, i.e. cutsFromDB is true
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       // check if it's single or double tower
1302       if (hcalDetId.ietaAbs() < theHcalTopology->firstHEDoublePhiRing()) {
1303         if (hcalCuts == nullptr) {  // this means cutsFromDB is false
1304           threshold = (depth == 1) ? theHESthreshold1 : theHESthreshold;
1305         } else {  // hcalCuts is not nullptr, i.e. cutsFromDB is true
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) {  // this means cutsFromDB is false
1316           threshold = (depth == 1) ? theHEDthreshold1 : theHEDthreshold;
1317         } else {  // hcalCuts is not nullptr, i.e. cutsFromDB is true
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       //check if it's ring 0 or +1 or +2 or -1 or -2
1331       if (hcalDetId.ietaAbs() <= 4)
1332         threshold = theHOthreshold0;
1333       else if (hcalDetId.ieta() < 0) {
1334         // set threshold for ring -1 or -2
1335         threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdMinus1 : theHOthresholdMinus2;
1336       } else {
1337         // set threshold for ring +1 or +2
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();  // face of the cell
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   // same code as above
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   // this is based on available RecHits, can lead to different actual depths if
1449   // hits in multi-depth towers are not all there
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     // do not use HO for deirection calculations for now
1468     if (HcalDetId(mc_it->first).subdet() == HcalOuter)
1469       continue;
1470     ++nConst;
1471 
1472     GlobalPoint p = hadSegmentShwrPos(mc_it->first, fracDepth);
1473 
1474     // longitudinal segmentation: do not weight by energy,
1475     // get the geometrical position
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   // set depth using geometry of cells that are associated with the
1486   // tower (regardless if they have non-zero energies)
1487 
1488   //  if (hadE <= 0) return GlobalPoint(0, 0, 0);
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     // forward, take the geometry for long fibers
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     //use constituents map
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     //fix for tower 28/29 - no tower 29 at highest depths
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   // uses the "front" and "back" cells
1570   // to determine the axis. point set by the predefined depth.
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;  // add crystals with E/E_EM > 1.5%
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;  // discretization (ns)
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 // Bad/anomolous cell handling
1676 
1677 void CaloTowersCreationAlgo::makeHcalDropChMap() {
1678   // This method fills the map of number of dead channels for the calotower,
1679   // The key of the map is CaloTowerDetId.
1680   // By definition these channels are not going to be in the RecHit collections.
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       // special case for tower 29: if HCAL hit is in depth 3 add to twr 29 as well
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   // now I know how many bad channels, but I also need to know if there's any good ones
1710   if (missingHcalRescaleFactorForEcal > 0) {
1711     for (auto& pair : hcalDropChMap) {
1712       if (pair.second.first == 0)
1713         continue;  // unexpected, but just in case
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;  // recount, since pair.second.first may include HO
1726         }
1727       }
1728       if (nbad > 0 && nbad >= ngood) {
1729         //uncomment for debug (may be useful to tune the criteria above)
1730         //CaloTowerDetId id(pair.first);
1731         //std::cout << "CaloTower at ieta = " << id.ieta() << ", iphi " << id.iphi() << ": set Hcal as not efficient (ngood =" << ngood << ", nbad = " << nbad << ")" << std::endl;
1732         pair.second.second = true;
1733       }
1734     }
1735   }
1736 }
1737 
1738 void CaloTowersCreationAlgo::makeEcalBadChs() {
1739   // std::cout << "VI making EcalBadChs ";
1740 
1741   // for ECAL the number of all bad channels is obtained here -----------------------
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     // this is utterly slow... (can be optmized if really needed)
1749 
1750     // get all possible constituents of the tower
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       // check if the Ecal severity is ok to keep
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     // if (0!=numBadEcalChan) std::cout << id << ":" << numBadEcalChan << ", ";
1768   }
1769 
1770   /*
1771   int tot=0;
1772   for (auto ind=0U; ind<theTowerTopology->sizeForDenseIndexing(); ++ind) {
1773     if (ecalBadChs[ind]!=0) ++tot; 
1774   }
1775   std::cout << " | " << tot << std::endl;
1776   */
1777 }
1778 
1779 //////  Get status of the channel contributing to the tower
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   // For use with hits rejected in the default reconstruction
1794   if (useRejectedHitsOnly) {
1795     if (!isRecovered) {
1796       if (severityLevel <= int(theHcalAcceptSeverityLevel) ||
1797           severityLevel > int(theHcalAcceptSeverityLevelForRejectedHit))
1798         return CaloTowersCreationAlgo::IgnoredChan;
1799       // this hit was either already accepted or is worse than
1800     } else {
1801       if (theRecoveredHcalHitsAreUsed || !useRejectedRecoveredHcalHits) {
1802         // skip recovered hits either because they were already used or because there was an explicit instruction
1803         return CaloTowersCreationAlgo::IgnoredChan;
1804       } else if (useRejectedRecoveredHcalHits) {
1805         return CaloTowersCreationAlgo::RecoveredChan;
1806       }
1807 
1808     }  // recovered channels
1809 
1810     // clasify channels as problematic: no good hits are supposed to be present in the
1811     // extra rechit collections
1812     return CaloTowersCreationAlgo::ProblematicChan;
1813 
1814   }  // treatment of rejected hits
1815 
1816   // this is for the regular reconstruction sequence
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   // const DetId id = hit->detid();
1834 
1835   //  uint16_t dbStatus = theEcalChStatus->find(id)->getStatusCode();
1836   //  uint32_t rhFlags  = hit->flags();
1837   //  int severityLevel = theEcalSevLvlAlgo->severityLevel(rhFlags, dbStatus);
1838   // The methods above will become private and cannot be usef for flagging ecal spikes.
1839   // Use the recommended interface - we leave the parameters for spilke removal to be specified by ECAL.
1840 
1841   //  int severityLevel = 999;
1842 
1843   EcalRecHit const& rh = *reinterpret_cast<EcalRecHit const*>(hit);
1844   int severityLevel = theEcalSevLvlAlgo->severityLevel(rh);
1845 
1846   //  if      (id.subdetId() == EcalBarrel) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEbHandle);//, *theEcalChStatus);
1847   //  else if (id.subdetId() == EcalEndcap) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEeHandle);//, *theEcalChStatus);
1848 
1849   // there should be no other ECAL types used in this reconstruction
1850 
1851   // The definition of ECAL severity levels uses categories that
1852   // are similar to the defined for CaloTower. (However, the categorization
1853   // for CaloTowers depends on the specified maximum acceptabel severity and therefore cannnot
1854   // be exact correspondence between the two. ECAL has additional categories describing modes of failure.)
1855   // This approach is different from the initial idea and from
1856   // the implementation for HCAL. Still make the logic similar to HCAL so that one has the ability to
1857   // exclude problematic channels as defined by ECAL.
1858   // For definitions of ECAL severity levels see RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h
1859 
1860   bool isBad = (severityLevel == EcalSeverityLevel::kBad);
1861 
1862   bool isRecovered = (severityLevel == EcalSeverityLevel::kRecovered);
1863 
1864   // check if the severity is compatible with our configuration
1865   // This applies to the "default" tower cleaning
1866   std::vector<int>::const_iterator sevit =
1867       std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), severityLevel);
1868   bool accepted = (sevit == theEcalSeveritiesToBeExcluded.end());
1869 
1870   // For use with hits that were rejected in the regular reconstruction:
1871   // This is for creating calotowers with lower level of cleaning by merging
1872   // the information from the default towers and a collection of towers created from
1873   // bad rechits
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       // this hit was either already accepted, or is not eligible for inclusion
1882     } else {
1883       if (theRecoveredEcalHitsAreUsed || !useRejectedRecoveredEcalHits) {
1884         // skip recovered hits either because they were already used or because there was an explicit instruction
1885         return std::make_tuple(CaloTowersCreationAlgo::IgnoredChan, isBad);
1886         ;
1887       } else if (useRejectedRecoveredEcalHits) {
1888         return std::make_tuple(CaloTowersCreationAlgo::RecoveredChan, isBad);
1889       }
1890 
1891     }  // recovered channels
1892 
1893     // clasify channels as problematic
1894     return std::make_tuple(CaloTowersCreationAlgo::ProblematicChan, isBad);
1895 
1896   }  // treatment of rejected hits
1897 
1898   // for normal reconstruction
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 }