File indexing completed on 2024-04-06 12:24:40
0001 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0004 #include "EcalBasicClusterLocalContCorrection.h"
0005
0006 #include "TVector2.h"
0007
0008 EcalBasicClusterLocalContCorrection::EcalBasicClusterLocalContCorrection(edm::ConsumesCollector &&cc)
0009 : paramsToken_{cc.esConsumes()}, caloGeometryToken_{cc.esConsumes()} {}
0010
0011 void EcalBasicClusterLocalContCorrection::init(const edm::EventSetup &es) {
0012 params_ = &es.getData(paramsToken_);
0013 caloGeometry_ = &es.getData(caloGeometryToken_);
0014 }
0015
0016 void EcalBasicClusterLocalContCorrection::checkInit() const {
0017 if (!params_) {
0018
0019 throw cms::Exception("EcalBasicClusterLocalContCorrection::checkInit()")
0020 << "Trying to access an uninitialized crack correction function.\n"
0021 "Please call `init( edm::EventSetup &)' before any use of the function.\n";
0022 }
0023 }
0024
0025 using namespace std;
0026 using namespace edm;
0027
0028 float EcalBasicClusterLocalContCorrection::operator()(const reco::BasicCluster &basicCluster,
0029 const EcalRecHitCollection &recHit) const {
0030 checkInit();
0031
0032
0033 constexpr size_t nparams = 24;
0034
0035
0036 double correction_factor = 1.;
0037 double fetacor = 1.;
0038 double fphicor = 1.;
0039
0040
0041 const CaloSubdetectorGeometry *geom =
0042 caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0043
0044 const math::XYZPoint &position_ = basicCluster.position();
0045 double Theta = -position_.theta() + 0.5 * M_PI;
0046 double Eta = position_.eta();
0047 double Phi = TVector2::Phi_mpi_pi(position_.phi());
0048
0049
0050
0051 constexpr float X0 = 0.89;
0052 constexpr float T0 = 7.4;
0053 double depth = X0 * (T0 + log(basicCluster.energy()));
0054
0055
0056
0057 std::vector<std::pair<DetId, float> > crystals_vector = basicCluster.hitsAndFractions();
0058 float dphimin = 999.;
0059 float detamin = 999.;
0060 int ietaclosest = 0;
0061 int iphiclosest = 0;
0062
0063 for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
0064 EBDetId crystal(crystals_vector[icry].first);
0065 auto cell = geom->getGeometry(crystal);
0066 GlobalPoint center_pos = cell->getPosition(depth);
0067 double EtaCentr = center_pos.eta();
0068 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0069 if (std::abs(EtaCentr - Eta) < detamin) {
0070 detamin = std::abs(EtaCentr - Eta);
0071 ietaclosest = crystal.ieta();
0072 }
0073 if (std::abs(TVector2::Phi_mpi_pi(PhiCentr - Phi)) < dphimin) {
0074 dphimin = std::abs(TVector2::Phi_mpi_pi(PhiCentr - Phi));
0075 iphiclosest = crystal.iphi();
0076 }
0077 }
0078
0079 EBDetId crystalseed(ietaclosest, iphiclosest);
0080
0081
0082 auto cell = geom->getGeometry(crystalseed);
0083 GlobalPoint center_pos = cell->getPosition(depth);
0084
0085
0086 double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0087 double PhiWidth = (M_PI / 180.);
0088 double PhiCry = (TVector2::Phi_mpi_pi(Phi - PhiCentr)) / PhiWidth;
0089 if (PhiCry > 0.5)
0090 PhiCry = 0.5;
0091 if (PhiCry < -0.5)
0092 PhiCry = -0.5;
0093
0094 if (ietaclosest < 0)
0095 PhiCry *= -1.;
0096
0097
0098 double ThetaCentr = -center_pos.theta() + 0.5 * M_PI;
0099 double ThetaWidth = (M_PI / 180.) * std::cos(ThetaCentr);
0100 double EtaCry = (Theta - ThetaCentr) / ThetaWidth;
0101 if (EtaCry > 0.5)
0102 EtaCry = 0.5;
0103 if (EtaCry < -0.5)
0104 EtaCry = -0.5;
0105
0106 if (ietaclosest < 0)
0107 EtaCry *= -1.;
0108
0109
0110
0111 size_t payloadsize = params_->params().size();
0112
0113 if (payloadsize < nparams)
0114 edm::LogError("Invalid Payload") << "Parametrization requires " << nparams << " parameters but only " << payloadsize
0115 << " are found in DB. Perhaps incompatible Global Tag" << std::endl;
0116
0117 if (payloadsize > nparams)
0118 edm::LogWarning("Size mismatch ") << "Parametrization requires " << nparams << " parameters but " << payloadsize
0119 << " are found in DB. Perhaps incompatible Global Tag" << std::endl;
0120
0121 std::pair<double, double> localPosition(EtaCry, PhiCry);
0122
0123
0124 float localEta = localPosition.first;
0125 float localPhi = localPosition.second;
0126
0127
0128 int imod = getEcalModule(basicCluster.seed());
0129
0130
0131 float pe[3], pp[3];
0132 pe[0] = (params_->params())[0 + imod * 3];
0133 pe[1] = (params_->params())[1 + imod * 3];
0134 pe[2] = (params_->params())[2 + imod * 3];
0135 pp[0] = (params_->params())[12 + imod * 3];
0136 pp[1] = (params_->params())[13 + imod * 3];
0137 pp[2] = (params_->params())[14 + imod * 3];
0138
0139
0140 fetacor = pe[0] + pe[1] * localEta + pe[2] * localEta * localEta;
0141
0142
0143 fphicor = pp[0] + pp[1] * localPhi + pp[2] * localPhi * localPhi;
0144
0145
0146 int iphimod20 = std::abs(iphiclosest % 20);
0147 if (iphimod20 <= 1)
0148 fphicor = 1.;
0149
0150 correction_factor = (1. / fetacor) * (1. / fphicor);
0151
0152
0153 return correction_factor;
0154 }
0155
0156
0157 int EcalBasicClusterLocalContCorrection::getEcalModule(DetId id) const {
0158 int mod = 0;
0159 int ieta = (EBDetId(id)).ieta();
0160
0161 if (fabs(ieta) <= 25)
0162 mod = 0;
0163 if (fabs(ieta) > 25 && fabs(ieta) <= 45)
0164 mod = 1;
0165 if (fabs(ieta) > 45 && fabs(ieta) <= 65)
0166 mod = 2;
0167 if (fabs(ieta) > 65 && fabs(ieta) <= 85)
0168 mod = 3;
0169
0170 return (mod);
0171 }