EgammaCoreTools/interface/EcalClusterLazyTools.h

0001 #ifndef RecoEcal_EgammaCoreTools_EcalClusterLazyTools_h
0002 #define RecoEcal_EgammaCoreTools_EcalClusterLazyTools_h
0003
0004 /** \class EcalClusterLazyTools
0005  *
0006  * various cluster tools (e.g. cluster shapes)
0007  *
0008  * \author Federico Ferri
0009  *
0010  * \version $Id:
0011  *
0012  */
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0018
0019 #include "FWCore/Framework/interface/ESHandle.h"
0020
0021 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0022 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
0023 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0026 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0027 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0028
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0031 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0032 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0033 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
0034 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0035 #include "FWCore/Framework/interface/ConsumesCollector.h"
0036 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0037 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0038 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
0039 #include <optional>
0040
0041 class CaloTopology;
0042 class CaloGeometry;
0043 class CaloSubdetectorTopology;
0044
0045 class EcalClusterLazyToolsBase {
0046 public:
0047   struct ESData {
0048     CaloGeometry const &caloGeometry;
0049     CaloTopology const &caloTopology;
0050     EcalIntercalibConstants const &ecalIntercalibConstants;
0051     EcalADCToGeVConstant const &ecalADCToGeV;
0052     EcalLaserDbService const &ecalLaserDbService;
0053   };
0054
0055   class ESGetTokens {
0056   public:
0057     ESGetTokens(edm::ConsumesCollector cc)
0058         : caloGeometryToken_{cc.esConsumes()},
0059           caloTopologyToken_{cc.esConsumes()},
0060           ecalIntercalibConstantsToken_{cc.esConsumes()},
0061           ecalADCToGeVConstantToken_{cc.esConsumes()},
0062           ecalLaserDbServiceToken_{cc.esConsumes()} {}
0063
0064     ESData get(edm::EventSetup const &eventSetup) const {
0065       return {.caloGeometry = eventSetup.getData(caloGeometryToken_),
0066               .caloTopology = eventSetup.getData(caloTopologyToken_),
0067               .ecalIntercalibConstants = eventSetup.getData(ecalIntercalibConstantsToken_),
0068               .ecalADCToGeV = eventSetup.getData(ecalADCToGeVConstantToken_),
0069               .ecalLaserDbService = eventSetup.getData(ecalLaserDbServiceToken_)};
0070     }
0071
0072   private:
0073     edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0074     edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0075     edm::ESGetToken<EcalIntercalibConstants, EcalIntercalibConstantsRcd> ecalIntercalibConstantsToken_;
0076     edm::ESGetToken<EcalADCToGeVConstant, EcalADCToGeVConstantRcd> ecalADCToGeVConstantToken_;
0077     edm::ESGetToken<EcalLaserDbService, EcalLaserDbRecord> ecalLaserDbServiceToken_;
0078   };
0079
0080   EcalClusterLazyToolsBase(const edm::Event &ev,
0081                            ESData const &esData,
0082                            edm::EDGetTokenT<EcalRecHitCollection> token1,
0083                            edm::EDGetTokenT<EcalRecHitCollection> token2,
0084                            std::optional<edm::EDGetTokenT<EcalRecHitCollection>> token3);
0085
0086   // get time of basic cluster seed crystal
0087   float BasicClusterSeedTime(const reco::BasicCluster &cluster);
0088   // error-weighted average of time from constituents of basic cluster
0089   float BasicClusterTime(const reco::BasicCluster &cluster, const edm::Event &ev);
0090   // get BasicClusterSeedTime of the seed basic cluser of the supercluster
0091   float SuperClusterSeedTime(const reco::SuperCluster &cluster);
0092   // get BasicClusterTime of the seed basic cluser of the supercluster
0093   inline float SuperClusterTime(const reco::SuperCluster &cluster, const edm::Event &ev) {
0094     return BasicClusterTime((*cluster.seed()), ev);
0095   }
0096
0097   // mapping for preshower rechits
0098   std::map<DetId, EcalRecHit> rechits_map_;
0099   // get Preshower hit array
0100   std::vector<float> getESHits(double X,
0101                                double Y,
0102                                double Z,
0103                                const std::map<DetId, EcalRecHit> &rechits_map,
0104                                const CaloGeometry *geometry,
0105                                CaloSubdetectorTopology const *topology_p,
0106                                int row = 0,
0107                                int plane = 1);
0108   // get Preshower hit shape
0109   float getESShape(const std::vector<float> &ESHits0);
0110   // get Preshower effective sigmaRR
0111   float eseffsirir(const reco::SuperCluster &cluster);
0112   float eseffsixix(const reco::SuperCluster &cluster);
0113   float eseffsiyiy(const reco::SuperCluster &cluster);
0114
0115 protected:
0116   EcalRecHitCollection const *getEcalRecHitCollection(const reco::BasicCluster &cluster) const;
0117
0118   const CaloGeometry *geometry_;
0119   const CaloTopology *topology_;
0120   const EcalRecHitCollection *ebRecHits_;
0121   const EcalRecHitCollection *eeRecHits_;
0122   const EcalRecHitCollection *esRecHits_;
0123
0124   std::unique_ptr<CaloSubdetectorTopology const> ecalPS_topology_ = nullptr;
0125
0126   EcalIntercalibConstants const *ical = nullptr;
0127   EcalIntercalibConstantMap const *icalMap = nullptr;
0128   EcalADCToGeVConstant const *agc = nullptr;
0129   EcalLaserDbService const *laser = nullptr;
0130
0131 private:
0132   void getESRecHits(const edm::Event &ev, edm::EDGetTokenT<EcalRecHitCollection> const &esRecHitsToken);
0133
0134 };  // class EcalClusterLazyToolsBase
0135
0136 template <class ClusterTools>
0137 class EcalClusterLazyToolsT : public EcalClusterLazyToolsBase {
0138 public:
0139   EcalClusterLazyToolsT(const edm::Event &ev,
0140                         ESData const &esData,
0141                         edm::EDGetTokenT<EcalRecHitCollection> token1,
0142                         edm::EDGetTokenT<EcalRecHitCollection> token2)
0143       : EcalClusterLazyToolsBase(ev, esData, token1, token2, std::nullopt) {}
0144
0145   EcalClusterLazyToolsT(const edm::Event &ev,
0146                         ESData const &esData,
0147                         edm::EDGetTokenT<EcalRecHitCollection> token1,
0148                         edm::EDGetTokenT<EcalRecHitCollection> token2,
0149                         edm::EDGetTokenT<EcalRecHitCollection> token3)
0150       : EcalClusterLazyToolsBase(ev, esData, token1, token2, token3) {}
0151
0152   // Get the rec hit energies in a rectangle matrix around the seed.
0153   std::vector<float> energyMatrix(reco::BasicCluster const &cluster, int size) const {
0154     auto recHits = getEcalRecHitCollection(cluster);
0155     DetId maxId = ClusterTools::getMaximum(cluster, recHits).first;
0156
0157     std::vector<float> energies;
0158     for (auto const &detId : CaloRectangleRange(size, maxId, *topology_)) {
0159       energies.push_back(ClusterTools::recHitEnergy(detId, recHits));
0160     }
0161
0162     return energies;
0163   }
0164
0165   // various energies in the matrix nxn surrounding the maximum energy crystal of the input cluster
0166   //
0167   // NOTE (29/10/08): we now use an eta/phi coordinate system rather than
0168   //                  phi/eta to minmise possible screwups, for now e5x1 isnt
0169   //                  defined all the majority of people who call it actually
0170   //                  want e1x5 and it is thought it is better that their code
0171   //                  doesnt compile rather than pick up the wrong function
0172   //                  therefore in this version and later e1x5 = e5x1 in the old
0173   //                  version so 1x5 is 1 crystal in eta and 5 crystals in phi
0174   //                  note e3x2 does not have a definate eta/phi geometry, it
0175   //                  takes the maximum 3x2 block containing the seed regardless
0176   //                  of whether that 3 in eta or phi
0177   float e1x3(const reco::BasicCluster &cluster) const {
0178     return ClusterTools::e1x3(cluster, getEcalRecHitCollection(cluster), topology_);
0179   }
0180   float e3x1(const reco::BasicCluster &cluster) const {
0181     return ClusterTools::e3x1(cluster, getEcalRecHitCollection(cluster), topology_);
0182   }
0183   float e1x5(const reco::BasicCluster &cluster) const {
0184     return ClusterTools::e1x5(cluster, getEcalRecHitCollection(cluster), topology_);
0185   }
0186   float e5x1(const reco::BasicCluster &cluster) const {
0187     return ClusterTools::e5x1(cluster, getEcalRecHitCollection(cluster), topology_);
0188   }
0189   float e2x2(const reco::BasicCluster &cluster) const {
0190     return ClusterTools::e2x2(cluster, getEcalRecHitCollection(cluster), topology_);
0191   }
0192   float e3x2(const reco::BasicCluster &cluster) const {
0193     return ClusterTools::e3x2(cluster, getEcalRecHitCollection(cluster), topology_);
0194   }
0195   float e3x3(const reco::BasicCluster &cluster) const {
0196     return ClusterTools::e3x3(cluster, getEcalRecHitCollection(cluster), topology_);
0197   }
0198   float e4x4(const reco::BasicCluster &cluster) const {
0199     return ClusterTools::e4x4(cluster, getEcalRecHitCollection(cluster), topology_);
0200   }
0201   float e5x5(const reco::BasicCluster &cluster) const {
0202     return ClusterTools::e5x5(cluster, getEcalRecHitCollection(cluster), topology_);
0203   }
0204   int n5x5(const reco::BasicCluster &cluster) const {
0205     return ClusterTools::n5x5(cluster, getEcalRecHitCollection(cluster), topology_);
0206   }
0207   // energy in the 2x5 strip right of the max crystal (does not contain max crystal)
0208   // 2 crystals wide in eta, 5 wide in phi.
0209   float e2x5Right(const reco::BasicCluster &cluster) const {
0210     return ClusterTools::e2x5Right(cluster, getEcalRecHitCollection(cluster), topology_);
0211   }
0212   // energy in the 2x5 strip left of the max crystal (does not contain max crystal)
0213   float e2x5Left(const reco::BasicCluster &cluster) const {
0214     return ClusterTools::e2x5Left(cluster, getEcalRecHitCollection(cluster), topology_);
0215   }
0216   // energy in the 5x2 strip above the max crystal (does not contain max crystal)
0217   // 5 crystals wide in eta, 2 wide in phi.
0218   float e2x5Top(const reco::BasicCluster &cluster) const {
0219     return ClusterTools::e2x5Top(cluster, getEcalRecHitCollection(cluster), topology_);
0220   }
0221
0222   // energy in the 5x2 strip below the max crystal (does not contain max crystal)
0223   float e2x5Bottom(const reco::BasicCluster &cluster) const {
0224     return ClusterTools::e2x5Bottom(cluster, getEcalRecHitCollection(cluster), topology_);
0225   }
0226   // energy in a 2x5 strip containing the seed (max) crystal.
0227   // 2 crystals wide in eta, 5 wide in phi.
0228   // it is the maximum of either (1x5left + 1x5center) or (1x5right + 1x5center)
0229   float e2x5Max(const reco::BasicCluster &cluster) const {
0230     return ClusterTools::e2x5Max(cluster, getEcalRecHitCollection(cluster), topology_);
0231   }
0232   // energies in the crystal left, right, top, bottom w.r.t. to the most energetic crystal
0233   float eLeft(const reco::BasicCluster &cluster) const {
0234     return ClusterTools::eLeft(cluster, getEcalRecHitCollection(cluster), topology_);
0235   }
0236   float eRight(const reco::BasicCluster &cluster) const {
0237     return ClusterTools::eRight(cluster, getEcalRecHitCollection(cluster), topology_);
0238   }
0239   float eTop(const reco::BasicCluster &cluster) const {
0240     return ClusterTools::eTop(cluster, getEcalRecHitCollection(cluster), topology_);
0241   }
0242   float eBottom(const reco::BasicCluster &cluster) const {
0243     return ClusterTools::eBottom(cluster, getEcalRecHitCollection(cluster), topology_);
0244   }
0245   // the energy of the most energetic crystal in the cluster
0246   float eMax(const reco::BasicCluster &cluster) const {
0247     return ClusterTools::eMax(cluster, getEcalRecHitCollection(cluster));
0248   }
0249   // the energy of the second most energetic crystal in the cluster
0250   float e2nd(const reco::BasicCluster &cluster) const {
0251     return ClusterTools::e2nd(cluster, getEcalRecHitCollection(cluster));
0252   }
0253
0254   // get the DetId and the energy of the maximum energy crystal of the input cluster
0255   std::pair<DetId, float> getMaximum(const reco::BasicCluster &cluster) const {
0256     return ClusterTools::getMaximum(cluster, getEcalRecHitCollection(cluster));
0257   }
0258   std::vector<float> energyBasketFractionEta(const reco::BasicCluster &cluster) const {
0259     return ClusterTools::energyBasketFractionEta(cluster, getEcalRecHitCollection(cluster));
0260   }
0261   std::vector<float> energyBasketFractionPhi(const reco::BasicCluster &cluster) const {
0262     return ClusterTools::energyBasketFractionPhi(cluster, getEcalRecHitCollection(cluster));
0263   }
0264
0265   // return a vector v with v[0] = etaLat, v[1] = phiLat, v[2] = lat
0266   std::vector<float> lat(const reco::BasicCluster &cluster, bool logW = true, float w0 = 4.7) const {
0267     return ClusterTools::lat(cluster, getEcalRecHitCollection(cluster), geometry_, logW, w0);
0268   }
0269
0270   // return an array v with v[0] = covEtaEta, v[1] = covEtaPhi, v[2] = covPhiPhi
0271   std::array<float, 3> covariances(const reco::BasicCluster &cluster, float w0 = 4.7) const {
0272     return ClusterTools::covariances(cluster, getEcalRecHitCollection(cluster), topology_, geometry_, w0);
0273   }
0274
0275   // return an array v with v[0] = covIEtaIEta, v[1] = covIEtaIPhi, v[2] = covIPhiIPhi
0276   // this function calculates differences in eta/phi in units of crystals not
0277   // global eta/phi this is gives better performance in the crack regions of
0278   // the calorimeter but gives otherwise identical results to covariances
0279   // function this is only defined for the barrel, it returns covariances when
0280   // the cluster is in the endcap Warning: covIEtaIEta has been studied by
0281   // egamma, but so far covIPhiIPhi hasnt been studied extensively so there
0282   // could be a bug in the covIPhiIEta or covIPhiIPhi calculations. I dont
0283   // think there is but as it hasnt been heavily used, there might be one
0284   std::array<float, 3> localCovariances(const reco::BasicCluster &cluster,
0285                                         float w0 = EgammaLocalCovParamDefaults::kRelEnCut,
0286                                         const EcalPFRecHitThresholds *rhthresholds = nullptr,
0287                                         float multEB = 0.0,
0288                                         float multEE = 0.0) const {
0289     return ClusterTools::localCovariances(
0290         cluster, getEcalRecHitCollection(cluster), topology_, w0, rhthresholds, multEB, multEE);
0291   }
0292   std::array<float, 3> scLocalCovariances(const reco::SuperCluster &cluster, float w0 = 4.7) const {
0293     return ClusterTools::scLocalCovariances(cluster, getEcalRecHitCollection(cluster), topology_, w0);
0294   }
0295   double zernike20(const reco::BasicCluster &cluster, double R0 = 6.6, bool logW = true, float w0 = 4.7) const {
0296     return ClusterTools::zernike20(cluster, getEcalRecHitCollection(cluster), geometry_, R0, logW, w0);
0297   }
0298   double zernike42(const reco::BasicCluster &cluster, double R0 = 6.6, bool logW = true, float w0 = 4.7) const {
0299     return ClusterTools::zernike42(cluster, getEcalRecHitCollection(cluster), geometry_, R0, logW, w0);
0300   }
0301
0302 };  // class EcalClusterLazyToolsT
0303
0304 namespace noZS {
0305   typedef EcalClusterLazyToolsT<noZS::EcalClusterTools> EcalClusterLazyTools;
0306 }
0307 typedef EcalClusterLazyToolsT<::EcalClusterTools> EcalClusterLazyTools;
0308
0309 #endif