CaloTowersCreator/src/CaloTowersCreationAlgo.h

0001 #ifndef RECOLOCALCALO_CALOTOWERSCREATOR_CALOTOWERSCREATIONALGO_H
0002 #define RECOLOCALCALO_CALOTOWERSCREATOR_CALOTOWERSCREATIONALGO_H 1
0003
0004 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0005 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0006 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0007 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0008 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0009
0010 // channel status
0011 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0012 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0013
0014 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0015 #include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
0016
0017 // severity level assignment for HCAL
0018 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0019 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0020
0021 // severity level assignment for ECAL
0022 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0023
0024 // need if we want to store the handles
0025 #include "FWCore/Framework/interface/ESHandle.h"
0026 #include <tuple>
0027
0028 #include <map>
0029
0030 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0031 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0032
0033 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0034
0035 class CaloTowerTopology;
0036 class HcalTopology;
0037 class CaloGeometry;
0038 class CaloSubdetectorGeometry;
0039 class CaloTowerConstituentsMap;
0040 class CaloRecHit;
0041 class DetId;
0042
0043 /** \class CaloTowersCreationAlgo
0044   *
0045   * \author R. Wilkinson - Caltech
0046   */
0047
0048 //
0049 // Modify MetaTower to save energy of rechits for use in tower 4-momentum assignment,
0050 // added containers for timing assignment and for holding status information.
0051 // Anton Anastassov (Northwestern)
0052 //
0053
0054 class CaloTowersCreationAlgo {
0055 public:
0056   int nalgo = -1;
0057
0058   CaloTowersCreationAlgo();
0059
0060   CaloTowersCreationAlgo(double EBthreshold,
0061                          double EEthreshold,
0062
0063                          bool useEtEBTreshold,
0064                          bool useEtEETreshold,
0065                          bool useSymEBTreshold,
0066                          bool useSymEETreshold,
0067
0068                          double HcalThreshold,
0069                          double HBthreshold,
0070                          double HBthreshold1,
0071                          double HBthreshold2,
0072                          double HESthreshold,
0073                          double HESthreshold1,
0074                          double HEDthreshold,
0075                          double HEDthreshold1,
0076                          double HOthreshold0,
0077                          double HOthresholdPlus1,
0078                          double HOthresholdMinus1,
0079                          double HOthresholdPlus2,
0080                          double HOthresholdMinus2,
0081                          double HF1threshold,
0082                          double HF2threshold,
0083                          double EBweight,
0084                          double EEweight,
0085                          double HBweight,
0086                          double HESweight,
0087                          double HEDweight,
0088                          double HOweight,
0089                          double HF1weight,
0090                          double HF2weight,
0091                          double EcutTower,
0092                          double EBSumThreshold,
0093                          double EESumThreshold,
0094                          bool useHO,
0095                          // (for momentum reconstruction algorithm)
0096                          int momConstrMethod,
0097                          double momHBDepth,
0098                          double momHEDepth,
0099                          double momEBDepth,
0100                          double momEEDepth,
0101                          int hcalPhase = 0);
0102
0103   CaloTowersCreationAlgo(double EBthreshold,
0104                          double EEthreshold,
0105
0106                          bool useEtEBTreshold,
0107                          bool useEtEETreshold,
0108                          bool useSymEBTreshold,
0109                          bool useSymEETreshold,
0110
0111                          double HcalThreshold,
0112                          double HBthreshold,
0113                          double HBthreshold1,
0114                          double HBthreshold2,
0115                          double HESthreshold,
0116                          double HESthreshold1,
0117                          double HEDthreshold,
0118                          double HEDthreshold1,
0119                          double HOthreshold0,
0120                          double HOthresholdPlus1,
0121                          double HOthresholdMinus1,
0122                          double HOthresholdPlus2,
0123                          double HOthresholdMinus2,
0124                          double HF1threshold,
0125                          double HF2threshold,
0126                          const std::vector<double>& EBGrid,
0127                          const std::vector<double>& EBWeights,
0128                          const std::vector<double>& EEGrid,
0129                          const std::vector<double>& EEWeights,
0130                          const std::vector<double>& HBGrid,
0131                          const std::vector<double>& HBWeights,
0132                          const std::vector<double>& HESGrid,
0133                          const std::vector<double>& HESWeights,
0134                          const std::vector<double>& HEDGrid,
0135                          const std::vector<double>& HEDWeights,
0136                          const std::vector<double>& HOGrid,
0137                          const std::vector<double>& HOWeights,
0138                          const std::vector<double>& HF1Grid,
0139                          const std::vector<double>& HF1Weights,
0140                          const std::vector<double>& HF2Grid,
0141                          const std::vector<double>& HF2Weights,
0142                          double EBweight,
0143                          double EEweight,
0144                          double HBweight,
0145                          double HESweight,
0146                          double HEDweight,
0147                          double HOweight,
0148                          double HF1weight,
0149                          double HF2weight,
0150                          double EcutTower,
0151                          double EBSumThreshold,
0152                          double EESumThreshold,
0153                          bool useHO,
0154                          // (for momentum reconstruction algorithm)
0155                          int momConstrMethod,
0156                          double momHBDepth,
0157                          double momHEDepth,
0158                          double momEBDepth,
0159                          double momEEDepth,
0160                          int hcalPhase = 0);
0161
0162   void setGeometry(const CaloTowerTopology* cttopo,
0163                    const CaloTowerConstituentsMap* ctmap,
0164                    const HcalTopology* htopo,
0165                    const CaloGeometry* geo);
0166
0167   void setThresFromDB(const EcalPFRecHitThresholds* EcalCuts, const HcalPFCuts* HcalCuts);
0168   // pass the containers of channels status from the event record (stored in DB)
0169   // these are called in  CaloTowersCreator
0170   void setHcalChStatusFromDB(const HcalChannelQuality* s) { theHcalChStatus = s; }
0171   void setEcalChStatusFromDB(const EcalChannelStatus* s) { theEcalChStatus = s; }
0172
0173   // Kake a map of number of channels not used in RecHit production.
0174   // The key is the calotower id.
0175   void makeHcalDropChMap();
0176
0177   void makeEcalBadChs();
0178
0179   void begin();
0180   void process(const HBHERecHitCollection& hbhe);
0181   void process(const HORecHitCollection& ho);
0182   void process(const HFRecHitCollection& hf);
0183   void process(const EcalRecHitCollection& ecal);
0184
0185   void process(const CaloTowerCollection& ctc);
0186
0187   void finish(CaloTowerCollection& destCollection);
0188
0189   // modified rescale method
0190   void rescaleTowers(const CaloTowerCollection& ctInput, CaloTowerCollection& ctResult);
0191
0192   void setEBEScale(double scale);
0193   void setEEEScale(double scale);
0194   void setHBEScale(double scale);
0195   void setHESEScale(double scale);
0196   void setHEDEScale(double scale);
0197   void setHOEScale(double scale);
0198   void setHF1EScale(double scale);
0199   void setHF2EScale(double scale);
0200
0201   // Assign to categories based on info from DB and RecHit status
0202   // Called in assignHit to check if the energy should be added to
0203   // calotower, and how to flag the channel
0204   unsigned int hcalChanStatusForCaloTower(const CaloRecHit* hit);
0205   std::tuple<unsigned int, bool> ecalChanStatusForCaloTower(const EcalRecHit* hit);
0206
0207   // Channel flagging is based on acceptable severity levels specified in the
0208   // configuration file. These methods are used to pass the values read in
0209   // CaloTowersCreator
0210   //
0211   // from DB
0212   void setHcalAcceptSeverityLevel(unsigned int level) { theHcalAcceptSeverityLevel = level; }
0213   void setEcalSeveritiesToBeExcluded(const std::vector<int>& ecalSev) { theEcalSeveritiesToBeExcluded = ecalSev; }
0214
0215   // flag to use recovered hits
0216   void setRecoveredHcalHitsAreUsed(bool flag) { theRecoveredHcalHitsAreUsed = flag; };
0217   void setRecoveredEcalHitsAreUsed(bool flag) { theRecoveredEcalHitsAreUsed = flag; };
0218
0219   //  severety level calculator for HCAL
0220   void setHcalSevLvlComputer(const HcalSeverityLevelComputer* c) { theHcalSevLvlComputer = c; };
0221
0222   // severity level calculator for ECAL
0223   void setEcalSevLvlAlgo(const EcalSeverityLevelAlgo* a) { theEcalSevLvlAlgo = a; }
0224
0225   // The following are needed for creating towers from rechits excluded from the  ------------------------------------
0226   // default reconstructions
0227
0228   // NB! Controls if rejected hits shold be used instead of the default!!!
0229   void setUseRejectedHitsOnly(bool flag) { useRejectedHitsOnly = flag; }
0230
0231   void setHcalAcceptSeverityLevelForRejectedHit(unsigned int level) {
0232     theHcalAcceptSeverityLevelForRejectedHit = level;
0233   }
0234   //  void setEcalAcceptSeverityLevelForRejectedHit(unsigned int level) {theEcalAcceptSeverityLevelForRejectedHit = level;}
0235   void SetEcalSeveritiesToBeUsedInBadTowers(const std::vector<int>& ecalSev) {
0236     theEcalSeveritiesToBeUsedInBadTowers = ecalSev;
0237   }
0238
0239   void setUseRejectedRecoveredHcalHits(bool flag) { useRejectedRecoveredHcalHits = flag; };
0240   void setUseRejectedRecoveredEcalHits(bool flag) { useRejectedRecoveredEcalHits = flag; };
0241   void setMissingHcalRescaleFactorForEcal(float factor) { missingHcalRescaleFactorForEcal = factor; };
0242
0243   //-------------------------------------------------------------------------------------------------------------------
0244
0245   // set the EE EB handles
0246
0247   void setEbHandle(const edm::Handle<EcalRecHitCollection> eb) { theEbHandle = eb; }
0248   void setEeHandle(const edm::Handle<EcalRecHitCollection> ee) { theEeHandle = ee; }
0249
0250   // Add methods to get the seperate positions for ECAL/HCAL
0251   // used in constructing the 4-vectors using new methods
0252   GlobalPoint emCrystalShwrPos(DetId detId, float fracDepth);
0253   GlobalPoint hadSegmentShwrPos(DetId detId, float fracDepth);
0254   // "effective" point for the EM/HAD shower in CaloTower
0255   //  position based on non-zero energy cells
0256   GlobalPoint hadShwrPos(const std::vector<std::pair<DetId, float>>& metaContains, float fracDepth, double hadE);
0257   GlobalPoint emShwrPos(const std::vector<std::pair<DetId, float>>& metaContains, float fracDepth, double totEmE);
0258
0259   // overloaded function to get had position based on all had cells in the tower
0260   GlobalPoint hadShwrPos(CaloTowerDetId id, float fracDepth);
0261   GlobalPoint hadShwPosFromCells(DetId frontCell, DetId backCell, float fracDepth);
0262
0263   // for Chris
0264   GlobalPoint emShwrLogWeightPos(const std::vector<std::pair<DetId, float>>& metaContains,
0265                                  float fracDepth,
0266                                  double totEmE);
0267
0268 private:
0269   struct MetaTower {
0270     MetaTower() {}
0271     bool empty() const { return metaConstituents.empty(); }
0272     // contains also energy of RecHit
0273     std::vector<std::pair<DetId, float>> metaConstituents;
0274     CaloTowerDetId id;
0275     float E = 0, E_em = 0, E_had = 0, E_outer = 0;
0276     float emSumTimeTimesE = 0, hadSumTimeTimesE = 0, emSumEForTime = 0,
0277           hadSumEForTime = 0;  // Sum(Energy x Timing) : intermediate container
0278
0279     // needed to set CaloTower status word
0280     int numBadEcalCells = 0, numRecEcalCells = 0, numProbEcalCells = 0, numBadHcalCells = 0, numRecHcalCells = 0,
0281         numProbHcalCells = 0;
0282   };
0283
0284   /// adds a single hit to the tower
0285   void assignHitEcal(const EcalRecHit* recHit);
0286   void assignHitHcal(const CaloRecHit* recHit);
0287
0288   void rescale(const CaloTower* ct);
0289
0290   /// looks for a given tower in the internal cache.  If it can't find it, it makes it.
0291   MetaTower& find(const CaloTowerDetId& id);
0292
0293   /// helper method to look up the appropriate threshold & weight
0294   void getThresholdAndWeight(const DetId& detId, double& threshold, double& weight) const;
0295
0296   double theEBthreshold, theEEthreshold;
0297   bool theUseEtEBTresholdFlag, theUseEtEETresholdFlag;
0298   bool theUseSymEBTresholdFlag, theUseSymEETresholdFlag;
0299
0300   double theHcalThreshold;
0301
0302   double theHBthreshold, theHBthreshold1, theHBthreshold2;
0303   double theHESthreshold, theHESthreshold1;
0304   double theHEDthreshold, theHEDthreshold1;
0305   double theHOthreshold0, theHOthresholdPlus1, theHOthresholdMinus1;
0306   double theHOthresholdPlus2, theHOthresholdMinus2, theHF1threshold, theHF2threshold;
0307   std::vector<double> theEBGrid, theEBWeights;
0308   std::vector<double> theEEGrid, theEEWeights;
0309   std::vector<double> theHBGrid, theHBWeights;
0310   std::vector<double> theHESGrid, theHESWeights;
0311   std::vector<double> theHEDGrid, theHEDWeights;
0312   std::vector<double> theHOGrid, theHOWeights;
0313   std::vector<double> theHF1Grid, theHF1Weights;
0314   std::vector<double> theHF2Grid, theHF2Weights;
0315   double theEBweight, theEEweight;
0316   double theHBweight, theHESweight, theHEDweight, theHOweight, theHF1weight, theHF2weight;
0317   double theEcutTower, theEBSumThreshold, theEESumThreshold;
0318
0319   double theEBEScale;
0320   double theEEEScale;
0321   double theHBEScale;
0322   double theHESEScale;
0323   double theHEDEScale;
0324   double theHOEScale;
0325   double theHF1EScale;
0326   double theHF2EScale;
0327   const EcalPFRecHitThresholds* ecalCuts;
0328   const HcalPFCuts* hcalCuts;
0329   const CaloTowerTopology* theTowerTopology;
0330   const HcalTopology* theHcalTopology;
0331   const CaloGeometry* theGeometry;
0332   const CaloTowerConstituentsMap* theTowerConstituentsMap;
0333   const CaloSubdetectorGeometry* theTowerGeometry;
0334
0335   // for checking the status of ECAL and HCAL channels stored in the DB
0336   const EcalChannelStatus* theEcalChStatus;
0337   const HcalChannelQuality* theHcalChStatus;
0338
0339   // calculator of severety level for HCAL
0340   const HcalSeverityLevelComputer* theHcalSevLvlComputer;
0341
0342   // calculator for severity level for ECAL
0343   const EcalSeverityLevelAlgo* theEcalSevLvlAlgo;
0344
0345   // fields that hold the information passed from the CaloTowersCreator configuration file:
0346   // controll what is considered bad/recovered/problematic channel for CaloTower purposes
0347   //
0348   unsigned int theHcalAcceptSeverityLevel;
0349   std::vector<int> theEcalSeveritiesToBeExcluded;
0350   // flag to use recovered hits
0351   bool theRecoveredHcalHitsAreUsed;
0352   bool theRecoveredEcalHitsAreUsed;
0353
0354   // controls the tower reconstruction from rejected hits
0355
0356   bool useRejectedHitsOnly;
0357   unsigned int theHcalAcceptSeverityLevelForRejectedHit;
0358   std::vector<int> theEcalSeveritiesToBeUsedInBadTowers;
0359
0360   unsigned int useRejectedRecoveredHcalHits;
0361   unsigned int useRejectedRecoveredEcalHits;
0362
0363   // if Hcal is missing, fudge it scaling Ecal by this factor (if > 0)
0364   float missingHcalRescaleFactorForEcal;
0365
0366   /// only affects energy and ET calculation.  HO is still recorded in the tower
0367   bool theHOIsUsed;
0368
0369   // Switches and paramters for CaloTower 4-momentum assignment
0370   // "depth" variables do not affect all algorithms
0371   int theMomConstrMethod;
0372   double theMomHBDepth;
0373   double theMomHEDepth;
0374   double theMomEBDepth;
0375   double theMomEEDepth;
0376
0377   // compactify timing info
0378   int compactTime(float time);
0379
0380   void convert(const CaloTowerDetId& id, const MetaTower& mt, CaloTowerCollection& collection);
0381
0382   // internal map
0383   typedef std::vector<MetaTower> MetaTowerMap;
0384   MetaTowerMap theTowerMap;
0385   unsigned int theTowerMapSize = 0;
0386
0387   // Number of channels in the tower that were not used in RecHit production (dead/off,...).
0388   // These channels are added to the other "bad" channels found in the recHit collection.
0389   typedef std::map<CaloTowerDetId, std::pair<short int, bool>> HcalDropChMap;
0390   HcalDropChMap hcalDropChMap;
0391
0392   // Number of bad Ecal channel in each tower
0393   //unsigned short ecalBadChs[CaloTowerDetId::kSizeForDenseIndexing];
0394   std::vector<unsigned short> ecalBadChs;
0395
0396   // clasification of channels in tower construction: the category definition is
0397   // affected by the setting in the configuration file
0398   //
0399   enum ctHitCategory { GoodChan = 0, BadChan = 1, RecoveredChan = 2, ProblematicChan = 3, IgnoredChan = 99 };
0400
0401   // the EE and EB collections for ecal anomalous cell info
0402
0403   edm::Handle<EcalRecHitCollection> theEbHandle;
0404   edm::Handle<EcalRecHitCollection> theEeHandle;
0405
0406   int theHcalPhase;
0407
0408   std::vector<HcalDetId> ids_;
0409 };
0410
0411 #endif