EgammaCoreTools/plugins/EcalClusterCrackCorrection.cc

0001 /** \class EcalClusterCrackCorrection
0002   *  Function to correct cluster for cracks in the calorimeter
0003   *
0004   *  $Id: EcalClusterCrackCorrection.h
0005   *  $Date:
0006   *  $Revision:
0007   *  \author Federico Ferri, CEA Saclay, November 2008
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   // get/set explicit methods for parameters
0030   const EcalClusterCrackCorrParameters *getParameters() const { return params_; }
0031   // check initialization
0032   void checkInit() const;
0033
0034   // compute the correction
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   // set parameters
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     // non initialized function parameters: throw exception
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   //correction factor to be returned, and to be calculated in this present function:
0068   double correction_factor = 1.;
0069   double fetacor = 1.;  //eta dependent part of the correction factor
0070   double fphicor = 1.;  //phi dependent part of the correction factor
0071
0072   //********************************************************************************************************************//
0073   //These ECAL barrel module and supermodule border corrections correct a photon energy for leakage outside a 5x5 crystal cluster. They  depend on the local position in the hit crystal. The hit crystal needs to be at the border of a barrel module. The local position coordinates, called later EtaCry and PhiCry in the code, are comprised between -0.5 and 0.5 and correspond to the distance between the photon supercluster position and the center of the hit crystal, expressed in number of  crystal widthes. The correction parameters (that should be filled in CalibCalorimetry/EcalTrivialCondModules/python/EcalTrivialCondRetriever_cfi.py) were calculated using simulaion and thus take into account the effect of the magnetic field. They  only apply to unconverted photons in the barrel, but a use for non brem electrons could be considered (not tested yet). For more details, cf the CMS internal note 2009-013 by S. Tourneur and C. Seez
0074
0075   //Beware: The user should make sure it only uses this correction factor for unconverted photons (or not breming electrons)
0076
0077   //const reco::CaloClusterPtr & seedbclus =  superCluster.seed();
0078
0079   //If not barrel, return 1:
0080   if (std::abs(seedbclus.eta()) > 1.4442)
0081     return 1.;
0082
0083   const CaloSubdetectorGeometry *geom = caloGeom_->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);  //EcalBarrel = 1
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   //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
0091   // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
0092   const float X0 = 0.89;
0093   const float T0 = 7.4;
0094   double depth = X0 * (T0 + log(seedbclus.energy()));
0095
0096   //search which crystal is closest to the cluster position and call it crystalseed:
0097   //std::vector<DetId> crystals_vector = seedbclus.getHitsByDetId();   //deprecated
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   // Get center cell position from shower depth
0121   auto cell = geom->getGeometry(crystalseed);
0122   GlobalPoint center_pos = cell->getPosition(depth);
0123
0124   //if the seed crystal isn't neighbourgh of a supermodule border, don't apply the phi dependent crack corrections, but use the smaller phi dependent local containment correction instead.
0125   if (ietaclosest < 0)
0126     iphiclosest = 361 - iphiclosest;  //inversion of phi 3 degree tilt
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     //flip to take into account ECAL barrel symmetries:
0140     if (ietaclosest < 0)
0141       PhiCry *= -1.;
0142
0143     //Fetching parameters of the polynomial (see  CMS IN-2009/013)
0144     double g[5];
0145     int offset = iphimod20 == 0 ? 10   //coefficients for one phi side of a SM
0146                                 : 15;  //coefficients for the other side
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   //if the seed crystal isn't neighbourgh of a module border, don't apply the eta dependent crack corrections, but use the smaller eta dependent local containment correction instead.
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     //flip to take into account ECAL barrel symmetries:
0169     if (ietaclosest < 0)
0170       EtaCry *= -1.;
0171
0172     //Fetching parameters of the polynomial (see  CMS IN-2009/013)
0173     double f[5];
0174     int offset = std::abs(ietamod20) == 5
0175                      ? 0   //coefficients for eta side of an intermodule gap closer to the interaction point
0176                      : 5;  //coefficients for the other eta side
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   //return the correction factor. Use it to multiply the cluster energy.
0189   return correction_factor;
0190 }
0191
0192 float EcalClusterCrackCorrection::getValue(const reco::SuperCluster &superCluster, const int mode) const {
0193   checkInit();
0194
0195   //********************************************************************************************************************//
0196   //These ECAL barrel module and supermodule border corrections correct a photon energy for leakage outside a 5x5 crystal cluster. They  depend on the local position in the hit crystal. The hit crystal needs to be at the border of a barrel module. The local position coordinates, called later EtaCry and PhiCry in the code, are comprised between -0.5 and 0.5 and correspond to the distance between the photon supercluster position and the center of the hit crystal, expressed in number of  crystal widthes. The correction parameters (that should be filled in CalibCalorimetry/EcalTrivialCondModules/python/EcalTrivialCondRetriever_cfi.py) were calculated using simulaion and thus take into account the effect of the magnetic field. They  only apply to unconverted photons in the barrel, but a use for non brem electrons could be considered (not tested yet). For more details, cf the CMS internal note 2009-013 by S. Tourneur and C. Seez
0197
0198   //Beware: The user should make sure it only uses this correction factor for unconverted photons (or not breming electrons)
0199
0200   return getValue(*(superCluster.seed()));
0201 }
0202
0203 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterFunctionFactory.h"
0204 DEFINE_EDM_PLUGIN(EcalClusterFunctionFactory, EcalClusterCrackCorrection, "EcalClusterCrackCorrection");