ECALPositionCalculator.cc

Go to the documentation of this file.

#include "RecoEgamma/EgammaTools/interface/ECALPositionCalculator.h"

// Framework includes
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/GeometryVector/interface/Pi.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
#include "MagneticField/Engine/interface/MagneticField.h"

static constexpr float R_ECAL = 136.5;
static constexpr float Z_Endcap = 328.0;
static constexpr float etaBarrelEndcap = 1.479;

namespace egammaTools {

  double ecalPhi(const MagneticField &magField,
                 const math::XYZVector &momentum,
                 const math::XYZPoint &vertex,
                 const int charge) {
    // Get kinematic variables

    float ptParticle = momentum.Rho();
    float etaParticle = momentum.eta();
    float phiParticle = momentum.phi();
    float vRho = vertex.Rho();

    // Magnetic field

    const float RBARM = 1.357;  // was 1.31 , updated on 16122003
    const float ZENDM = 3.186;  // was 3.15 , updated on 16122003

    float rbend = RBARM - (vRho / 100.0);  //Assumed vRho in cm
    float bend = 0.3 * magField.inTesla(GlobalPoint(0., 0., 0.)).z() * rbend / 2.0;
    float phi = 0.0;

    if (fabs(etaParticle) <= etaBarrelEndcap) {
      if (fabs(bend / ptParticle) <= 1.) {
        phi = phiParticle - asin(bend / ptParticle) * charge;
        if (phi > Geom::pi())
          phi = phi - Geom::twoPi();
        if (phi < -Geom::pi())
          phi = phi + Geom::twoPi();
      } else {
        edm::LogWarning("") << "[EcalPositionFromTrack::phiTransformation] Warning: "
                            << "Too low Pt, giving up";
        return phiParticle;
      }

    }  // end if in the barrel

    if (fabs(etaParticle) > etaBarrelEndcap) {
      float rHit = 0.0;
      rHit = ZENDM / sinh(fabs(etaParticle));
      if (fabs(((rHit - (vRho / 100.0)) / rbend) * bend / ptParticle) <= 1.0) {
        phi = phiParticle - asin(((rHit - (vRho / 100.0)) / rbend) * bend / ptParticle) * charge;
        if (phi > Geom::pi())
          phi = phi - Geom::twoPi();
        if (phi < -Geom::pi())
          phi = phi + Geom::twoPi();
      } else {
        edm::LogWarning("") << "[EcalPositionFromTrack::phiTransformation] Warning: "
                            << "Too low Pt, giving up";
        return phiParticle;
      }

    }  // end if in the endcap

    return phi;
  }

  double ecalEta(const math::XYZVector &momentum, const math::XYZPoint &vertex) {
    // Get kinematic variables

    float etaParticle = momentum.eta();
    float vZ = vertex.z();
    float vRho = vertex.Rho();

    if (etaParticle != 0.0) {
      float theta = 0.0;
      float zEcal = (R_ECAL - vRho) * sinh(etaParticle) + vZ;

      if (zEcal != 0.0)
        theta = atan(R_ECAL / zEcal);
      if (theta < 0.0)
        theta = theta + Geom::pi();

      float ETA = -log(tan(0.5 * theta));

      if (fabs(ETA) > etaBarrelEndcap) {
        float Zend = Z_Endcap;
        if (etaParticle < 0.0)
          Zend = -Zend;
        float Zlen = Zend - vZ;
        float RR = Zlen / sinh(etaParticle);
        theta = atan((RR + vRho) / Zend);
        if (theta < 0.0)
          theta = theta + Geom::pi();
        ETA = -log(tan(0.5 * theta));
      }
      return ETA;

    } else {
      edm::LogWarning("") << "[EcalPositionFromTrack::etaTransformation] Warning: "
                          << "Eta equals to zero, not correcting";
      return etaParticle;
    }
  }

}  // namespace egammaTools