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