EgammaTools/src/ECALPositionCalculator.cc

0001
0002 #include "RecoEgamma/EgammaTools/interface/ECALPositionCalculator.h"
0003
0004 // Framework includes
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DataFormats/GeometryVector/interface/Pi.h"
0007 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0008 #include "MagneticField/Engine/interface/MagneticField.h"
0009
0010 static constexpr float R_ECAL = 136.5;
0011 static constexpr float Z_Endcap = 328.0;
0012 static constexpr float etaBarrelEndcap = 1.479;
0013
0014 namespace egammaTools {
0015
0016   double ecalPhi(const MagneticField &magField,
0017                  const math::XYZVector &momentum,
0018                  const math::XYZPoint &vertex,
0019                  const int charge) {
0020     // Get kinematic variables
0021
0022     float ptParticle = momentum.Rho();
0023     float etaParticle = momentum.eta();
0024     float phiParticle = momentum.phi();
0025     float vRho = vertex.Rho();
0026
0027     // Magnetic field
0028
0029     const float RBARM = 1.357;  // was 1.31 , updated on 16122003
0030     const float ZENDM = 3.186;  // was 3.15 , updated on 16122003
0031
0032     float rbend = RBARM - (vRho / 100.0);  //Assumed vRho in cm
0033     float bend = 0.3 * magField.inTesla(GlobalPoint(0., 0., 0.)).z() * rbend / 2.0;
0034     float phi = 0.0;
0035
0036     if (fabs(etaParticle) <= etaBarrelEndcap) {
0037       if (fabs(bend / ptParticle) <= 1.) {
0038         phi = phiParticle - asin(bend / ptParticle) * charge;
0039         if (phi > Geom::pi())
0040           phi = phi - Geom::twoPi();
0041         if (phi < -Geom::pi())
0042           phi = phi + Geom::twoPi();
0043       } else {
0044         edm::LogWarning("") << "[EcalPositionFromTrack::phiTransformation] Warning: "
0045                             << "Too low Pt, giving up";
0046         return phiParticle;
0047       }
0048
0049     }  // end if in the barrel
0050
0051     if (fabs(etaParticle) > etaBarrelEndcap) {
0052       float rHit = 0.0;
0053       rHit = ZENDM / sinh(fabs(etaParticle));
0054       if (fabs(((rHit - (vRho / 100.0)) / rbend) * bend / ptParticle) <= 1.0) {
0055         phi = phiParticle - asin(((rHit - (vRho / 100.0)) / rbend) * bend / ptParticle) * charge;
0056         if (phi > Geom::pi())
0057           phi = phi - Geom::twoPi();
0058         if (phi < -Geom::pi())
0059           phi = phi + Geom::twoPi();
0060       } else {
0061         edm::LogWarning("") << "[EcalPositionFromTrack::phiTransformation] Warning: "
0062                             << "Too low Pt, giving up";
0063         return phiParticle;
0064       }
0065
0066     }  // end if in the endcap
0067
0068     return phi;
0069   }
0070
0071   double ecalEta(const math::XYZVector &momentum, const math::XYZPoint &vertex) {
0072     // Get kinematic variables
0073
0074     float etaParticle = momentum.eta();
0075     float vZ = vertex.z();
0076     float vRho = vertex.Rho();
0077
0078     if (etaParticle != 0.0) {
0079       float theta = 0.0;
0080       float zEcal = (R_ECAL - vRho) * sinh(etaParticle) + vZ;
0081
0082       if (zEcal != 0.0)
0083         theta = atan(R_ECAL / zEcal);
0084       if (theta < 0.0)
0085         theta = theta + Geom::pi();
0086
0087       float ETA = -log(tan(0.5 * theta));
0088
0089       if (fabs(ETA) > etaBarrelEndcap) {
0090         float Zend = Z_Endcap;
0091         if (etaParticle < 0.0)
0092           Zend = -Zend;
0093         float Zlen = Zend - vZ;
0094         float RR = Zlen / sinh(etaParticle);
0095         theta = atan((RR + vRho) / Zend);
0096         if (theta < 0.0)
0097           theta = theta + Geom::pi();
0098         ETA = -log(tan(0.5 * theta));
0099       }
0100       return ETA;
0101
0102     } else {
0103       edm::LogWarning("") << "[EcalPositionFromTrack::etaTransformation] Warning: "
0104                           << "Eta equals to zero, not correcting";
0105       return etaParticle;
0106     }
0107   }
0108
0109 }  // namespace egammaTools