EgammaTools/src/EcalClusterLocal.cc

0001 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
0002 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0003 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0004 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0005 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0009 #include "TVector2.h"
0010
0011 namespace egammaTools {
0012
0013   void localEcalClusterCoordsEB(const reco::CaloCluster &bclus,
0014                                 const CaloGeometry &caloGeometry,
0015                                 float &etacry,
0016                                 float &phicry,
0017                                 int &ieta,
0018                                 int &iphi,
0019                                 float &thetatilt,
0020                                 float &phitilt) {
0021     assert(bclus.hitsAndFractions().at(0).first.subdetId() == EcalBarrel);
0022
0023     const CaloSubdetectorGeometry *geom =
0024         caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);  //EcalBarrel = 1
0025
0026     const math::XYZPoint &position_ = bclus.position();
0027     double Theta = -position_.theta() + 0.5 * M_PI;
0028     double Eta = position_.eta();
0029     double Phi = TVector2::Phi_mpi_pi(position_.phi());
0030
0031     //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
0032     // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
0033     const float X0 = 0.89;
0034     const float T0 = 7.4;
0035     double depth = X0 * (T0 + log(bclus.energy()));
0036
0037     //find max energy crystal
0038     std::vector<std::pair<DetId, float> > crystals_vector = bclus.hitsAndFractions();
0039     float drmin = 999.;
0040     EBDetId crystalseed;
0041     //printf("starting loop over crystals, etot = %5f:\n",bclus.energy());
0042     for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
0043       EBDetId crystal(crystals_vector[icry].first);
0044
0045       auto cell = geom->getGeometry(crystal);
0046       const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());
0047       GlobalPoint center_pos = cpyr->getPosition(depth);
0048       double EtaCentr = center_pos.eta();
0049       double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0050
0051       float dr = reco::deltaR(Eta, Phi, EtaCentr, PhiCentr);
0052       if (dr < drmin) {
0053         drmin = dr;
0054         crystalseed = crystal;
0055       }
0056     }
0057
0058     ieta = crystalseed.ieta();
0059     iphi = crystalseed.iphi();
0060
0061     // Get center cell position from shower depth
0062     auto cell = geom->getGeometry(crystalseed);
0063     const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());
0064
0065     thetatilt = cpyr->getThetaAxis();
0066     phitilt = cpyr->getPhiAxis();
0067
0068     GlobalPoint center_pos = cpyr->getPosition(depth);
0069
0070     double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0071     double PhiWidth = (M_PI / 180.);
0072     phicry = (TVector2::Phi_mpi_pi(Phi - PhiCentr)) / PhiWidth;
0073     //Some flips to take into account ECAL barrel symmetries:
0074     if (ieta < 0)
0075       phicry *= -1.;
0076
0077     double ThetaCentr = -center_pos.theta() + 0.5 * M_PI;
0078     double ThetaWidth = (M_PI / 180.) * std::cos(ThetaCentr);
0079     etacry = (Theta - ThetaCentr) / ThetaWidth;
0080     //flip to take into account ECAL barrel symmetries:
0081     if (ieta < 0)
0082       etacry *= -1.;
0083
0084     return;
0085   }
0086
0087   void localEcalClusterCoordsEE(const reco::CaloCluster &bclus,
0088                                 const CaloGeometry &caloGeometry,
0089                                 float &xcry,
0090                                 float &ycry,
0091                                 int &ix,
0092                                 int &iy,
0093                                 float &thetatilt,
0094                                 float &phitilt) {
0095     assert(bclus.hitsAndFractions().at(0).first.subdetId() == EcalEndcap);
0096
0097     const CaloSubdetectorGeometry *geom =
0098         caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);  //EcalBarrel = 1
0099
0100     const math::XYZPoint &position_ = bclus.position();
0101     //double Theta = -position_.theta()+0.5*M_PI;
0102     double Eta = position_.eta();
0103     double Phi = TVector2::Phi_mpi_pi(position_.phi());
0104     double X = position_.x();
0105     double Y = position_.y();
0106
0107     //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
0108     // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
0109     const float X0 = 0.89;
0110     float T0 = 1.2;
0111     //different T0 value if outside of preshower coverage
0112     if (std::abs(bclus.eta()) < 1.653)
0113       T0 = 3.1;
0114
0115     double depth = X0 * (T0 + log(bclus.energy()));
0116
0117     //find max energy crystal
0118     std::vector<std::pair<DetId, float> > crystals_vector = bclus.hitsAndFractions();
0119     float drmin = 999.;
0120     EEDetId crystalseed;
0121     //printf("starting loop over crystals, etot = %5f:\n",bclus.energy());
0122     for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
0123       EEDetId crystal(crystals_vector[icry].first);
0124
0125       auto cell = geom->getGeometry(crystal);
0126       const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());
0127       GlobalPoint center_pos = cpyr->getPosition(depth);
0128       double EtaCentr = center_pos.eta();
0129       double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
0130
0131       float dr = reco::deltaR(Eta, Phi, EtaCentr, PhiCentr);
0132       if (dr < drmin) {
0133         drmin = dr;
0134         crystalseed = crystal;
0135       }
0136     }
0137
0138     ix = crystalseed.ix();
0139     iy = crystalseed.iy();
0140
0141     // Get center cell position from shower depth
0142     auto cell = geom->getGeometry(crystalseed);
0143     const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());
0144
0145     thetatilt = cpyr->getThetaAxis();
0146     phitilt = cpyr->getPhiAxis();
0147
0148     GlobalPoint center_pos = cpyr->getPosition(depth);
0149
0150     double XCentr = center_pos.x();
0151     double XWidth = 2.59;
0152     xcry = (X - XCentr) / XWidth;
0153
0154     double YCentr = center_pos.y();
0155     double YWidth = 2.59;
0156     ycry = (Y - YCentr) / YWidth;
0157
0158     return;
0159   }
0160
0161 }  // namespace egammaTools