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