File indexing completed on 2023-10-25 09:59:03
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "CondFormats/DataRecord/interface/EcalClusterCrackCorrParametersRcd.h"
0011 #include "CondFormats/EcalObjects/interface/EcalClusterCrackCorrParameters.h"
0012 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0013 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0014 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0015 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0016 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionBaseClass.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020 #include "FWCore/Utilities/interface/ESGetToken.h"
0021
0022 #include "TVector2.h"
0023
0024 class EcalClusterCrackCorrection : public EcalClusterFunctionBaseClass {
0025 public:
0026 EcalClusterCrackCorrection(const edm::ParameterSet &, edm::ConsumesCollector iC)
0027 : paramsToken_(iC.esConsumes()), geomToken_(iC.esConsumes()) {}
0028
0029
0030 const EcalClusterCrackCorrParameters *getParameters() const { return params_; }
0031
0032 void checkInit() const;
0033
0034
0035 float getValue(const reco::BasicCluster &, const EcalRecHitCollection &) const override { return 1.f; }
0036 float getValue(const reco::SuperCluster &, const int mode) const override;
0037
0038 float getValue(const reco::CaloCluster &) const override;
0039
0040
0041 void init(const edm::EventSetup &es) override;
0042
0043 private:
0044 const edm::ESGetToken<EcalClusterCrackCorrParameters, EcalClusterCrackCorrParametersRcd> paramsToken_;
0045 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0046 const EcalClusterCrackCorrParameters *params_ = nullptr;
0047 const CaloGeometry *caloGeom_ = nullptr;
0048 };
0049
0050 void EcalClusterCrackCorrection::init(const edm::EventSetup &es) {
0051 params_ = &es.getData(paramsToken_);
0052 caloGeom_ = &es.getData(geomToken_);
0053 }
0054
0055 void EcalClusterCrackCorrection::checkInit() const {
0056 if (!params_) {
0057
0058 throw cms::Exception("EcalClusterCrackCorrection::checkInit()")
0059 << "Trying to access an uninitialized crack correction function.\n"
0060 "Please call `init( edm::EventSetup &)' before any use of the function.\n";
0061 }
0062 }
0063
0064 float EcalClusterCrackCorrection::getValue(const reco::CaloCluster &seedbclus) const {
0065 checkInit();
0066
0067
0068 double correction_factor = 1.;
0069 double fetacor = 1.;
0070 double fphicor = 1.;
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080 if (std::abs(seedbclus.eta()) > 1.4442)
0081 return 1.;
0082
0083 const CaloSubdetectorGeometry *geom = caloGeom_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0084
0085 const math::XYZPoint &position_ = seedbclus.position();
0086 double Theta = -position_.theta() + 0.5 * M_PI;
0087 double Eta = position_.eta();
0088 double Phi = TVector2::Phi_mpi_pi(position_.phi());
0089
0090
0091
0092 const float X0 = 0.89;
0093 const float T0 = 7.4;
0094 double depth = X0 * (T0 + log(seedbclus.energy()));
0095
0096
0097
0098 std::vector<std::pair<DetId, float> > crystals_vector = seedbclus.hitsAndFractions();
0099 float dphimin = 999.;
0100 float detamin = 999.;
0101 int ietaclosest = 0;
0102 int iphiclosest = 0;
0103 for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
0104 EBDetId crystal(crystals_vector[icry].first);
0105 auto cell = geom->getGeometry(crystal);
0106 GlobalPoint center_pos = cell->getPosition(depth);
0107 double EtaCentr = center_pos.eta();
0108 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0109 if (std::abs(EtaCentr - Eta) < detamin) {
0110 detamin = std::abs(EtaCentr - Eta);
0111 ietaclosest = crystal.ieta();
0112 }
0113 if (std::abs(TVector2::Phi_mpi_pi(PhiCentr - Phi)) < dphimin) {
0114 dphimin = std::abs(TVector2::Phi_mpi_pi(PhiCentr - Phi));
0115 iphiclosest = crystal.iphi();
0116 }
0117 }
0118 EBDetId crystalseed(ietaclosest, iphiclosest);
0119
0120
0121 auto cell = geom->getGeometry(crystalseed);
0122 GlobalPoint center_pos = cell->getPosition(depth);
0123
0124
0125 if (ietaclosest < 0)
0126 iphiclosest = 361 - iphiclosest;
0127 int iphimod20 = iphiclosest % 20;
0128 if (iphimod20 > 1)
0129 fphicor = 1.;
0130
0131 else {
0132 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0133 double PhiWidth = (M_PI / 180.);
0134 double PhiCry = (TVector2::Phi_mpi_pi(Phi - PhiCentr)) / PhiWidth;
0135 if (PhiCry > 0.5)
0136 PhiCry = 0.5;
0137 if (PhiCry < -0.5)
0138 PhiCry = -0.5;
0139
0140 if (ietaclosest < 0)
0141 PhiCry *= -1.;
0142
0143
0144 double g[5];
0145 int offset = iphimod20 == 0 ? 10
0146 : 15;
0147 for (int k = 0; k != 5; ++k)
0148 g[k] = (params_->params())[k + offset];
0149
0150 fphicor = 0.;
0151 for (int k = 0; k != 5; ++k)
0152 fphicor += g[k] * std::pow(PhiCry, k);
0153 }
0154
0155
0156 int ietamod20 = ietaclosest % 20;
0157 if (std::abs(ietaclosest) < 25 || (std::abs(ietamod20) != 5 && std::abs(ietamod20) != 6))
0158 fetacor = 1.;
0159
0160 else {
0161 double ThetaCentr = -center_pos.theta() + 0.5 * M_PI;
0162 double ThetaWidth = (M_PI / 180.) * std::cos(ThetaCentr);
0163 double EtaCry = (Theta - ThetaCentr) / ThetaWidth;
0164 if (EtaCry > 0.5)
0165 EtaCry = 0.5;
0166 if (EtaCry < -0.5)
0167 EtaCry = -0.5;
0168
0169 if (ietaclosest < 0)
0170 EtaCry *= -1.;
0171
0172
0173 double f[5];
0174 int offset = std::abs(ietamod20) == 5
0175 ? 0
0176 : 5;
0177 for (int k = 0; k != 5; ++k)
0178 f[k] = (params_->params())[k + offset];
0179
0180 fetacor = 0.;
0181 for (int k = 0; k != 5; ++k)
0182 fetacor += f[k] * std::pow(EtaCry, k);
0183 }
0184
0185 correction_factor = 1. / (fetacor * fphicor);
0186
0187
0188
0189 return correction_factor;
0190 }
0191
0192 float EcalClusterCrackCorrection::getValue(const reco::SuperCluster &superCluster, const int mode) const {
0193 checkInit();
0194
0195
0196
0197
0198
0199
0200 return getValue(*(superCluster.seed()));
0201 }
0202
0203 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
0204 DEFINE_EDM_PLUGIN(EcalClusterFunctionFactory, EcalClusterCrackCorrection, "EcalClusterCrackCorrection");