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