EgammaClusterProducers/src/EcalBasicClusterLocalContCorrection.cc

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     // non-initialized function parameters: throw exception
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   // number of parameters needed by this parametrization
0033   constexpr size_t nparams = 24;
0034
0035   //correction factor to be returned, and to be calculated in this present function:
0036   double correction_factor = 1.;
0037   double fetacor = 1.;  //eta dependent part of the correction factor
0038   double fphicor = 1.;  //phi dependent part of the correction factor
0039
0040   //--------------if barrel calculate local position wrt xtal center -------------------
0041   const CaloSubdetectorGeometry *geom =
0042       caloGeometry_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);  //EcalBarrel = 1
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   //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
0050   // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
0051   constexpr float X0 = 0.89;
0052   constexpr float T0 = 7.4;
0053   double depth = X0 * (T0 + log(basicCluster.energy()));
0054
0055   //search which crystal is closest to the cluster position and call it crystalseed:
0056   //std::vector<DetId> crystals_vector = *scRef.getHitsByDetId();   //deprecated
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);  // problema qui
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   // Get center cell position from shower depth
0082   auto cell = geom->getGeometry(crystalseed);
0083   GlobalPoint center_pos = cell->getPosition(depth);
0084
0085   //PHI
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   //flip to take into account ECAL barrel symmetries:
0094   if (ietaclosest < 0)
0095     PhiCry *= -1.;
0096
0097   //ETA
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   //flip to take into account ECAL barrel symmetries:
0106   if (ietaclosest < 0)
0107     EtaCry *= -1.;
0108
0109   //-------------- end calculate local position -------------
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   //--- local cluster coordinates
0124   float localEta = localPosition.first;
0125   float localPhi = localPosition.second;
0126
0127   //--- ecal module
0128   int imod = getEcalModule(basicCluster.seed());
0129
0130   //-- corrections parameters
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   //--- correction vs local eta
0140   fetacor = pe[0] + pe[1] * localEta + pe[2] * localEta * localEta;
0141
0142   //--- correction vs local phi
0143   fphicor = pp[0] + pp[1] * localPhi + pp[2] * localPhi * localPhi;
0144
0145   //if the seed crystal is neighbourgh of a supermodule border, don't apply the phi dependent  containment corrections, but use the larger crack corrections instead.
0146   int iphimod20 = std::abs(iphiclosest % 20);
0147   if (iphimod20 <= 1)
0148     fphicor = 1.;
0149
0150   correction_factor = (1. / fetacor) * (1. / fphicor);
0151
0152   //return the correction factor. Use it to multiply the cluster energy.
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 }