EcalClusterLocal.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "DataFormats/CaloRecHit/interface/CaloCluster.h"
#include "TMath.h"
#include "TVector2.h"

namespace egammaTools {

  void localEcalClusterCoordsEB(const reco::CaloCluster &bclus,
                                const CaloGeometry &caloGeometry,
                                float &etacry,
                                float &phicry,
                                int &ieta,
                                int &iphi,
                                float &thetatilt,
                                float &phitilt) {
    assert(bclus.hitsAndFractions().at(0).first.subdetId() == EcalBarrel);

    const CaloSubdetectorGeometry *geom =
        caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalBarrel);  //EcalBarrel = 1

    const math::XYZPoint &position_ = bclus.position();
    double Theta = -position_.theta() + 0.5 * TMath::Pi();
    double Eta = position_.eta();
    double Phi = TVector2::Phi_mpi_pi(position_.phi());

    //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
    // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
    const float X0 = 0.89;
    const float T0 = 7.4;
    double depth = X0 * (T0 + log(bclus.energy()));

    //find max energy crystal
    std::vector<std::pair<DetId, float> > crystals_vector = bclus.hitsAndFractions();
    float drmin = 999.;
    EBDetId crystalseed;
    //printf("starting loop over crystals, etot = %5f:\n",bclus.energy());
    for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
      EBDetId crystal(crystals_vector[icry].first);

      auto cell = geom->getGeometry(crystal);
      const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());
      GlobalPoint center_pos = cpyr->getPosition(depth);
      double EtaCentr = center_pos.eta();
      double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());

      float dr = reco::deltaR(Eta, Phi, EtaCentr, PhiCentr);
      if (dr < drmin) {
        drmin = dr;
        crystalseed = crystal;
      }
    }

    ieta = crystalseed.ieta();
    iphi = crystalseed.iphi();

    // Get center cell position from shower depth
    auto cell = geom->getGeometry(crystalseed);
    const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());

    thetatilt = cpyr->getThetaAxis();
    phitilt = cpyr->getPhiAxis();

    GlobalPoint center_pos = cpyr->getPosition(depth);

    double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());
    double PhiWidth = (TMath::Pi() / 180.);
    phicry = (TVector2::Phi_mpi_pi(Phi - PhiCentr)) / PhiWidth;
    //Some flips to take into account ECAL barrel symmetries:
    if (ieta < 0)
      phicry *= -1.;

    double ThetaCentr = -center_pos.theta() + 0.5 * TMath::Pi();
    double ThetaWidth = (TMath::Pi() / 180.) * TMath::Cos(ThetaCentr);
    etacry = (Theta - ThetaCentr) / ThetaWidth;
    //flip to take into account ECAL barrel symmetries:
    if (ieta < 0)
      etacry *= -1.;

    return;
  }

  void localEcalClusterCoordsEE(const reco::CaloCluster &bclus,
                                const CaloGeometry &caloGeometry,
                                float &xcry,
                                float &ycry,
                                int &ix,
                                int &iy,
                                float &thetatilt,
                                float &phitilt) {
    assert(bclus.hitsAndFractions().at(0).first.subdetId() == EcalEndcap);

    const CaloSubdetectorGeometry *geom =
        caloGeometry.getSubdetectorGeometry(DetId::Ecal, EcalEndcap);  //EcalBarrel = 1

    const math::XYZPoint &position_ = bclus.position();
    //double Theta = -position_.theta()+0.5*TMath::Pi();
    double Eta = position_.eta();
    double Phi = TVector2::Phi_mpi_pi(position_.phi());
    double X = position_.x();
    double Y = position_.y();

    //Calculate expected depth of the maximum shower from energy (like in PositionCalc::Calculate_Location()):
    // The parameters X0 and T0 are hardcoded here because these values were used to calculate the corrections:
    const float X0 = 0.89;
    float T0 = 1.2;
    //different T0 value if outside of preshower coverage
    if (TMath::Abs(bclus.eta()) < 1.653)
      T0 = 3.1;

    double depth = X0 * (T0 + log(bclus.energy()));

    //find max energy crystal
    std::vector<std::pair<DetId, float> > crystals_vector = bclus.hitsAndFractions();
    float drmin = 999.;
    EEDetId crystalseed;
    //printf("starting loop over crystals, etot = %5f:\n",bclus.energy());
    for (unsigned int icry = 0; icry != crystals_vector.size(); ++icry) {
      EEDetId crystal(crystals_vector[icry].first);

      auto cell = geom->getGeometry(crystal);
      const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());
      GlobalPoint center_pos = cpyr->getPosition(depth);
      double EtaCentr = center_pos.eta();
      double PhiCentr = TVector2::Phi_mpi_pi(center_pos.phi());

      float dr = reco::deltaR(Eta, Phi, EtaCentr, PhiCentr);
      if (dr < drmin) {
        drmin = dr;
        crystalseed = crystal;
      }
    }

    ix = crystalseed.ix();
    iy = crystalseed.iy();

    // Get center cell position from shower depth
    auto cell = geom->getGeometry(crystalseed);
    const TruncatedPyramid *cpyr = dynamic_cast<const TruncatedPyramid *>(cell.get());

    thetatilt = cpyr->getThetaAxis();
    phitilt = cpyr->getPhiAxis();

    GlobalPoint center_pos = cpyr->getPosition(depth);

    double XCentr = center_pos.x();
    double XWidth = 2.59;
    xcry = (X - XCentr) / XWidth;

    double YCentr = center_pos.y();
    double YWidth = 2.59;
    ycry = (Y - YCentr) / YWidth;

    return;
  }

}  // namespace egammaTools