CaloTowersCreator/src/CaloTowersCreationAlgo.cc

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::setGeometry(const CaloTowerTopology* cttopo,
0357                                          const CaloTowerConstituentsMap* ctmap,
0358                                          const HcalTopology* htopo,
0359                                          const CaloGeometry* geo) {
0360   theTowerTopology = cttopo;
0361   theTowerConstituentsMap = ctmap;
0362   theHcalTopology = htopo;
0363   theGeometry = geo;
0364   theTowerGeometry = geo->getSubdetectorGeometry(DetId::Calo, CaloTowerDetId::SubdetId);
0365
0366   //initialize ecal bad channel map
0367   ecalBadChs.resize(theTowerTopology->sizeForDenseIndexing(), 0);
0368 }
0369
0370 void CaloTowersCreationAlgo::begin() {
0371   theTowerMap.clear();
0372   theTowerMapSize = 0;
0373   //hcalDropChMap.clear();
0374 }
0375
0376 void CaloTowersCreationAlgo::process(const HBHERecHitCollection& hbhe) {
0377   for (HBHERecHitCollection::const_iterator hbheItr = hbhe.begin(); hbheItr != hbhe.end(); ++hbheItr)
0378     assignHitHcal(&(*hbheItr));
0379 }
0380
0381 void CaloTowersCreationAlgo::process(const HORecHitCollection& ho) {
0382   for (HORecHitCollection::const_iterator hoItr = ho.begin(); hoItr != ho.end(); ++hoItr)
0383     assignHitHcal(&(*hoItr));
0384 }
0385
0386 void CaloTowersCreationAlgo::process(const HFRecHitCollection& hf) {
0387   for (HFRecHitCollection::const_iterator hfItr = hf.begin(); hfItr != hf.end(); ++hfItr)
0388     assignHitHcal(&(*hfItr));
0389 }
0390
0391 void CaloTowersCreationAlgo::process(const EcalRecHitCollection& ec) {
0392   for (EcalRecHitCollection::const_iterator ecItr = ec.begin(); ecItr != ec.end(); ++ecItr)
0393     assignHitEcal(&(*ecItr));
0394 }
0395
0396 // this method should not be used any more as the towers in the changed format
0397 // can not be properly rescaled with the "rescale" method.
0398 // "rescale was replaced by "rescaleTowers"
0399 //
0400 void CaloTowersCreationAlgo::process(const CaloTowerCollection& ctc) {
0401   for (CaloTowerCollection::const_iterator ctcItr = ctc.begin(); ctcItr != ctc.end(); ++ctcItr) {
0402     rescale(&(*ctcItr));
0403   }
0404 }
0405
0406 void CaloTowersCreationAlgo::finish(CaloTowerCollection& result) {
0407   // now copy this map into the final collection
0408   result.reserve(theTowerMapSize);
0409   // auto k=0U;
0410   //  if (!theEbHandle.isValid()) std::cout << "VI ebHandle not valid" << std::endl;
0411   // if (!theEeHandle.isValid()) std::cout << "VI eeHandle not valid" << std::endl;
0412
0413   for (auto const& mt : theTowerMap) {
0414     // Convert only if there is at least one constituent in the metatower.
0415     // The check of constituents size in the coverted tower is still needed!
0416     if (!mt.empty()) {
0417       convert(mt.id, mt, result);
0418     }  // ++k;}
0419   }
0420
0421   // assert(k==theTowerMapSize);
0422   // std::cout << "VI TowerMap " << theTowerMapSize << " " << k << std::endl;
0423
0424   theTowerMap.clear();  // save the memory
0425   theTowerMapSize = 0;
0426 }
0427
0428 void CaloTowersCreationAlgo::rescaleTowers(const CaloTowerCollection& ctc, CaloTowerCollection& ctcResult) {
0429   for (CaloTowerCollection::const_iterator ctcItr = ctc.begin(); ctcItr != ctc.end(); ++ctcItr) {
0430     CaloTowerDetId twrId = ctcItr->id();
0431     double newE_em = ctcItr->emEnergy();
0432     double newE_had = ctcItr->hadEnergy();
0433     double newE_outer = ctcItr->outerEnergy();
0434
0435     double threshold = 0.0;  // not used: we do not change thresholds
0436     double weight = 1.0;
0437
0438     // HF
0439     if (ctcItr->ietaAbs() >= theTowerTopology->firstHFRing()) {
0440       double E_short = 0.5 * newE_had;           // from the definitions for HF
0441       double E_long = newE_em + 0.5 * newE_had;  //
0442       // scale
0443       E_long *= theHF1weight;
0444       E_short *= theHF2weight;
0445       // convert
0446       newE_em = E_long - E_short;
0447       newE_had = 2.0 * E_short;
0448     }
0449
0450     else {  // barrel/endcap
0451
0452       // find if its in EB, or EE; determine from first ecal constituent found
0453       for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0454         DetId constId = ctcItr->constituent(iConst);
0455         if (constId.det() != DetId::Ecal)
0456           continue;
0457         getThresholdAndWeight(constId, threshold, weight);
0458         newE_em *= weight;
0459         break;
0460       }
0461       // HO
0462       for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0463         DetId constId = ctcItr->constituent(iConst);
0464         if (constId.det() != DetId::Hcal)
0465           continue;
0466         if (HcalDetId(constId).subdet() != HcalOuter)
0467           continue;
0468         getThresholdAndWeight(constId, threshold, weight);
0469         newE_outer *= weight;
0470         break;
0471       }
0472       // HB/HE
0473       for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0474         DetId constId = ctcItr->constituent(iConst);
0475         if (constId.det() != DetId::Hcal)
0476           continue;
0477         if (HcalDetId(constId).subdet() == HcalOuter)
0478           continue;
0479         getThresholdAndWeight(constId, threshold, weight);
0480         newE_had *= weight;
0481         if (ctcItr->ietaAbs() > theTowerTopology->firstHERing())
0482           newE_outer *= weight;
0483         break;
0484       }
0485
0486     }  // barrel/endcap region
0487
0488     // now make the new tower
0489
0490     double newE_hadTot =
0491         (theHOIsUsed && twrId.ietaAbs() <= theTowerTopology->lastHORing()) ? newE_had + newE_outer : newE_had;
0492
0493     GlobalPoint emPoint = ctcItr->emPosition();
0494     GlobalPoint hadPoint = ctcItr->emPosition();
0495
0496     double f_em = 1.0 / cosh(emPoint.eta());
0497     double f_had = 1.0 / cosh(hadPoint.eta());
0498
0499     CaloTower::PolarLorentzVector towerP4;
0500
0501     if (ctcItr->ietaAbs() < theTowerTopology->firstHFRing()) {
0502       if (newE_em > 0)
0503         towerP4 += CaloTower::PolarLorentzVector(newE_em * f_em, emPoint.eta(), emPoint.phi(), 0);
0504       if (newE_hadTot > 0)
0505         towerP4 += CaloTower::PolarLorentzVector(newE_hadTot * f_had, hadPoint.eta(), hadPoint.phi(), 0);
0506     } else {
0507       double newE_tot = newE_em + newE_had;
0508       // for HF we use common point for ecal, hcal shower positions regardless of the method
0509       if (newE_tot > 0)
0510         towerP4 += CaloTower::PolarLorentzVector(newE_tot * f_had, hadPoint.eta(), hadPoint.phi(), 0);
0511     }
0512
0513     CaloTower rescaledTower(twrId, newE_em, newE_had, newE_outer, -1, -1, towerP4, emPoint, hadPoint);
0514     // copy the timings, have to convert back to int, 1 unit = 0.01 ns
0515     rescaledTower.setEcalTime(int(ctcItr->ecalTime() * 100.0 + 0.5));
0516     rescaledTower.setHcalTime(int(ctcItr->hcalTime() * 100.0 + 0.5));
0517     //add topology info
0518     rescaledTower.setHcalSubdet(theTowerTopology->lastHBRing(),
0519                                 theTowerTopology->lastHERing(),
0520                                 theTowerTopology->lastHFRing(),
0521                                 theTowerTopology->lastHORing());
0522
0523     std::vector<DetId> contains;
0524     for (unsigned int iConst = 0; iConst < ctcItr->constituentsSize(); ++iConst) {
0525       contains.push_back(ctcItr->constituent(iConst));
0526     }
0527     rescaledTower.addConstituents(contains);
0528
0529     rescaledTower.setCaloTowerStatus(ctcItr->towerStatusWord());
0530
0531     ctcResult.push_back(rescaledTower);
0532
0533   }  // end of loop over towers
0534 }
0535
0536 void CaloTowersCreationAlgo::assignHitHcal(const CaloRecHit* recHit) {
0537   DetId detId = recHit->detid();
0538   DetId detIdF(detId);
0539   if (detId.det() == DetId::Hcal && theHcalTopology->getMergePositionFlag()) {
0540     detIdF = theHcalTopology->idFront(HcalDetId(detId));
0541 #ifdef EDM_ML_DEBUG
0542     std::cout << "AssignHitHcal DetId " << HcalDetId(detId) << " Front " << HcalDetId(detIdF) << std::endl;
0543 #endif
0544   }
0545
0546   unsigned int chStatusForCT = hcalChanStatusForCaloTower(recHit);
0547
0548   // this is for skipping channls: mostly needed for the creation of
0549   // bad towers from hits i the bad channel collections.
0550   if (chStatusForCT == CaloTowersCreationAlgo::IgnoredChan)
0551     return;
0552
0553   double threshold, weight;
0554   getThresholdAndWeight(detId, threshold, weight);
0555
0556   double energy = recHit->energy();  // original RecHit energy is used to apply thresholds
0557   double e = energy * weight;        // energies scaled by user weight: used in energy assignments
0558
0559   // SPECIAL handling of tower 28 merged depths --> half into tower 28 and half into tower 29
0560   bool merge(false);
0561   if (detIdF.det() == DetId::Hcal && HcalDetId(detIdF).subdet() == HcalEndcap &&
0562       (theHcalPhase == 0 || theHcalPhase == 1) &&
0563       //HcalDetId(detId).depth()==3 &&
0564       HcalDetId(detIdF).ietaAbs() == theHcalTopology->lastHERing() - 1) {
0565     merge = theHcalTopology->mergedDepth29(HcalDetId(detIdF));
0566 #ifdef EDM_ML_DEBUG
0567     std::cout << "Merge " << HcalDetId(detIdF) << ":" << merge << std::endl;
0568 #endif
0569   }
0570   if (merge) {
0571     //////////////////////////////    unsigned int chStatusForCT = hcalChanStatusForCaloTower(recHit);
0572
0573     // bad channels are counted regardless of energy threshold
0574
0575     if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0576       CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0577       if (towerDetId.null())
0578         return;
0579       MetaTower& tower28 = find(towerDetId);
0580       CaloTowerDetId towerDetId29(towerDetId.ieta() + towerDetId.zside(), towerDetId.iphi());
0581       MetaTower& tower29 = find(towerDetId29);
0582       tower28.numBadHcalCells += 1;
0583       tower29.numBadHcalCells += 1;
0584     }
0585
0586     else if (0.5 * energy >= threshold) {  // not bad channel: use energy if above threshold
0587
0588       CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0589       if (towerDetId.null())
0590         return;
0591       MetaTower& tower28 = find(towerDetId);
0592       CaloTowerDetId towerDetId29(towerDetId.ieta() + towerDetId.zside(), towerDetId.iphi());
0593       MetaTower& tower29 = find(towerDetId29);
0594
0595       if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0596         tower28.numRecHcalCells += 1;
0597         tower29.numRecHcalCells += 1;
0598       } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0599         tower28.numProbHcalCells += 1;
0600         tower29.numProbHcalCells += 1;
0601       }
0602
0603       // NOTE DIVIDE BY 2!!!
0604       double e28 = 0.5 * e;
0605       double e29 = 0.5 * e;
0606
0607       tower28.E_had += e28;
0608       tower28.E += e28;
0609       std::pair<DetId, float> mc(detId, e28);
0610       tower28.metaConstituents.push_back(mc);
0611
0612       tower29.E_had += e29;
0613       tower29.E += e29;
0614       tower29.metaConstituents.push_back(mc);
0615
0616       // time info: do not use in averaging if timing error is found: need
0617       // full set of status info to implement: use only "good" channels for now
0618
0619       if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0620         tower28.hadSumTimeTimesE += (e28 * recHit->time());
0621         tower28.hadSumEForTime += e28;
0622         tower29.hadSumTimeTimesE += (e29 * recHit->time());
0623         tower29.hadSumEForTime += e29;
0624       }
0625
0626       // store the energy in layer 3 also in E_outer
0627       tower28.E_outer += e28;
0628       tower29.E_outer += e29;
0629     }  // not a "bad" hit
0630   }    // end of special case
0631
0632   else {
0633     HcalDetId hcalDetId(detId);
0634
0635     ///////////////////////      unsigned int chStatusForCT = hcalChanStatusForCaloTower(recHit);
0636
0637     if (hcalDetId.subdet() == HcalOuter) {
0638       CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0639       if (towerDetId.null())
0640         return;
0641       MetaTower& tower = find(towerDetId);
0642
0643       if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0644         if (theHOIsUsed)
0645           tower.numBadHcalCells += 1;
0646       }
0647
0648       else if (energy >= threshold) {
0649         tower.E_outer += e;  // store HO energy even if HO is not used
0650         // add energy of the tower and/or flag if theHOIsUsed
0651         if (theHOIsUsed) {
0652           tower.E += e;
0653
0654           if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0655             tower.numRecHcalCells += 1;
0656           } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0657             tower.numProbHcalCells += 1;
0658           }
0659         }  // HO is used
0660
0661         // add HO to constituents even if it is not used: JetMET wants to keep these towers
0662         std::pair<DetId, float> mc(detId, e);
0663         tower.metaConstituents.push_back(mc);
0664
0665       }  // not a bad channel, energy above threshold
0666
0667     }  // HO hit
0668
0669     // HF calculates EM fraction differently
0670     else if (hcalDetId.subdet() == HcalForward) {
0671       if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0672         CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0673         if (towerDetId.null())
0674           return;
0675         MetaTower& tower = find(towerDetId);
0676         tower.numBadHcalCells += 1;
0677       }
0678
0679       else if (energy >= threshold) {
0680         CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0681         if (towerDetId.null())
0682           return;
0683         MetaTower& tower = find(towerDetId);
0684
0685         if (hcalDetId.depth() == 1) {
0686           // long fiber, so E_EM = E(Long) - E(Short)
0687           tower.E_em += e;
0688         } else {
0689           // short fiber, EHAD = 2 * E(Short)
0690           tower.E_em -= e;
0691           tower.E_had += 2. * e;
0692         }
0693         tower.E += e;
0694         if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0695           tower.numRecHcalCells += 1;
0696         } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0697           tower.numProbHcalCells += 1;
0698         }
0699
0700         // put the timing in HCAL -> have to check timing errors when available
0701         // for now use only good channels
0702         if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0703           tower.hadSumTimeTimesE += (e * recHit->time());
0704           tower.hadSumEForTime += e;
0705         }
0706
0707         std::pair<DetId, float> mc(detId, e);
0708         tower.metaConstituents.push_back(mc);
0709
0710       }  // not a bad HF channel, energy above threshold
0711
0712     }  // HF hit
0713
0714     else {
0715       // HCAL situation normal in HB/HE
0716       if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0717         CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0718         if (towerDetId.null())
0719           return;
0720         MetaTower& tower = find(towerDetId);
0721         tower.numBadHcalCells += 1;
0722       } else if (energy >= threshold) {
0723         CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0724         if (towerDetId.null())
0725           return;
0726         MetaTower& tower = find(towerDetId);
0727         tower.E_had += e;
0728         tower.E += e;
0729         if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0730           tower.numRecHcalCells += 1;
0731         } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0732           tower.numProbHcalCells += 1;
0733         }
0734
0735         // Timing information: need specific accessors
0736         // for now use only good channels
0737         if (chStatusForCT == CaloTowersCreationAlgo::GoodChan) {
0738           tower.hadSumTimeTimesE += (e * recHit->time());
0739           tower.hadSumEForTime += e;
0740         }
0741         // store energy in highest depth for towers 18-27 (for electron,photon ID in endcap)
0742         // also, store energy in HE part of tower 16 (for JetMET cleanup)
0743         HcalDetId hcalDetId(detId);
0744         if (hcalDetId.subdet() == HcalEndcap) {
0745           if (theHcalPhase == 0) {
0746             if ((hcalDetId.depth() == 2 && hcalDetId.ietaAbs() >= 18 && hcalDetId.ietaAbs() < 27) ||
0747                 (hcalDetId.depth() == 3 && hcalDetId.ietaAbs() == 27) ||
0748                 (hcalDetId.depth() == 3 && hcalDetId.ietaAbs() == 16)) {
0749               tower.E_outer += e;
0750             }
0751           }
0752           //combine depths in phase0-like way
0753           else if (theHcalPhase == 1) {
0754             if ((hcalDetId.depth() >= 3 && hcalDetId.ietaAbs() >= 18 && hcalDetId.ietaAbs() < 26) ||
0755                 (hcalDetId.depth() >= 4 && (hcalDetId.ietaAbs() == 26 || hcalDetId.ietaAbs() == 27)) ||
0756                 (hcalDetId.depth() == 3 && hcalDetId.ietaAbs() == 17) ||
0757                 (hcalDetId.depth() == 4 && hcalDetId.ietaAbs() == 16)) {
0758               tower.E_outer += e;
0759             }
0760           }
0761         }
0762
0763         std::pair<DetId, float> mc(detId, e);
0764         tower.metaConstituents.push_back(mc);
0765
0766       }  // not a "bad" channel, energy above threshold
0767
0768     }  // channel in HBHE (excluding twrs 28,29)
0769
0770   }  // recHit normal case (not in HE towers 28,29)
0771
0772 }  // end of assignHitHcal method
0773
0774 void CaloTowersCreationAlgo::assignHitEcal(const EcalRecHit* recHit) {
0775   DetId detId = recHit->detid();
0776
0777   unsigned int chStatusForCT;
0778   bool ecalIsBad = false;
0779   std::tie(chStatusForCT, ecalIsBad) = ecalChanStatusForCaloTower(recHit);
0780
0781   // this is for skipping channls: mostly needed for the creation of
0782   // bad towers from hits i the bad channel collections.
0783   if (chStatusForCT == CaloTowersCreationAlgo::IgnoredChan)
0784     return;
0785
0786   double threshold, weight;
0787   getThresholdAndWeight(detId, threshold, weight);
0788
0789   double energy = recHit->energy();  // original RecHit energy is used to apply thresholds
0790   double e = energy * weight;        // energies scaled by user weight: used in energy assignments
0791
0792   /////////////////////////////      unsigned int chStatusForCT = ecalChanStatusForCaloTower(recHit);
0793
0794   // For ECAL we count all bad channels after the metatower is complete
0795
0796   // Include options for symmetric thresholds and cut on Et
0797   // for ECAL RecHits
0798
0799   bool passEmThreshold = false;
0800
0801   if (detId.subdetId() == EcalBarrel) {
0802     if (theUseEtEBTresholdFlag)
0803       energy /= cosh((theGeometry->getGeometry(detId)->getPosition()).eta());
0804     if (theUseSymEBTresholdFlag)
0805       passEmThreshold = (fabs(energy) >= threshold);
0806     else
0807       passEmThreshold = (energy >= threshold);
0808
0809   } else if (detId.subdetId() == EcalEndcap) {
0810     if (theUseEtEETresholdFlag)
0811       energy /= cosh((theGeometry->getGeometry(detId)->getPosition()).eta());
0812     if (theUseSymEETresholdFlag)
0813       passEmThreshold = (fabs(energy) >= threshold);
0814     else
0815       passEmThreshold = (energy >= threshold);
0816   }
0817
0818   CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(detId);
0819   if (towerDetId.null())
0820     return;
0821   MetaTower& tower = find(towerDetId);
0822
0823   // count bad cells and avoid double counting with those from DB (Recovered are counted bad)
0824
0825   // somehow misses some
0826   // if ( (chStatusForCT == CaloTowersCreationAlgo::BadChan) & (!ecalIsBad) ) ++tower.numBadEcalCells;
0827
0828   // a bit slower...
0829   if (chStatusForCT == CaloTowersCreationAlgo::BadChan) {
0830     auto thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel(detId);
0831     // check if the Ecal severity is ok to keep
0832     auto sevit = std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), thisEcalSevLvl);
0833     if (sevit == theEcalSeveritiesToBeExcluded.end())
0834       ++tower.numBadEcalCells;  // notinDB
0835   }
0836
0837   //      if (chStatusForCT != CaloTowersCreationAlgo::BadChan && energy >= threshold) {
0838   if (chStatusForCT != CaloTowersCreationAlgo::BadChan && passEmThreshold) {
0839     tower.E_em += e;
0840     tower.E += e;
0841
0842     if (chStatusForCT == CaloTowersCreationAlgo::RecoveredChan) {
0843       tower.numRecEcalCells += 1;
0844     } else if (chStatusForCT == CaloTowersCreationAlgo::ProblematicChan) {
0845       tower.numProbEcalCells += 1;
0846     }
0847
0848     // change when full status info is available
0849     // for now use only good channels
0850
0851     // add e>0 check (new options allow e<0)
0852     if (chStatusForCT == CaloTowersCreationAlgo::GoodChan && e > 0) {
0853       tower.emSumTimeTimesE += (e * recHit->time());
0854       tower.emSumEForTime += e;  // see above
0855     }
0856
0857     std::pair<DetId, float> mc(detId, e);
0858     tower.metaConstituents.push_back(mc);
0859   }
0860 }  // end of assignHitEcal method
0861
0862 // This method is not flexible enough for the new CaloTower format.
0863 // For now make a quick compatibility "fix" : WILL NOT WORK CORRECTLY with anything
0864 // except the default simple p4 assignment!!!
0865 // Must be rewritten for full functionality.
0866 void CaloTowersCreationAlgo::rescale(const CaloTower* ct) {
0867   double threshold, weight;
0868   CaloTowerDetId towerDetId = theTowerConstituentsMap->towerOf(ct->id());
0869   if (towerDetId.null())
0870     return;
0871   MetaTower& tower = find(towerDetId);
0872
0873   tower.E_em = 0.;
0874   tower.E_had = 0.;
0875   tower.E_outer = 0.;
0876   for (unsigned int i = 0; i < ct->constituentsSize(); i++) {
0877     DetId detId = ct->constituent(i);
0878     getThresholdAndWeight(detId, threshold, weight);
0879     DetId::Detector det = detId.det();
0880     if (det == DetId::Ecal) {
0881       tower.E_em = ct->emEnergy() * weight;
0882     } else {
0883       HcalDetId hcalDetId(detId);
0884       if (hcalDetId.subdet() == HcalForward) {
0885         if (hcalDetId.depth() == 1)
0886           tower.E_em = ct->emEnergy() * weight;
0887         if (hcalDetId.depth() == 2)
0888           tower.E_had = ct->hadEnergy() * weight;
0889       } else if (hcalDetId.subdet() == HcalOuter) {
0890         tower.E_outer = ct->outerEnergy() * weight;
0891       } else {
0892         tower.E_had = ct->hadEnergy() * weight;
0893       }
0894     }
0895     tower.E = tower.E_had + tower.E_em + tower.E_outer;
0896
0897     // this is to be compliant with the new MetaTower setup
0898     // used only for the default simple vector assignment
0899     std::pair<DetId, float> mc(detId, 0);
0900     tower.metaConstituents.push_back(mc);
0901   }
0902
0903   // preserve time inforamtion
0904   tower.emSumTimeTimesE = ct->ecalTime();
0905   tower.hadSumTimeTimesE = ct->hcalTime();
0906   tower.emSumEForTime = 1.0;
0907   tower.hadSumEForTime = 1.0;
0908 }
0909
0910 CaloTowersCreationAlgo::MetaTower& CaloTowersCreationAlgo::find(const CaloTowerDetId& detId) {
0911   if (theTowerMap.empty()) {
0912     theTowerMap.resize(theTowerTopology->sizeForDenseIndexing());
0913   }
0914
0915   auto& mt = theTowerMap[theTowerTopology->denseIndex(detId)];
0916
0917   if (mt.empty()) {
0918     mt.id = detId;
0919     mt.metaConstituents.reserve(detId.ietaAbs() < theTowerTopology->firstHFRing() ? 12 : 2);
0920     ++theTowerMapSize;
0921   }
0922
0923   return mt;
0924 }
0925
0926 void CaloTowersCreationAlgo::convert(const CaloTowerDetId& id, const MetaTower& mt, CaloTowerCollection& collection) {
0927   assert(id.rawId() != 0);
0928
0929   double ecalThres = (id.ietaAbs() <= 17) ? (theEBSumThreshold) : (theEESumThreshold);
0930   double E = mt.E;
0931   double E_em = mt.E_em;
0932   double E_had = mt.E_had;
0933   double E_outer = mt.E_outer;
0934
0935   // Note: E_outer is used to save HO energy OR energy in the outermost depths in endcap region
0936   // In the methods with separate treatment of EM and HAD components:
0937   //  - HO is not used to determine direction, however HO energy is added to get "total had energy"
0938   //  => Check if the tower is within HO coverage before adding E_outer to the "total had" energy
0939   //     else the energy will be double counted
0940   // When summing up the energy of the tower these checks are performed in the loops over RecHits
0941
0942   std::vector<std::pair<DetId, float> > metaContains = mt.metaConstituents;
0943   if (id.ietaAbs() < theTowerTopology->firstHFRing() && E_em < ecalThres) {  // ignore EM threshold in HF
0944     E -= E_em;
0945     E_em = 0;
0946     std::vector<std::pair<DetId, float> > metaContains_noecal;
0947
0948     for (std::vector<std::pair<DetId, float> >::iterator i = metaContains.begin(); i != metaContains.end(); ++i)
0949       if (i->first.det() != DetId::Ecal)
0950         metaContains_noecal.push_back(*i);
0951     metaContains.swap(metaContains_noecal);
0952   }
0953   if (id.ietaAbs() < theTowerTopology->firstHFRing() && E_had < theHcalThreshold) {
0954     E -= E_had;
0955
0956     if (theHOIsUsed && id.ietaAbs() <= theTowerTopology->lastHORing())
0957       E -= E_outer;  // not subtracted before, think it should be done
0958
0959     E_had = 0;
0960     E_outer = 0;
0961     std::vector<std::pair<DetId, float> > metaContains_nohcal;
0962
0963     for (std::vector<std::pair<DetId, float> >::iterator i = metaContains.begin(); i != metaContains.end(); ++i)
0964       if (i->first.det() != DetId::Hcal)
0965         metaContains_nohcal.push_back(*i);
0966     metaContains.swap(metaContains_nohcal);
0967   }
0968
0969   if (metaContains.empty())
0970     return;
0971
0972   if (missingHcalRescaleFactorForEcal > 0 && E_had == 0 && E_em > 0) {
0973     auto match = hcalDropChMap.find(id);
0974     if (match != hcalDropChMap.end() && match->second.second) {
0975       E_had = missingHcalRescaleFactorForEcal * E_em;
0976       E += E_had;
0977     }
0978   }
0979
0980   double E_had_tot = (theHOIsUsed && id.ietaAbs() <= theTowerTopology->lastHORing()) ? E_had + E_outer : E_had;
0981
0982   // create CaloTower using the selected algorithm
0983
0984   GlobalPoint emPoint, hadPoint;
0985
0986   // this is actually a 4D vector
0987   Basic3DVectorF towerP4;
0988   bool massless = true;
0989   // float mass1=0;
0990   float mass2 = 0;
0991
0992   // conditional assignment of depths for barrel/endcap
0993   // Some additional tuning may be required in the transitional region
0994   // 14<|iEta|<19
0995   double momEmDepth = 0.;
0996   double momHadDepth = 0.;
0997   if (id.ietaAbs() <= 17) {
0998     momHadDepth = theMomHBDepth;
0999     momEmDepth = theMomEBDepth;
1000   } else {
1001     momHadDepth = theMomHEDepth;
1002     momEmDepth = theMomEEDepth;
1003   }
1004
1005   switch (theMomConstrMethod) {
1006       // FIXME  : move to simple cartesian algebra
1007     case 0: {  // Simple 4-momentum assignment
1008       GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1009       towerP4 = p.basicVector().unit();
1010       towerP4[3] = 1.f;  // energy
1011       towerP4 *= E;
1012
1013       // double pf=1.0/cosh(p.eta());
1014       // if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);
1015
1016       emPoint = p;
1017       hadPoint = p;
1018     }  // end case 0
1019     break;
1020
1021     case 1: {  // separate 4-vectors for ECAL, HCAL, add to get the 4-vector of the tower (=>tower has mass!)
1022       if (id.ietaAbs() < theTowerTopology->firstHFRing()) {
1023         Basic3DVectorF emP4;
1024         if (E_em > 0) {
1025           emPoint = emShwrPos(metaContains, momEmDepth, E_em);
1026           emP4 = emPoint.basicVector().unit();
1027           emP4[3] = 1.f;  // energy
1028           towerP4 = emP4 * E_em;
1029
1030           // double emPf = 1.0/cosh(emPoint.eta());
1031           // towerP4 += CaloTower::PolarLorentzVector(E_em*emPf, emPoint.eta(), emPoint.phi(), 0);
1032         }
1033         if ((E_had + E_outer) > 0) {
1034           massless = (E_em <= 0);
1035           hadPoint = hadShwrPos(id, momHadDepth);
1036           auto lP4 = hadPoint.basicVector().unit();
1037           lP4[3] = 1.f;  // energy
1038           if (!massless) {
1039             auto diff = lP4 - emP4;
1040             mass2 = std::sqrt(E_em * E_had_tot * diff.mag2());
1041           }
1042           lP4 *= E_had_tot;
1043           towerP4 += lP4;
1044           /*
1045             if (!massless) {
1046                auto p = towerP4;
1047                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;
1048             }
1049             */
1050           // double hadPf = 1.0/cosh(hadPoint.eta());
1051           // if (E_had_tot>0) {
1052           //  towerP4 += CaloTower::PolarLorentzVector(E_had_tot*hadPf, hadPoint.eta(), hadPoint.phi(), 0);
1053           // }
1054         }
1055       } else {  // forward detector: use the CaloTower position
1056         GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1057         towerP4 = p.basicVector().unit();
1058         towerP4[3] = 1.f;  // energy
1059         towerP4 *= E;
1060         // double pf=1.0/cosh(p.eta());
1061         // if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);  // simple momentum assignment, same position
1062         emPoint = p;
1063         hadPoint = p;
1064       }
1065     }  // end case 1
1066     break;
1067
1068     case 2: {  // use ECAL position for the tower (when E_cal>0), else default CaloTower position (massless tower)
1069       if (id.ietaAbs() < theTowerTopology->firstHFRing()) {
1070         if (E_em > 0)
1071           emPoint = emShwrLogWeightPos(metaContains, momEmDepth, E_em);
1072         else
1073           emPoint = theTowerGeometry->getGeometry(id)->getPosition();
1074         towerP4 = emPoint.basicVector().unit();
1075         towerP4[3] = 1.f;  // energy
1076         towerP4 *= E;
1077
1078         // double sumPf = 1.0/cosh(emPoint.eta());
1079         /// if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*sumPf, emPoint.eta(), emPoint.phi(), 0);
1080
1081         hadPoint = emPoint;
1082       } else {  // forward detector: use the CaloTower position
1083         GlobalPoint p = theTowerGeometry->getGeometry(id)->getPosition();
1084         towerP4 = p.basicVector().unit();
1085         towerP4[3] = 1.f;  // energy
1086         towerP4 *= E;
1087
1088         // double pf=1.0/cosh(p.eta());
1089         // if (E>0) towerP4 = CaloTower::PolarLorentzVector(E*pf, p.eta(), p.phi(), 0);  // simple momentum assignment, same position
1090         emPoint = p;
1091         hadPoint = p;
1092       }
1093     }  // end case 2
1094     break;
1095
1096   }  // end of decision on p4 reconstruction method
1097
1098   // insert in collection (remove and return if below threshold)
1099   if UNLIKELY ((towerP4[3] == 0) & (E_outer > 0)) {
1100     float val = theHOIsUsed ? 0 : 1E-9;  // to keep backwards compatibility for theHOIsUsed == true
1101     collection.emplace_back(id,
1102                             E_em,
1103                             E_had,
1104                             E_outer,
1105                             -1,
1106                             -1,
1107                             CaloTower::PolarLorentzVector(val, hadPoint.eta(), hadPoint.phi(), 0),
1108                             emPoint,
1109                             hadPoint);
1110   } else {
1111     collection.emplace_back(
1112         id, E_em, E_had, E_outer, -1, -1, GlobalVector(towerP4), towerP4[3], mass2, emPoint, hadPoint);
1113   }
1114   auto& caloTower = collection.back();
1115
1116   // if (!massless) std::cout << "massive " << id <<' ' << mass1 <<' ' << mass2 <<' ' <<  caloTower.mass() << std::endl;
1117   // std::cout << "CaloTowerVI " <<theMomConstrMethod <<' ' << id <<' '<< E_em <<' '<< E_had <<' '<< E_outer <<' '<< GlobalVector(towerP4) <<' '<< towerP4[3] <<' '<< emPoint <<' '<< hadPoint << std::endl;
1118   //if (towerP4[3]==0) std::cout << "CaloTowerVIzero " << theEcutTower << ' ' << collection.back().eta() <<' '<< collection.back().phi() << std::endl;
1119
1120   if (caloTower.energy() < theEcutTower) {
1121     collection.pop_back();
1122     return;
1123   }
1124
1125   // set the timings
1126   float ecalTime = (mt.emSumEForTime > 0) ? mt.emSumTimeTimesE / mt.emSumEForTime : -9999;
1127   float hcalTime = (mt.hadSumEForTime > 0) ? mt.hadSumTimeTimesE / mt.hadSumEForTime : -9999;
1128   caloTower.setEcalTime(compactTime(ecalTime));
1129   caloTower.setHcalTime(compactTime(hcalTime));
1130   //add topology info
1131   caloTower.setHcalSubdet(theTowerTopology->lastHBRing(),
1132                           theTowerTopology->lastHERing(),
1133                           theTowerTopology->lastHFRing(),
1134                           theTowerTopology->lastHORing());
1135
1136   // set the CaloTower status word =====================================
1137   // Channels must be counter exclusively in the defined cathegories
1138   // "Bad" channels (not used in energy assignment) can be flagged during
1139   // CaloTower creation only if specified in the configuration file
1140
1141   unsigned int numBadHcalChan = mt.numBadHcalCells;
1142   //    unsigned int numBadEcalChan  = mt.numBadEcalCells;
1143   unsigned int numBadEcalChan = 0;  //
1144
1145   unsigned int numRecHcalChan = mt.numRecHcalCells;
1146   unsigned int numRecEcalChan = mt.numRecEcalCells;
1147   unsigned int numProbHcalChan = mt.numProbHcalCells;
1148   unsigned int numProbEcalChan = mt.numProbEcalCells;
1149
1150   // now add dead/off/... channels not used in RecHit reconstruction for HCAL
1151   HcalDropChMap::iterator dropChItr = hcalDropChMap.find(id);
1152   if (dropChItr != hcalDropChMap.end())
1153     numBadHcalChan += dropChItr->second.first;
1154
1155   // for ECAL the number of all bad channels is obtained here -----------------------
1156
1157   /*
1158     // old hyper slow algorithm
1159     // get all possible constituents of the tower
1160     std::vector<DetId> allConstituents = theTowerConstituentsMap->constituentsOf(id);
1161
1162     for (std::vector<DetId>::iterator ac_it=allConstituents.begin();
1163      ac_it!=allConstituents.end(); ++ac_it) {
1164
1165       if (ac_it->det()!=DetId::Ecal) continue;
1166
1167       int thisEcalSevLvl = -999;
1168
1169       if (ac_it->subdetId() == EcalBarrel && theEbHandle.isValid()) {
1170     thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEbHandle);//, *theEcalChStatus);
1171       }
1172       else if (ac_it->subdetId() == EcalEndcap && theEeHandle.isValid()) {
1173     thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel( *ac_it, *theEeHandle);//, *theEcalChStatus);
1174       }
1175
1176       // check if the Ecal severity is ok to keep
1177       std::vector<int>::const_iterator sevit = std::find(theEcalSeveritiesToBeExcluded.begin(),
1178                              theEcalSeveritiesToBeExcluded.end(),
1179                              thisEcalSevLvl);
1180       if (sevit!=theEcalSeveritiesToBeExcluded.end()) {
1181     ++numBadEcalChan;
1182       }
1183
1184      }
1185
1186      // compare with fast version
1187
1188    //  hcal:
1189     int inEcals[2] = {0,0};
1190     for (std::vector<std::pair<DetId,float> >::iterator i=metaContains.begin(); i!=metaContains.end(); ++i) {
1191       DetId detId = i->first;
1192       if(detId.det() == DetId::Ecal){
1193     if( detId.subdetId()==EcalBarrel ) inEcals[0] =1;
1194     else if( detId.subdetId()==EcalEndcap ) inEcals[1] =1;
1195       }
1196     }
1197
1198      auto  numBadEcalChanNew = ecalBadChs[theTowerTopology->denseIndex(id)]+mt.numBadEcalCells; // - mt.numRecEcalCells
1199      if (int(numBadEcalChanNew)!=int(numBadEcalChan)) {
1200        std::cout << "VI wrong " << ((inEcals[1]==1) ? "EE" : "" ) << id << " " << numBadEcalChanNew << " " << numBadEcalChan
1201          << " " << mt.numBadEcalCells << " " << mt.numRecEcalCells << std::endl;
1202      }
1203      */
1204
1205   numBadEcalChan = ecalBadChs[theTowerTopology->denseIndex(id)] + mt.numBadEcalCells;  // - mt.numRecEcalCells
1206
1207   //--------------------------------------------------------------------------------------
1208
1209   caloTower.setCaloTowerStatus(
1210       numBadHcalChan, numBadEcalChan, numRecHcalChan, numRecEcalChan, numProbHcalChan, numProbEcalChan);
1211
1212   double maxCellE = -999.0;  // for storing the hottest cell E in the calotower
1213
1214   std::vector<DetId> contains;
1215   contains.reserve(metaContains.size());
1216   for (std::vector<std::pair<DetId, float> >::iterator i = metaContains.begin(); i != metaContains.end(); ++i) {
1217     contains.push_back(i->first);
1218
1219     if (maxCellE < i->second) {
1220       // need an extra check because of the funny towers that are empty except for the presence of an HO
1221       // hit in the constituents (JetMET wanted them saved)
1222       // This constituent is only used for storing the tower, but should not be concidered as a hot cell canditate for
1223       // configurations with useHO = false
1224
1225       if (i->first.det() == DetId::Ecal) {  // ECAL
1226         maxCellE = i->second;
1227       } else {  // HCAL
1228         if (HcalDetId(i->first).subdet() != HcalOuter)
1229           maxCellE = i->second;
1230         else if (theHOIsUsed)
1231           maxCellE = i->second;
1232       }
1233
1234     }  // found higher E cell
1235
1236   }  // loop over matacontains
1237
1238   caloTower.setConstituents(std::move(contains));
1239   caloTower.setHottestCellE(maxCellE);
1240
1241   // std::cout << "CaloTowerVI " << nalgo << ' ' << caloTower.id() << ((inEcals[1]==1) ? "EE " : " " )  << caloTower.pt() << ' ' << caloTower.et() << ' ' << caloTower.mass() << ' '
1242   //           << caloTower.constituentsSize() <<' '<< caloTower.towerStatusWord() << std::endl;
1243 }
1244
1245 void CaloTowersCreationAlgo::getThresholdAndWeight(const DetId& detId, double& threshold, double& weight) const {
1246   DetId::Detector det = detId.det();
1247   weight = 0;  // in case the hit is not identified
1248
1249   if (det == DetId::Ecal) {
1250     // may or may not be EB.  We'll find out.
1251
1252     EcalSubdetector subdet = (EcalSubdetector)(detId.subdetId());
1253     if (subdet == EcalBarrel) {
1254       threshold = theEBthreshold;
1255       weight = theEBweight;
1256       if (weight <= 0.) {
1257         ROOT::Math::Interpolator my(theEBGrid, theEBWeights, ROOT::Math::Interpolation::kAKIMA);
1258         weight = my.Eval(theEBEScale);
1259       }
1260     } else if (subdet == EcalEndcap) {
1261       threshold = theEEthreshold;
1262       weight = theEEweight;
1263       if (weight <= 0.) {
1264         ROOT::Math::Interpolator my(theEEGrid, theEEWeights, ROOT::Math::Interpolation::kAKIMA);
1265         weight = my.Eval(theEEEScale);
1266       }
1267     }
1268   } else if (det == DetId::Hcal) {
1269     HcalDetId hcalDetId(detId);
1270     HcalSubdetector subdet = hcalDetId.subdet();
1271     int depth = hcalDetId.depth();
1272
1273     if (subdet == HcalBarrel) {
1274       threshold = (depth == 1) ? theHBthreshold1 : (depth == 2) ? theHBthreshold2 : theHBthreshold;
1275       weight = theHBweight;
1276       if (weight <= 0.) {
1277         ROOT::Math::Interpolator my(theHBGrid, theHBWeights, ROOT::Math::Interpolation::kAKIMA);
1278         weight = my.Eval(theHBEScale);
1279       }
1280     }
1281
1282     else if (subdet == HcalEndcap) {
1283       // check if it's single or double tower
1284       if (hcalDetId.ietaAbs() < theHcalTopology->firstHEDoublePhiRing()) {
1285         threshold = (depth == 1) ? theHESthreshold1 : theHESthreshold;
1286         weight = theHESweight;
1287         if (weight <= 0.) {
1288           ROOT::Math::Interpolator my(theHESGrid, theHESWeights, ROOT::Math::Interpolation::kAKIMA);
1289           weight = my.Eval(theHESEScale);
1290         }
1291       } else {
1292         threshold = (depth == 1) ? theHEDthreshold1 : theHEDthreshold;
1293         weight = theHEDweight;
1294         if (weight <= 0.) {
1295           ROOT::Math::Interpolator my(theHEDGrid, theHEDWeights, ROOT::Math::Interpolation::kAKIMA);
1296           weight = my.Eval(theHEDEScale);
1297         }
1298       }
1299     }
1300
1301     else if (subdet == HcalOuter) {
1302       //check if it's ring 0 or +1 or +2 or -1 or -2
1303       if (hcalDetId.ietaAbs() <= 4)
1304         threshold = theHOthreshold0;
1305       else if (hcalDetId.ieta() < 0) {
1306         // set threshold for ring -1 or -2
1307         threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdMinus1 : theHOthresholdMinus2;
1308       } else {
1309         // set threshold for ring +1 or +2
1310         threshold = (hcalDetId.ietaAbs() <= 10) ? theHOthresholdPlus1 : theHOthresholdPlus2;
1311       }
1312       weight = theHOweight;
1313       if (weight <= 0.) {
1314         ROOT::Math::Interpolator my(theHOGrid, theHOWeights, ROOT::Math::Interpolation::kAKIMA);
1315         weight = my.Eval(theHOEScale);
1316       }
1317     }
1318
1319     else if (subdet == HcalForward) {
1320       if (hcalDetId.depth() == 1) {
1321         threshold = theHF1threshold;
1322         weight = theHF1weight;
1323         if (weight <= 0.) {
1324           ROOT::Math::Interpolator my(theHF1Grid, theHF1Weights, ROOT::Math::Interpolation::kAKIMA);
1325           weight = my.Eval(theHF1EScale);
1326         }
1327       } else {
1328         threshold = theHF2threshold;
1329         weight = theHF2weight;
1330         if (weight <= 0.) {
1331           ROOT::Math::Interpolator my(theHF2Grid, theHF2Weights, ROOT::Math::Interpolation::kAKIMA);
1332           weight = my.Eval(theHF2EScale);
1333         }
1334       }
1335     }
1336   } else {
1337     edm::LogError("CaloTowersCreationAlgo") << "Bad cell: " << det << std::endl;
1338   }
1339 }
1340
1341 void CaloTowersCreationAlgo::setEBEScale(double scale) {
1342   if (scale > 0.00001)
1343     *&theEBEScale = scale;
1344   else
1345     *&theEBEScale = 50.;
1346 }
1347
1348 void CaloTowersCreationAlgo::setEEEScale(double scale) {
1349   if (scale > 0.00001)
1350     *&theEEEScale = scale;
1351   else
1352     *&theEEEScale = 50.;
1353 }
1354
1355 void CaloTowersCreationAlgo::setHBEScale(double scale) {
1356   if (scale > 0.00001)
1357     *&theHBEScale = scale;
1358   else
1359     *&theHBEScale = 50.;
1360 }
1361
1362 void CaloTowersCreationAlgo::setHESEScale(double scale) {
1363   if (scale > 0.00001)
1364     *&theHESEScale = scale;
1365   else
1366     *&theHESEScale = 50.;
1367 }
1368
1369 void CaloTowersCreationAlgo::setHEDEScale(double scale) {
1370   if (scale > 0.00001)
1371     *&theHEDEScale = scale;
1372   else
1373     *&theHEDEScale = 50.;
1374 }
1375
1376 void CaloTowersCreationAlgo::setHOEScale(double scale) {
1377   if (scale > 0.00001)
1378     *&theHOEScale = scale;
1379   else
1380     *&theHOEScale = 50.;
1381 }
1382
1383 void CaloTowersCreationAlgo::setHF1EScale(double scale) {
1384   if (scale > 0.00001)
1385     *&theHF1EScale = scale;
1386   else
1387     *&theHF1EScale = 50.;
1388 }
1389
1390 void CaloTowersCreationAlgo::setHF2EScale(double scale) {
1391   if (scale > 0.00001)
1392     *&theHF2EScale = scale;
1393   else
1394     *&theHF2EScale = 50.;
1395 }
1396
1397 GlobalPoint CaloTowersCreationAlgo::emCrystalShwrPos(DetId detId, float fracDepth) {
1398   auto cellGeometry = theGeometry->getGeometry(detId);
1399   GlobalPoint point = cellGeometry->getPosition();  // face of the cell
1400
1401   if (fracDepth <= 0)
1402     return point;
1403   if (fracDepth > 1)
1404     fracDepth = 1;
1405
1406   const GlobalPoint& backPoint = cellGeometry->getBackPoint();
1407   point += fracDepth * (backPoint - point);
1408
1409   return point;
1410 }
1411
1412 GlobalPoint CaloTowersCreationAlgo::hadSegmentShwrPos(DetId detId, float fracDepth) {
1413   // same code as above
1414   return emCrystalShwrPos(detId, fracDepth);
1415 }
1416
1417 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(const std::vector<std::pair<DetId, float> >& metaContains,
1418                                                float fracDepth,
1419                                                double hadE) {
1420   // this is based on available RecHits, can lead to different actual depths if
1421   // hits in multi-depth towers are not all there
1422 #ifdef EDM_ML_DEBUG
1423   std::cout << "hadShwrPos called with " << metaContains.size() << " elements and energy " << hadE << ":" << fracDepth
1424             << std::endl;
1425 #endif
1426   if (hadE <= 0)
1427     return GlobalPoint(0, 0, 0);
1428
1429   double hadX = 0.0;
1430   double hadY = 0.0;
1431   double hadZ = 0.0;
1432
1433   int nConst = 0;
1434
1435   std::vector<std::pair<DetId, float> >::const_iterator mc_it = metaContains.begin();
1436   for (; mc_it != metaContains.end(); ++mc_it) {
1437     if (mc_it->first.det() != DetId::Hcal)
1438       continue;
1439     // do not use HO for deirection calculations for now
1440     if (HcalDetId(mc_it->first).subdet() == HcalOuter)
1441       continue;
1442     ++nConst;
1443
1444     GlobalPoint p = hadSegmentShwrPos(mc_it->first, fracDepth);
1445
1446     // longitudinal segmentation: do not weight by energy,
1447     // get the geometrical position
1448     hadX += p.x();
1449     hadY += p.y();
1450     hadZ += p.z();
1451   }
1452
1453   return GlobalPoint(hadX / nConst, hadY / nConst, hadZ / nConst);
1454 }
1455
1456 GlobalPoint CaloTowersCreationAlgo::hadShwrPos(CaloTowerDetId towerId, float fracDepth) {
1457   // set depth using geometry of cells that are associated with the
1458   // tower (regardless if they have non-zero energies)
1459
1460   //  if (hadE <= 0) return GlobalPoint(0, 0, 0);
1461 #ifdef EDM_ML_DEBUG
1462   std::cout << "hadShwrPos " << towerId << " frac " << fracDepth << std::endl;
1463 #endif
1464   if (fracDepth < 0)
1465     fracDepth = 0;
1466   else if (fracDepth > 1)
1467     fracDepth = 1;
1468
1469   GlobalPoint point(0, 0, 0);
1470
1471   int iEta = towerId.ieta();
1472   int iPhi = towerId.iphi();
1473
1474   HcalDetId frontCellId, backCellId;
1475
1476   if (towerId.ietaAbs() >= theTowerTopology->firstHFRing()) {
1477     // forward, take the geometry for long fibers
1478     frontCellId = HcalDetId(HcalForward, towerId.zside() * theTowerTopology->convertCTtoHcal(abs(iEta)), iPhi, 1);
1479     backCellId = HcalDetId(HcalForward, towerId.zside() * theTowerTopology->convertCTtoHcal(abs(iEta)), iPhi, 1);
1480   } else {
1481     //use constituents map
1482     std::vector<DetId> items = theTowerConstituentsMap->constituentsOf(towerId);
1483     int frontDepth = 1000;
1484     int backDepth = -1000;
1485     for (unsigned i = 0; i < items.size(); i++) {
1486       if (items[i].det() != DetId::Hcal)
1487         continue;
1488       HcalDetId hid(items[i]);
1489       if (hid.subdet() == HcalOuter)
1490         continue;
1491       if (!theHcalTopology->validHcal(hid, 2))
1492         continue;
1493
1494       if (theHcalTopology->idFront(hid).depth() < frontDepth) {
1495         frontCellId = hid;
1496         frontDepth = theHcalTopology->idFront(hid).depth();
1497       }
1498       if (theHcalTopology->idBack(hid).depth() > backDepth) {
1499         backCellId = hid;
1500         backDepth = theHcalTopology->idBack(hid).depth();
1501       }
1502     }
1503 #ifdef EDM_ML_DEBUG
1504     std::cout << "Front " << frontCellId << " Back " << backCellId << " Depths " << frontDepth << ":" << backDepth
1505               << std::endl;
1506 #endif
1507     //fix for tower 28/29 - no tower 29 at highest depths
1508     if (towerId.ietaAbs() == theTowerTopology->lastHERing() && (theHcalPhase == 0 || theHcalPhase == 1)) {
1509       CaloTowerDetId towerId28(towerId.ieta() - towerId.zside(), towerId.iphi());
1510       std::vector<DetId> items28 = theTowerConstituentsMap->constituentsOf(towerId28);
1511 #ifdef EDM_ML_DEBUG
1512       std::cout << towerId28 << " with " << items28.size() << " constituents:";
1513       for (unsigned k = 0; k < items28.size(); ++k)
1514         if (items28[k].det() == DetId::Hcal)
1515           std::cout << " " << HcalDetId(items28[k]);
1516       std::cout << std::endl;
1517 #endif
1518       for (unsigned i = 0; i < items28.size(); i++) {
1519         if (items28[i].det() != DetId::Hcal)
1520           continue;
1521         HcalDetId hid(items28[i]);
1522         if (hid.subdet() == HcalOuter)
1523           continue;
1524
1525         if (theHcalTopology->idBack(hid).depth() > backDepth) {
1526           backCellId = hid;
1527           backDepth = theHcalTopology->idBack(hid).depth();
1528         }
1529       }
1530     }
1531 #ifdef EDM_ML_DEBUG
1532     std::cout << "Back " << backDepth << " ID " << backCellId << std::endl;
1533 #endif
1534   }
1535   point = hadShwPosFromCells(DetId(frontCellId), DetId(backCellId), fracDepth);
1536
1537   return point;
1538 }
1539
1540 GlobalPoint CaloTowersCreationAlgo::hadShwPosFromCells(DetId frontCellId, DetId backCellId, float fracDepth) {
1541   // uses the "front" and "back" cells
1542   // to determine the axis. point set by the predefined depth.
1543
1544   HcalDetId hid1(frontCellId), hid2(backCellId);
1545   if (theHcalTopology->getMergePositionFlag()) {
1546     hid1 = theHcalTopology->idFront(frontCellId);
1547 #ifdef EDM_ML_DEBUG
1548     std::cout << "Front " << HcalDetId(frontCellId) << " " << hid1 << "\n";
1549 #endif
1550     hid2 = theHcalTopology->idBack(backCellId);
1551 #ifdef EDM_ML_DEBUG
1552     std::cout << "Back  " << HcalDetId(backCellId) << " " << hid2 << "\n";
1553 #endif
1554   }
1555
1556   auto frontCellGeometry = theGeometry->getGeometry(DetId(hid1));
1557   auto backCellGeometry = theGeometry->getGeometry(DetId(hid2));
1558
1559   GlobalPoint point = frontCellGeometry->getPosition();
1560   const GlobalPoint& backPoint = backCellGeometry->getBackPoint();
1561
1562   point += fracDepth * (backPoint - point);
1563
1564   return point;
1565 }
1566
1567 GlobalPoint CaloTowersCreationAlgo::emShwrPos(const std::vector<std::pair<DetId, float> >& metaContains,
1568                                               float fracDepth,
1569                                               double emE) {
1570   if (emE <= 0)
1571     return GlobalPoint(0, 0, 0);
1572
1573   double emX = 0.0;
1574   double emY = 0.0;
1575   double emZ = 0.0;
1576
1577   double eSum = 0;
1578
1579   std::vector<std::pair<DetId, float> >::const_iterator mc_it = metaContains.begin();
1580   for (; mc_it != metaContains.end(); ++mc_it) {
1581     if (mc_it->first.det() != DetId::Ecal)
1582       continue;
1583     GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
1584     double e = mc_it->second;
1585
1586     if (e > 0) {
1587       emX += p.x() * e;
1588       emY += p.y() * e;
1589       emZ += p.z() * e;
1590       eSum += e;
1591     }
1592   }
1593
1594   return GlobalPoint(emX / eSum, emY / eSum, emZ / eSum);
1595 }
1596
1597 GlobalPoint CaloTowersCreationAlgo::emShwrLogWeightPos(const std::vector<std::pair<DetId, float> >& metaContains,
1598                                                        float fracDepth,
1599                                                        double emE) {
1600   double emX = 0.0;
1601   double emY = 0.0;
1602   double emZ = 0.0;
1603
1604   double weight = 0;
1605   double sumWeights = 0;
1606   double sumEmE = 0;  // add crystals with E/E_EM > 1.5%
1607   double crystalThresh = 0.015 * emE;
1608
1609   std::vector<std::pair<DetId, float> >::const_iterator mc_it = metaContains.begin();
1610   for (; mc_it != metaContains.end(); ++mc_it) {
1611     if (mc_it->second < 0)
1612       continue;
1613     if (mc_it->first.det() == DetId::Ecal && mc_it->second > crystalThresh)
1614       sumEmE += mc_it->second;
1615   }
1616
1617   for (mc_it = metaContains.begin(); mc_it != metaContains.end(); ++mc_it) {
1618     if (mc_it->first.det() != DetId::Ecal || mc_it->second < crystalThresh)
1619       continue;
1620
1621     GlobalPoint p = emCrystalShwrPos(mc_it->first, fracDepth);
1622
1623     weight = 4.2 + log(mc_it->second / sumEmE);
1624     sumWeights += weight;
1625
1626     emX += p.x() * weight;
1627     emY += p.y() * weight;
1628     emZ += p.z() * weight;
1629   }
1630
1631   return GlobalPoint(emX / sumWeights, emY / sumWeights, emZ / sumWeights);
1632 }
1633
1634 int CaloTowersCreationAlgo::compactTime(float time) {
1635   const float timeUnit = 0.01;  // discretization (ns)
1636
1637   if (time > 300.0)
1638     return 30000;
1639   if (time < -300.0)
1640     return -30000;
1641
1642   return int(time / timeUnit + 0.5);
1643 }
1644
1645 //========================================================
1646 //
1647 // Bad/anomolous cell handling
1648
1649 void CaloTowersCreationAlgo::makeHcalDropChMap() {
1650   // This method fills the map of number of dead channels for the calotower,
1651   // The key of the map is CaloTowerDetId.
1652   // By definition these channels are not going to be in the RecHit collections.
1653   hcalDropChMap.clear();
1654   std::vector<DetId> allChanInStatusCont = theHcalChStatus->getAllChannels();
1655
1656 #ifdef EDM_ML_DEBUG
1657   std::cout << "DropChMap with " << allChanInStatusCont.size() << " channels" << std::endl;
1658 #endif
1659   for (std::vector<DetId>::iterator it = allChanInStatusCont.begin(); it != allChanInStatusCont.end(); ++it) {
1660     const uint32_t dbStatusFlag = theHcalChStatus->getValues(*it)->getValue();
1661     if (theHcalSevLvlComputer->dropChannel(dbStatusFlag)) {
1662       DetId id = theHcalTopology->mergedDepthDetId(HcalDetId(*it));
1663
1664       CaloTowerDetId twrId = theTowerConstituentsMap->towerOf(id);
1665
1666       hcalDropChMap[twrId].first += 1;
1667
1668       HcalDetId hid(*it);
1669
1670       // special case for tower 29: if HCAL hit is in depth 3 add to twr 29 as well
1671       if (hid.subdet() == HcalEndcap && (theHcalPhase == 0 || theHcalPhase == 1) &&
1672           hid.ietaAbs() == theHcalTopology->lastHERing() - 1) {
1673         bool merge = theHcalTopology->mergedDepth29(hid);
1674         if (merge) {
1675           CaloTowerDetId twrId29(twrId.ieta() + twrId.zside(), twrId.iphi());
1676           hcalDropChMap[twrId29].first += 1;
1677         }
1678       }
1679     }
1680   }
1681   // now I know how many bad channels, but I also need to know if there's any good ones
1682   if (missingHcalRescaleFactorForEcal > 0) {
1683     for (auto& pair : hcalDropChMap) {
1684       if (pair.second.first == 0)
1685         continue;  // unexpected, but just in case
1686       int ngood = 0, nbad = 0;
1687       for (DetId id : theTowerConstituentsMap->constituentsOf(pair.first)) {
1688         if (id.det() != DetId::Hcal)
1689           continue;
1690         HcalDetId hid(id);
1691         if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap)
1692           continue;
1693         const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1694         if (dbStatusFlag == 0 || !theHcalSevLvlComputer->dropChannel(dbStatusFlag)) {
1695           ngood += 1;
1696         } else {
1697           nbad += 1;  // recount, since pair.second.first may include HO
1698         }
1699       }
1700       if (nbad > 0 && nbad >= ngood) {
1701         //uncomment for debug (may be useful to tune the criteria above)
1702         //CaloTowerDetId id(pair.first);
1703         //std::cout << "CaloTower at ieta = " << id.ieta() << ", iphi " << id.iphi() << ": set Hcal as not efficient (ngood =" << ngood << ", nbad = " << nbad << ")" << std::endl;
1704         pair.second.second = true;
1705       }
1706     }
1707   }
1708 }
1709
1710 void CaloTowersCreationAlgo::makeEcalBadChs() {
1711   // std::cout << "VI making EcalBadChs ";
1712
1713   // for ECAL the number of all bad channels is obtained here -----------------------
1714
1715   for (auto ind = 0U; ind < theTowerTopology->sizeForDenseIndexing(); ++ind) {
1716     auto& numBadEcalChan = ecalBadChs[ind];
1717     numBadEcalChan = 0;
1718     auto id = theTowerTopology->detIdFromDenseIndex(ind);
1719
1720     // this is utterly slow... (can be optmized if really needed)
1721
1722     // get all possible constituents of the tower
1723     std::vector<DetId> allConstituents = theTowerConstituentsMap->constituentsOf(id);
1724
1725     for (std::vector<DetId>::iterator ac_it = allConstituents.begin(); ac_it != allConstituents.end(); ++ac_it) {
1726       if (ac_it->det() != DetId::Ecal)
1727         continue;
1728
1729       auto thisEcalSevLvl = theEcalSevLvlAlgo->severityLevel(*ac_it);
1730
1731       // check if the Ecal severity is ok to keep
1732       std::vector<int>::const_iterator sevit =
1733           std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), thisEcalSevLvl);
1734       if (sevit != theEcalSeveritiesToBeExcluded.end()) {
1735         ++numBadEcalChan;
1736       }
1737     }
1738
1739     // if (0!=numBadEcalChan) std::cout << id << ":" << numBadEcalChan << ", ";
1740   }
1741
1742   /*
1743   int tot=0;
1744   for (auto ind=0U; ind<theTowerTopology->sizeForDenseIndexing(); ++ind) {
1745     if (ecalBadChs[ind]!=0) ++tot;
1746   }
1747   std::cout << " | " << tot << std::endl;
1748   */
1749 }
1750
1751 //////  Get status of the channel contributing to the tower
1752
1753 unsigned int CaloTowersCreationAlgo::hcalChanStatusForCaloTower(const CaloRecHit* hit) {
1754   HcalDetId hid(hit->detid());
1755   DetId id = theHcalTopology->idFront(hid);
1756 #ifdef EDM_ML_DEBUG
1757   std::cout << "ChanStatusForCaloTower for " << hid << " to " << HcalDetId(id) << std::endl;
1758 #endif
1759   const uint32_t recHitFlag = hit->flags();
1760   const uint32_t dbStatusFlag = theHcalChStatus->getValues(id)->getValue();
1761
1762   int severityLevel = theHcalSevLvlComputer->getSeverityLevel(id, recHitFlag, dbStatusFlag);
1763   bool isRecovered = theHcalSevLvlComputer->recoveredRecHit(id, recHitFlag);
1764
1765   // For use with hits rejected in the default reconstruction
1766   if (useRejectedHitsOnly) {
1767     if (!isRecovered) {
1768       if (severityLevel <= int(theHcalAcceptSeverityLevel) ||
1769           severityLevel > int(theHcalAcceptSeverityLevelForRejectedHit))
1770         return CaloTowersCreationAlgo::IgnoredChan;
1771       // this hit was either already accepted or is worse than
1772     } else {
1773       if (theRecoveredHcalHitsAreUsed || !useRejectedRecoveredHcalHits) {
1774         // skip recovered hits either because they were already used or because there was an explicit instruction
1775         return CaloTowersCreationAlgo::IgnoredChan;
1776       } else if (useRejectedRecoveredHcalHits) {
1777         return CaloTowersCreationAlgo::RecoveredChan;
1778       }
1779
1780     }  // recovered channels
1781
1782     // clasify channels as problematic: no good hits are supposed to be present in the
1783     // extra rechit collections
1784     return CaloTowersCreationAlgo::ProblematicChan;
1785
1786   }  // treatment of rejected hits
1787
1788   // this is for the regular reconstruction sequence
1789
1790   if (severityLevel == 0)
1791     return CaloTowersCreationAlgo::GoodChan;
1792
1793   if (isRecovered) {
1794     return (theRecoveredHcalHitsAreUsed) ? CaloTowersCreationAlgo::RecoveredChan : CaloTowersCreationAlgo::BadChan;
1795   } else {
1796     if (severityLevel > int(theHcalAcceptSeverityLevel)) {
1797       return CaloTowersCreationAlgo::BadChan;
1798     } else {
1799       return CaloTowersCreationAlgo::ProblematicChan;
1800     }
1801   }
1802 }
1803
1804 std::tuple<unsigned int, bool> CaloTowersCreationAlgo::ecalChanStatusForCaloTower(const EcalRecHit* hit) {
1805   // const DetId id = hit->detid();
1806
1807   //  uint16_t dbStatus = theEcalChStatus->find(id)->getStatusCode();
1808   //  uint32_t rhFlags  = hit->flags();
1809   //  int severityLevel = theEcalSevLvlAlgo->severityLevel(rhFlags, dbStatus);
1810   // The methods above will become private and cannot be usef for flagging ecal spikes.
1811   // Use the recommended interface - we leave the parameters for spilke removal to be specified by ECAL.
1812
1813   //  int severityLevel = 999;
1814
1815   EcalRecHit const& rh = *reinterpret_cast<EcalRecHit const*>(hit);
1816   int severityLevel = theEcalSevLvlAlgo->severityLevel(rh);
1817
1818   //  if      (id.subdetId() == EcalBarrel) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEbHandle);//, *theEcalChStatus);
1819   //  else if (id.subdetId() == EcalEndcap) severityLevel = theEcalSevLvlAlgo->severityLevel( id, *theEeHandle);//, *theEcalChStatus);
1820
1821   // there should be no other ECAL types used in this reconstruction
1822
1823   // The definition of ECAL severity levels uses categories that
1824   // are similar to the defined for CaloTower. (However, the categorization
1825   // for CaloTowers depends on the specified maximum acceptabel severity and therefore cannnot
1826   // be exact correspondence between the two. ECAL has additional categories describing modes of failure.)
1827   // This approach is different from the initial idea and from
1828   // the implementation for HCAL. Still make the logic similar to HCAL so that one has the ability to
1829   // exclude problematic channels as defined by ECAL.
1830   // For definitions of ECAL severity levels see RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h
1831
1832   bool isBad = (severityLevel == EcalSeverityLevel::kBad);
1833
1834   bool isRecovered = (severityLevel == EcalSeverityLevel::kRecovered);
1835
1836   // check if the severity is compatible with our configuration
1837   // This applies to the "default" tower cleaning
1838   std::vector<int>::const_iterator sevit =
1839       std::find(theEcalSeveritiesToBeExcluded.begin(), theEcalSeveritiesToBeExcluded.end(), severityLevel);
1840   bool accepted = (sevit == theEcalSeveritiesToBeExcluded.end());
1841
1842   // For use with hits that were rejected in the regular reconstruction:
1843   // This is for creating calotowers with lower level of cleaning by merging
1844   // the information from the default towers and a collection of towers created from
1845   // bad rechits
1846
1847   if (useRejectedHitsOnly) {
1848     if (!isRecovered) {
1849       if (accepted || std::find(theEcalSeveritiesToBeUsedInBadTowers.begin(),
1850                                 theEcalSeveritiesToBeUsedInBadTowers.end(),
1851                                 severityLevel) == theEcalSeveritiesToBeUsedInBadTowers.end())
1852         return std::make_tuple(CaloTowersCreationAlgo::IgnoredChan, isBad);
1853       // this hit was either already accepted, or is not eligible for inclusion
1854     } else {
1855       if (theRecoveredEcalHitsAreUsed || !useRejectedRecoveredEcalHits) {
1856         // skip recovered hits either because they were already used or because there was an explicit instruction
1857         return std::make_tuple(CaloTowersCreationAlgo::IgnoredChan, isBad);
1858         ;
1859       } else if (useRejectedRecoveredEcalHits) {
1860         return std::make_tuple(CaloTowersCreationAlgo::RecoveredChan, isBad);
1861       }
1862
1863     }  // recovered channels
1864
1865     // clasify channels as problematic
1866     return std::make_tuple(CaloTowersCreationAlgo::ProblematicChan, isBad);
1867
1868   }  // treatment of rejected hits
1869
1870   // for normal reconstruction
1871   if (severityLevel == EcalSeverityLevel::kGood)
1872     return std::make_tuple(CaloTowersCreationAlgo::GoodChan, false);
1873
1874   if (isRecovered) {
1875     return std::make_tuple(
1876         (theRecoveredEcalHitsAreUsed) ? CaloTowersCreationAlgo::RecoveredChan : CaloTowersCreationAlgo::BadChan, true);
1877   } else {
1878     return std::make_tuple(accepted ? CaloTowersCreationAlgo::ProblematicChan : CaloTowersCreationAlgo::BadChan, isBad);
1879   }
1880 }