FastCaloSimCaloExtrapolation.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
 
/* Athena includes */
#include "AthenaBaseComps/AthAlgTool.h"
#include "GaudiKernel/ToolHandle.h"
#include "GaudiKernel/IPartPropSvc.h"
 
#include "CxxUtils/inline_hints.h"
/* Header include */
#include "FastCaloSimCaloExtrapolation.h"
 
/* ISF includes */
#include "ISF_FastCaloSimEvent/FastCaloSim_CaloCell_ID.h"
#include "ISF_FastCaloSimParametrization/IFastCaloSimCaloTransportation.h"
#include "ISF_FastCaloSimEvent/TFCSTruthState.h"
 
/* Geometry primitives */
#include "GeoPrimitives/GeoPrimitivesHelpers.h"
#include "GeoPrimitives/CLHEPtoEigenConverter.h"
 
/* Particle data */
#include "HepPDT/ParticleDataTable.hh"
 
/* Tracking includes */
#include "TrkGeometry/TrackingGeometry.h"
 
/* G4FieldTrack used to store transportation steps */
#include "G4FieldTrack.hh"
 
 
/* Preprocessor macro to use
   -- DEBUG   if CONDITION is True
   -- WARNING if CONDITION is False
 */
#define ATH_MSG_COND(MSG, CONDITION) \
do {                                 \
    if (CONDITION) {                 \
      ATH_MSG_DEBUG(MSG);            \
    } else {                         \
      ATH_MSG_WARNING(MSG);          \
    }                                \
} while (0)
 
 
FastCaloSimCaloExtrapolation::FastCaloSimCaloExtrapolation(const std::string& t, const std::string& n, const IInterface* p)
  : base_class(t,n,p)
{
}
 
StatusCode FastCaloSimCaloExtrapolation::initialize()
{
  ATH_MSG_INFO( "Initializing FastCaloSimCaloExtrapolation" );
 
  // Retrieve the fast calo sim geometry helper
  ATH_CHECK(m_CaloGeometryHelper.retrieve());
  // Retrieve the tool to transport particles through calorimeter with ATLAS tracking tools
  ATH_CHECK(m_CaloTransportation.retrieve());
 
  return StatusCode::SUCCESS;
 
 
}
 
StatusCode FastCaloSimCaloExtrapolation::finalize(){
  ATH_MSG_INFO( "Finalizing FastCaloSimCaloExtrapolation" );
  return StatusCode::SUCCESS;
}
 
 
void FastCaloSimCaloExtrapolation::extrapolate(TFCSExtrapolationState& result, const TFCSTruthState* truth, const std::vector<G4FieldTrack>& caloSteps) const{ 
  
  ATH_MSG_DEBUG("[extrapolate] Initializing extrapolation to ID-Calo boundary");
  extrapolateToID(result, caloSteps, truth);
 
  ATH_MSG_DEBUG("[extrapolate] Initializing extrapolation to calorimeter layers");
  extrapolateToLayers(result, caloSteps, truth);
 
  ATH_MSG_DEBUG("[extrapolate] Extrapolation done");
 
}
 
void FastCaloSimCaloExtrapolation::extrapolate(TFCSExtrapolationState& result, const TFCSTruthState* truth) const{ 
  
  ATH_MSG_DEBUG("[extrapolate] Initializing transport of track through calorimeter system with ATLAS tracking tools.");
  std::vector<G4FieldTrack> caloSteps = m_CaloTransportation -> transport(truth, false);
  ATH_MSG_DEBUG("[extrapolate] Finalized transport of track through calorimeter system with ATLAS tracking tools.");
  
  extrapolate(result, truth, caloSteps);
}
 
 
void FastCaloSimCaloExtrapolation::extrapolateToID(TFCSExtrapolationState& result, const std::vector<G4FieldTrack>& caloSteps, const TFCSTruthState* truth) const{
 
  ATH_MSG_DEBUG("Start extrapolateToID()");
 
  //pT threshold of truth particles over which extrapolation failures will be printed as warnings
  const float transverseMomWarningLimit = 500;
 
  //initialize values
  result.set_IDCaloBoundary_eta(-999.);
  result.set_IDCaloBoundary_phi(-999.);
  result.set_IDCaloBoundary_r(0);
  result.set_IDCaloBoundary_z(0);
  result.set_IDCaloBoundary_AngleEta(-999.);
  result.set_IDCaloBoundary_Angle3D(-999.);
  
  //magnitude of extrapolated position
  double extPosDist = -1;
 
  for (unsigned int surfID = 0; surfID<3; surfID++){
 
    double R = m_CaloBoundaryR.value().at(surfID);
    double Z = m_CaloBoundaryZ.value().at(surfID);
 
    ATH_MSG_DEBUG("[ExtrapolateToID] Extrapolating to ID-Calo boundary with ID="<<surfID<<" R="<<R<<" Z="<<Z);
 
    //extrapolated position and momentum direction at IDCaloBoundary
    Amg::Vector3D extPos, momDir;
    
    //main extrapolation call
    if(!extrapolateToCylinder(caloSteps, R, Z, extPos, momDir)) continue;
    
    double tolerance = 0.001;
    
    //test if z inside previous cylinder within some tolerance
    ATH_MSG_DEBUG("[ExtrapolateToID] Testing condition 1: hit z="<< extPos[Amg::z]);
    if(surfID > 0 && std::abs(extPos[Amg::z]) < m_CaloBoundaryZ[surfID-1] - tolerance) continue;
    ATH_MSG_DEBUG("[ExtrapolateToID] Passed condition 1.");
 
    //test if r inside next cylinder within some tolerance
    ATH_MSG_DEBUG("[ExtrapolateToID] Testing condition 2: hit r="<< extPos.perp());
    if(surfID < m_CaloBoundaryR.size()-1 && extPos.perp() < m_CaloBoundaryR[surfID + 1] - tolerance) continue;
    ATH_MSG_DEBUG("[ExtrapolateToID] Passed condition 2.");
 
    ATH_MSG_DEBUG("[ExtrapolateToID] Testing condition 3: hit magnitude="<< extPos.mag());
    if(extPosDist >= 0 && extPos.mag() > extPosDist) continue;
    ATH_MSG_DEBUG("[ExtrapolateToID] Passed condition 3.");
 
    extPosDist = extPos.mag();
 
    result.set_IDCaloBoundary_eta(extPos.eta());
    result.set_IDCaloBoundary_phi(extPos.phi());
    result.set_IDCaloBoundary_r(extPos.perp());
    result.set_IDCaloBoundary_z(extPos[Amg::z]);
 
    ATH_MSG_DEBUG("[ExtrapolateToID] Setting IDCaloBoundary to eta="<<extPos.eta()<<" phi="<<extPos.phi()<< " r="<<extPos.perp()<<" z="<<extPos.z());
 
    //compute angle between extrapolated position vector and momentum at IDCaloBoundary
    //can be used to correct shower shapes for particles which do not originate from {0,0,0}
    double Angle3D  = Amg::angle(extPos, momDir);
    double AngleEta = extPos.theta() - momDir.theta();
    result.set_IDCaloBoundary_AngleEta(AngleEta);
    result.set_IDCaloBoundary_Angle3D(Angle3D);
 
  } //end of loop over surfaces
  
  if(result.IDCaloBoundary_eta() == -999) ATH_MSG_COND("[ExtrapolateToID] Failed extrapolation to ID-Calo boundary. \n[ExtrapolateToID] Particle with truth vertex at (" << truth->vertex().X() <<","<<truth->vertex().Y()<<","<<truth->vertex().Z()<<")"<<" with"<<" PdgId="<<truth->pdgid()<<" pT="<<truth->Pt()<<" eta="<<truth->Eta()<<" phi="<<truth->Phi()<<" E="<<truth->E()<<" Ekin_off="<<truth->Ekin_off(), truth->Pt() < transverseMomWarningLimit);
  
  ATH_MSG_DEBUG("[ExtrapolateToID] End extrapolateToID()");
 
} 
 
 
void FastCaloSimCaloExtrapolation::extrapolateToLayers(TFCSExtrapolationState& result, const std::vector<G4FieldTrack>& caloSteps, const TFCSTruthState* truth) const
{
  ATH_MSG_DEBUG("[extrapolateToLayers] Start extrapolate");
 
  //pT threshold of truth particles over which extrapolation failures will be printed as warnings
  const float transverseMomWarningLimit = 500;
 
  // Start calo extrapolation
 
  //only continue if inside the calo
  if(std::abs(result.IDCaloBoundary_eta()) < 6){
      //now try to extrapolate to all calo layers that contain energy
      for(int sample=CaloCell_ID_FCS::FirstSample; sample<CaloCell_ID_FCS::MaxSample; ++sample){   
          for(int subpos=SUBPOS_MID; subpos<=SUBPOS_EXT; ++subpos){ 
            
            float cylR, cylZ;
            if(isCaloBarrel(sample)){
              cylR = std::abs(rpos(sample, result.IDCaloBoundary_eta(), subpos));
              //EMB0 - EMB3 use z position of EME1 front end surface for extrapolation
              //else extrapolate to cylinder with symmetrized maximum Z bounds
              //set eta to a dummy value of 1000 and -1000 to force detector side
              if(sample < 4) cylZ = result.IDCaloBoundary_eta() > 0 ? std::abs(zpos(5, 1000, 1)) : std::abs(zpos(5, -1000, 1));
              else cylZ = 0.5*(std::abs(zpos(sample, 1000, subpos)) + std::abs(zpos(sample, -1000, subpos)));
            }
            else{
              //if we are not at barrel surface, extrapolate to cylinder with maximum R to reduce extrapolation length
              cylZ = std::abs(zpos(sample, result.IDCaloBoundary_eta(), subpos));
              //calculate radius of cylinder we will extrapolate to
              double mineta, maxeta, eta;
              minmaxeta(sample, result.IDCaloBoundary_eta(), mineta, maxeta);
              //get eta where we will look up the layer radius
              eta = result.IDCaloBoundary_eta() > 0 ? mineta : maxeta;
              //calculate azimuthal angle from pseudorapidity
              double theta = 2*std::atan(std::exp(-eta));
              //calculate maximum R of last cell of layer from z and theta
              cylR = std::abs(cylZ*std::sqrt((1/(std::cos(theta)*std::cos(theta))) - 1));
            }
 
            Amg::Vector3D extPos, momDir;
            if(extrapolateToCylinder(caloSteps, cylR, cylZ, extPos, momDir)){
 
              //scale the extrapolation to fit the radius of the cylinder in the case of barrel and scale extrapolation to fit z component in case of endcap layer
              //scale is only non-unitary in case we extrapolate to the endcaps of the cylinder for barrel and in case we extrapolate to cover for endcaps
              //this will keep phi, eta intact and only scale r and z to fit a sensible position on the cylinder
              double scale = 1;
              if (isCaloBarrel(sample) && std::abs(extPos.perp()) > 1e-6)     scale = cylR / extPos.perp();
              else if (!isCaloBarrel(sample) && std::abs(extPos.z()) > 1e-6)  scale = cylZ / std::abs(extPos.z());
              //scale extrapolated position accordingly
              extPos = scale * extPos; 
 
              result.set_OK(sample, subpos,  true);
              result.set_phi(sample, subpos, extPos.phi());
              result.set_z  (sample, subpos, extPos.z());
              result.set_eta(sample, subpos, extPos.eta());
              result.set_r  (sample, subpos, extPos.perp());
            }
            else{
              ATH_MSG_COND(" [extrapolateToLayers] Extrapolation to cylinder failed. Sample="<<sample<<" subpos="<<subpos<<" eta="<<result.IDCaloBoundary_eta()<<" phi="<<result.IDCaloBoundary_phi()<<" r="<<result.IDCaloBoundary_r()<<" z="<<result.IDCaloBoundary_z(), truth->Pt() < transverseMomWarningLimit);
            }   
          } //for pos
        } //for sample
    } //inside calo
  
  else ATH_MSG_COND("[extrapolateToLayers] Ups. Not inside calo. result.IDCaloBoundary_eta()="<<result.IDCaloBoundary_eta()<< "\n[extrapolateToLayers] Particle with truth vertex at (" << truth->vertex().X() <<","<<truth->vertex().Y()<<","<<truth->vertex().Z()<<")"<<" with"<<" PdgId="<<truth->pdgid()<<" pT="<<truth->Pt()<<" eta="<<truth->Eta()<<" phi="<<truth->Phi()<<" E="<<truth->E()<<" Ekin_off="<<truth->Ekin_off(), truth->Pt() < transverseMomWarningLimit);
  
 
  ATH_MSG_DEBUG("[extrapolateToLayers] End extrapolateToLayers()");
}
 
bool FastCaloSimCaloExtrapolation::extrapolateToCylinder(const std::vector<G4FieldTrack>& caloSteps, float cylR, float cylZ, Amg::Vector3D& extPos, Amg::Vector3D& momDir) const{
 
  if(caloSteps.size() == 1){
    Amg::Vector3D hitPos = Amg::Hep3VectorToEigen(caloSteps.at(0).GetPosition());
    ATH_MSG_DEBUG("[extrapolateWithPCA(R="<<cylR<<",Z="<<cylZ<<")] Extrapolating single hit position to surface.");
    extPos = projectOnCylinder(cylR, cylZ, hitPos);
    momDir = Amg::Hep3VectorToEigen(caloSteps.at(0).GetMomentum());
    return true;
  }
 
  //if we do not find any good intersections, extrapolate to closest point on surface 
  bool foundHit = extrapolateWithIntersection(caloSteps, cylR, cylZ, extPos, momDir) ? true : extrapolateWithPCA(caloSteps, cylR, cylZ, extPos, momDir);
 
  if(foundHit){
    ATH_MSG_DEBUG("[extrapolateToCylinder(R="<<cylR<<",Z="<<cylZ<<")::END] Extrapolated to cylinder with R="<<cylR<<" and Z="<<cylZ<<" at ("<< extPos[Amg::x]<<","<<extPos[Amg::y]<<","<<extPos[Amg::z]<<")");
  }
  else{ 
    //this is not expected to ever happen
    ATH_MSG_DEBUG("(R="<<cylR<<", Z="<<cylZ<<"::END) Extrapolation to cylinder surface failed!"); 
  }
 
 
  return foundHit;
 
}
 
 
bool FastCaloSimCaloExtrapolation::extrapolateWithIntersection(const std::vector<G4FieldTrack>& caloSteps, float cylR, float cylZ, Amg::Vector3D& extPos, Amg::Vector3D& momDir) const{
 
  ATH_MSG_DEBUG("[extrapolateWithIntersection(R="<<cylR<<",Z="<<cylZ<<")] Checking for cylinder intersections of line segments.");
  
  //counter for number of computed extrapolations, does not count cases of rejected extrapolations due to close by hit positions
  unsigned int nExtrapolations = 0;
  for (size_t hitID = 1; hitID < caloSteps.size(); hitID++){   
    //initialize intersection result variables
    //get current and consecutive hit position and build hitLine
    Amg::Vector3D hitPos1 = Amg::Hep3VectorToEigen(caloSteps.at(hitID-1).GetPosition());
    Amg::Vector3D hitPos2 = Amg::Hep3VectorToEigen(caloSteps.at(hitID).GetPosition());
    Amg::Vector3D hitDir  = hitPos2 - hitPos1; 
 
    ATH_MSG_DEBUG("[extrapolateWithIntersection(R="<<cylR<<",Z="<<cylZ<<")] Considering line segment between ("<<hitPos1[Amg::x]<<","<<hitPos1[Amg::y]<<","<<hitPos1[Amg::z]<<") and ("
                                                                                                         <<hitPos2[Amg::x]<<","<<hitPos2[Amg::y]<<","<<hitPos2[Amg::z]<<")");
    //get position of the hit positions on the cylinder
    HITPOSITION cylPosHit1 = whereOnCylinder(cylR, cylZ, hitPos1);    
    HITPOSITION cylPosHit2 = whereOnCylinder(cylR, cylZ, hitPos2);    
    
    //check if one of the hit positions already lays on the cylinder surface    
    if(cylPosHit1 == ON || cylPosHit2 == ON){   
      extPos = cylPosHit1 == ON ? hitPos1 : hitPos2;    
      momDir = cylPosHit1 == ON ? Amg::Hep3VectorToEigen(caloSteps.at(hitID-1).GetMomentum()) : Amg::Hep3VectorToEigen(caloSteps.at(hitID).GetMomentum());     
      ATH_MSG_DEBUG("[extrapolateWithIntersection(R="<<cylR<<",Z="<<cylZ<<")] Hit position already on cylinder surface.");    
      return true;    
    }
 
    //do not try to extrapolate with intersections if the hit position are very close together
    if(hitDir.norm() < 0.01) continue;
    
    //get intersections through cylinder
    CylinderIntersections intersections = getCylinderIntersections(cylR, cylZ, hitPos1, hitPos2);
    nExtrapolations++;
    
    Amg::Vector3D selectedIntersection(0, 0, 0); 
 
    //select the best intersection
    if(intersections.number == 1)      selectedIntersection = intersections.first;
    else if(intersections.number > 1)  selectedIntersection = whichIntersection(cylR, cylZ, hitPos1, hitPos2, intersections.first, intersections.second) == 0 ? 
                                                              intersections.first : intersections.second; 
                                     
    if(intersections.number > 0){
      
      bool isForwardExtrapolation = (selectedIntersection[Amg::x] - hitPos1[Amg::x]) /  (hitPos2[Amg::x] -  hitPos1[Amg::x]) >= 0;
      bool travelThroughSurface = doesTravelThroughSurface(cylR, cylZ, hitPos1, hitPos2);
      
      //do not allow for backward extrapolation except in the case of first two (distinguishable) hit positions outside cylinder 
      //and in the case we detect a travel though the surface
      if(nExtrapolations > 1 && !isForwardExtrapolation && !travelThroughSurface) continue;
      
      //check if the intersection between infinite line and cylinder lays on segment spanned by hit positions
      bool intersectionOnSegment = isOnSegment(selectedIntersection, hitPos1, hitPos2);
      //check if both hit positions lay outside of the cylinder
      bool hitPosOutside = cylPosHit1 == OUTSIDE && cylPosHit2 == OUTSIDE;
      
      //we found our extrapolated hit position in case that either
      //we detect that the line segment crosses the surface of the cylinder
      //the intersection between the infinite lines and the cylinder lays on the line segment
      //both hit positions are outside of the cylinder and there is a backwards extrapolation for the first two hit positions
      //if this is not the case for any of the hit position pairs we will use the last two hit position for the linear extrapolation
      //if these do not have any intersection, then we will pass back to extrapolateWithPCA
      if(travelThroughSurface || intersectionOnSegment || (hitPosOutside && !isForwardExtrapolation && nExtrapolations == 1) || caloSteps.size()-1 == hitID){
        //take momentum direction of hit position closest to cylinder surface
        //alternatively one could also take the extrapolated direction normDir = hitPos2 - hitPos1
        double distHitPos1 = (hitPos1 - projectOnCylinder(cylR, cylZ, hitPos1)).norm(); 
        double distHitPos2 = (hitPos2 - projectOnCylinder(cylR, cylZ, hitPos2)).norm(); 
        momDir = distHitPos1 < distHitPos2 ? Amg::Hep3VectorToEigen(caloSteps.at(hitID-1).GetMomentum()) : Amg::Hep3VectorToEigen(caloSteps.at(hitID).GetMomentum());
        extPos = selectedIntersection; 
        return true;
      } 
      ATH_MSG_DEBUG("[extrapolateWithIntersection(R="<<cylR<<",Z="<<cylZ<<")] Extrapolated position at ("<<selectedIntersection[Amg::x]<<","<<selectedIntersection[Amg::y]<<","<<selectedIntersection[Amg::z]<<")");
    }
  } //end of loop over hit positions
 
  return false;
}
 
 
bool FastCaloSimCaloExtrapolation::extrapolateWithPCA(const std::vector<G4FieldTrack>& caloSteps, float cylR, float cylZ, Amg::Vector3D& extPos, Amg::Vector3D& momDir) const{
 
  bool foundHit = false;
  ATH_MSG_DEBUG("[extrapolateWithPCA(R="<<cylR<<",Z="<<cylZ<<")] No forward intersections with cylinder surface. Extrapolating to closest point on surface.");
 
  //here we also need to consider distances from line segments to the cylinder 
  double minDistToSurface = 100000;
  for (size_t hitID = 1; hitID < caloSteps.size(); hitID++){   
   
    Amg::Vector3D hitPos1 = Amg::Hep3VectorToEigen(caloSteps.at(hitID-1).GetPosition());
    Amg::Vector3D hitPos2 = Amg::Hep3VectorToEigen(caloSteps.at(hitID).GetPosition());
 
    ATH_MSG_DEBUG("[extrapolateWithPCA(R="<<cylR<<",Z="<<cylZ<<")] Considering line segment between ("<<hitPos1[Amg::x]<<","<<hitPos1[Amg::y]<<","<<hitPos1[Amg::z]<<") and ("<<hitPos2[Amg::x]<<","<<hitPos2[Amg::y]<<","<<hitPos2[Amg::z]<<")");
 
    Amg::Vector3D PCA;
    //find the point of closest approach (PCA) to the cylinder on the line segment 
    findPCA(cylR, cylZ, hitPos1, hitPos2, PCA);
    //compute distance between PCA and cylinder
    Amg::Vector3D cylinderSurfacePCA = projectOnCylinder(cylR, cylZ, PCA); 
    double tmpMinDistToSurface = (PCA - cylinderSurfacePCA).norm();
    
    ATH_MSG_DEBUG("[extrapolateWithPCA(R="<<cylR<<",Z="<<cylZ<<")] Extrapolated line segment to ("<<cylinderSurfacePCA[Amg::x]<<","<<cylinderSurfacePCA[Amg::y]<<","<<cylinderSurfacePCA[Amg::z]<<") with distance "<<tmpMinDistToSurface);
    
    if(tmpMinDistToSurface < minDistToSurface){
      foundHit = true;
      extPos = cylinderSurfacePCA;
      //take momentum direction of hit position closest to cylinder surface
      //alternatively one could also take the extrapolated direction normDir = hitPos2 - hitPos1
      double distHitPos1 = (hitPos1 - projectOnCylinder(cylR, cylZ, hitPos1)).norm(); 
      double distHitPos2 = (hitPos2 - projectOnCylinder(cylR, cylZ, hitPos2)).norm(); 
      momDir = distHitPos1 < distHitPos2 ? Amg::Hep3VectorToEigen(caloSteps.at(hitID-1).GetMomentum()) : Amg::Hep3VectorToEigen(caloSteps.at(hitID).GetMomentum());
 
      minDistToSurface = tmpMinDistToSurface;
    }
  } //end over loop of hit postions
 
  return foundHit;
} 
 
 
void FastCaloSimCaloExtrapolation::findPCA(float cylR, float cylZ, Amg::Vector3D& hitPos1, Amg::Vector3D& hitPos2, Amg::Vector3D& PCA) const{
  //in the following we will try to find the closest point-of-approach (PCA) to the cylinder on the line segment
  //hit direction
  Amg::Vector3D hitDir = hitPos2 - hitPos1;
 
  //project both hit positions onto the cylinder
  Amg::Vector3D projCylinderHitPos1 = projectOnCylinder(cylR, cylZ, hitPos1);
  Amg::Vector3D projCylinderHitPos2 = projectOnCylinder(cylR, cylZ, hitPos2);
  //direction of line spanned by the two projected points on the cylinder surface
  Amg::Vector3D cylinderProjDir = projCylinderHitPos2 - projCylinderHitPos1;
 
  //CASE A: projections on the cylinder are close enough, take one of the hit positions as PCA 
  if(cylinderProjDir.norm() < 0.0001) {PCA = hitPos1; return;};   
 
  //CASE B: we are outside the Z bounds of the cylinder
  if((hitPos1[Amg::z] > cylZ || hitPos1[Amg::z] < -cylZ) || (hitPos2[Amg::z] > cylZ || hitPos2[Amg::z] < -cylZ)){
           
    //calculate PCA to point on endcap
    Amg::Vector3D cylZEndcap(0, 0, cylZ);
    bool isParallelToEndcap = std::abs(hitPos1[Amg::z] - hitPos2[Amg::z]) < 0.00001;
 
    //Check if parallel to endcap plane
    if(isParallelToEndcap){
 
      //if both inside there are infinite solutions take one in the middle
      Amg::Vector3D intersectA, intersectB;
      intersectA.setZero();
      intersectB.setZero();
      int nIntersections = circleLineIntersection2D(cylR, hitPos1, hitPos2, intersectA, intersectB);
      
      if(nIntersections == 2){
        
        bool IntAOnSegment = isOnSegment(intersectA, hitPos1, hitPos2);
        bool IntBOnSegment = isOnSegment(intersectB, hitPos1, hitPos2);
 
        if(IntAOnSegment && IntBOnSegment) PCA = intersectA + 0.5*(intersectB-intersectA);
        else if(IntAOnSegment) PCA = hitPos1.perp() <= cylR ?  intersectA + 0.5*(hitPos1 - intersectA) : intersectA + 0.5*(hitPos2 - intersectA);
        else if(IntBOnSegment) PCA = hitPos1.perp() <= cylR ?  intersectB + 0.5*(hitPos1 - intersectB) : intersectB + 0.5*(hitPos2 - intersectB);
        //intersections are not on line segment, i.e. line segment is within extended cylinder
        else PCA = hitPos1 + 0.5*hitDir;
      
      }
      else if(!intersectA.isZero() || !intersectB.isZero()){
        //this can only happen if the extended line is tangetial to the cylinder
        //if intersection lays on segment PCA will be intersection, if not it will be the corresponding end points
        Amg::Vector3D intersect = intersectA.isZero() ? intersectB : intersectA;
        Amg::Vector3D hitPos = (hitPos1 - intersect).norm() < (hitPos2 - intersect).norm() ? hitPos1 : hitPos2;
        bool IntOnSegment = isOnSegment(intersectA, hitPos1, hitPos2);
        PCA = IntOnSegment ? intersect : hitPos;
      
      }
      else{
        //line segment is outside extended cylinder
        //PCA corresponds to closest distance to center {0, 0, cylZ}
        Amg::Vector3D infLinePCA = hitPos1 + ((cylZEndcap-hitPos1).dot(hitDir)/hitDir.dot(hitDir))*(hitDir);
 
        if(isOnSegment(infLinePCA, hitPos1, hitPos2)) PCA = infLinePCA;
        else PCA = (hitPos1 - infLinePCA).norm() < (hitPos2 - infLinePCA).norm() ? hitPos1 : hitPos2;    
        
      }
    }
 
    else{
 
      //figure out all other cases iteratively beginning with BoundA and BoundB
      Amg::Vector3D BoundA, BoundB;
      //this is point on line closest to {0, 0, cylZ}, always on segment
      double t =  ((cylZEndcap-hitPos1).dot(hitDir)/hitDir.dot(hitDir));
      BoundA = t <= 0 ? hitPos1 : (t >= 1 ? hitPos2 : hitPos1 + t*hitDir);
      
      //calculate intersection point of line segment and endcap plane and project intersection onto cylinder
      //check if t is between 0 and 1, if not, take hitpos as starting bound
      t = (cylZ-hitPos1[Amg::z]) / hitDir[Amg::z];
      BoundB = t <= 0 ? hitPos1 : (t >= 1 ? hitPos2 : hitPos1 + t*hitDir);
      //looks for the PCA iteratively in cases there is no easy analytical solution
      getIterativePCA(cylR, cylZ, BoundA, BoundB, PCA);
      
    }
 
    return;
  }
 
  //CASE C: we are inside the Z bounds of the cylinder
  //construct Z axis as straight line surface
  Trk::StraightLineSurface line(Amg::Transform3D(Trk::s_idTransform), 0, cylZ);
  //compute point of closest approach to z axis
  //this is analogous to finding the PCA of two 3D lines
  Trk::Intersection PCACylBounds = line.straightLineIntersection(hitPos1, hitDir.unit(), false, true);
 
  double distSurfHit1 = (projCylinderHitPos1 - hitPos1).norm();
  double distSurfHit2 = (projCylinderHitPos2 - hitPos2).norm();
  
  //take PCA on line in case it lays on segment, otherwise take closest hit position to surface
  PCA = isOnSegment(PCACylBounds.position, hitPos1, hitPos2) ? PCACylBounds.position : (distSurfHit1 < distSurfHit2 ? hitPos1 : hitPos2);
  
}
 
 
#if defined(FLATTEN)
// We compile this package with optimization, even in debug builds; otherwise,
// the heavy use of Eigen makes it too slow.  However, from here we may call
// to out-of-line Eigen code that is linked from other DSOs; in that case,
// it would not be optimized.  Avoid this by forcing all Eigen code
// to be inlined here if possible.
ATH_FLATTEN
#endif
void FastCaloSimCaloExtrapolation::getIterativePCA(float cylR, float cylZ, Amg::Vector3D& BoundA, Amg::Vector3D& BoundB, Amg::Vector3D& PCA) const{
 
    ATH_MSG_DEBUG("[getIterativePCA] Finding PCA iteratively.");
    
    Amg::Vector3D boundDir = BoundB - BoundA;
    double distBounds = boundDir.norm();
 
    //this sets the precision of the iterative finding procedure
    const double stepSize = 0.01;
 
    //if bounds are close enough together, there is nothing to do
    //should make sure that nHalfDivisions >= 2
    if (distBounds <= 4*stepSize){PCA = BoundA + 0.5*(BoundB - BoundA); return;} 
 
    Amg::Vector3D tmpBoundA, tmpBoundB, tmpOnCylinderBoundA, tmpOnCylinderBoundB;
    Amg::Vector3D resBoundA, resBoundB, resOnCylinderBoundA, resOnCylinderBoundB;
 
    
    //initial positions on cylinder and distance to line segment
    Amg::Vector3D OnCylinderBoundA = projectOnCylinder(cylR, cylZ, BoundA);
    Amg::Vector3D OnCylinderBoundB = projectOnCylinder(cylR, cylZ, BoundB);
 
    double minDistA = (BoundA - OnCylinderBoundA).norm();
    double minDistB = (BoundB - OnCylinderBoundB).norm();
 
    //initialize result bounds with closest input bounds as fall back option
    if(minDistA < minDistB){
      resBoundA = BoundA;
      resBoundB = BoundA;
    }
    else{
      resBoundA = BoundB;
      resBoundB = BoundB;
    }
    double tmpMinDistA, tmpMinDistB;
    unsigned int nHalfDivisions = (distBounds/stepSize)/2;
    for(unsigned int step = 0; step < nHalfDivisions; step++){
 
      //temporary bounds on line segment
      tmpBoundA = BoundA + (step+1)*stepSize*(boundDir/distBounds);
      tmpBoundB = BoundB - (step+1)*stepSize*(boundDir/distBounds);
 
      //temporary projected bounds on cylinder
      tmpOnCylinderBoundA = projectOnCylinder(cylR, cylZ, tmpBoundA);
      tmpOnCylinderBoundB = projectOnCylinder(cylR, cylZ, tmpBoundB);
      
      //temporary minimum distance between bound on segment and bound on cylinder
      tmpMinDistA = (tmpBoundA - tmpOnCylinderBoundA).norm();
      tmpMinDistB = (tmpBoundB - tmpOnCylinderBoundB).norm();
 
      if(minDistA >= tmpMinDistA){
        minDistA = tmpMinDistA;
      }
      else{
        double t = (step*stepSize)/distBounds;
        resBoundA = BoundA + t*boundDir;
        resBoundB = tmpBoundA;
        break;
      }
 
      if(minDistB >= tmpMinDistB){
        minDistB = tmpMinDistB;
      }
      else{
        double t = (step*stepSize)/distBounds;
        resBoundB = BoundB - t*boundDir;
        resBoundA = tmpBoundB;
        break;
      }
    }
 
    //return middle of best bounds
    PCA = resBoundA + 0.5*(resBoundB - resBoundA);
 
}
 
 
int FastCaloSimCaloExtrapolation::circleLineIntersection2D(float circR, Amg::Vector3D& pointA, Amg::Vector3D& pointB, Amg::Vector3D& intersectA, Amg::Vector3D& intersectB) const{
    //find intersections intA and intB with line spanned by pointA and pointB
    //returns number of intersections
    //assumes circle lays in xy plane
 
    double dx, dy, A, B, C, det, t;
 
    dx = pointB[Amg::x] - pointA[Amg::x];
    dy = pointB[Amg::y] - pointA[Amg::y];
 
    A = dx * dx + dy * dy;
    B = 2 * (dx * pointA[Amg::x] + dy * pointA[Amg::y]);
    C = pointA[Amg::x] * pointA[Amg::x] + pointA[Amg::y] * pointA[Amg::y] - circR * circR;
 
    det = B * B - 4 * A * C;
    
    if (A <= 0.0000001 || det < 0){   
        ATH_MSG_DEBUG("[circleLineIntersection2D] No intersections.");
        return 0;
    }
    else if (std::abs(det) < 0.00001){
        //one solution, tangential case.
        t = -B / (2 * A);
        intersectA = {pointA[Amg::x] + t * dx, pointA[Amg::y] + t * dy, pointA[Amg::z]};
        ATH_MSG_DEBUG("[circleLineIntersection2D] One intersection at ("<<intersectA[Amg::x]<<","<<intersectA[Amg::y]<<","<<intersectA[Amg::z]<<").");
        return 1;
    }
    else{
        // two solutions
        t = (-B + std::sqrt(det)) / (2 * A);
        intersectA = {pointA[Amg::x] + t * dx, pointA[Amg::y] + t * dy, pointA[Amg::z]};
        t = (-B - std::sqrt(det)) / (2 * A);
        intersectB = {pointA[Amg::x] + t * dx, pointA[Amg::y] + t * dy, pointB[Amg::z]};
        ATH_MSG_DEBUG("[circleLineIntersection2D] Two intersections at ("<<intersectA[Amg::x]<<","<<intersectA[Amg::y]<<","<<intersectA[Amg::z]<<") and at ("<<intersectB[Amg::x]<<","<<intersectB[Amg::y]<<","<<intersectB[Amg::z]<<").");
        return 2;
    }
 
 
}
 
 
#if defined(FLATTEN)
// We compile this package with optimization, even in debug builds; otherwise,
// the heavy use of Eigen makes it too slow.  However, from here we may call
// to out-of-line Eigen code that is linked from other DSOs; in that case,
// it would not be optimized.  Avoid this by forcing all Eigen code
// to be inlined here if possible.
ATH_FLATTEN
#endif
Amg::Vector3D FastCaloSimCaloExtrapolation::projectOnCylinder(float cylR, float cylZ, Amg::Vector3D& hitPos) {
        
  Amg::Vector3D closestPointOnCylinder;
  Amg::Vector3D cylAxis(0, 0, cylZ);
 
  //positive side
  if(hitPos[Amg::z] >= cylZ){
    //project hit position on x-y plane at positive side
    Amg::Vector3D projHitPos(hitPos[Amg::x], hitPos[Amg::y], cylZ);
    
    //if r of hit position outside cylinder, closest hit is always on edge  
    if(hitPos.perp() > cylR) closestPointOnCylinder = cylAxis + cylR * (projHitPos - cylAxis).unit();
    else closestPointOnCylinder = cylAxis + hitPos.perp() * (projHitPos - cylAxis).unit();
 
  } 
  //negative side
  else if (hitPos[Amg::z] <= -cylZ){
    //project hit position on x-y plane at negative side
    Amg::Vector3D projHitPos(hitPos[Amg::x], hitPos[Amg::y], -cylZ);
    
    if(hitPos.perp() > cylR) closestPointOnCylinder = -cylAxis + cylR * (projHitPos + cylAxis).unit();
    else closestPointOnCylinder = -cylAxis + hitPos.perp() * (projHitPos + cylAxis).unit();
 
  }
  else{
    Amg::Vector3D hitPosZ(0, 0, hitPos[Amg::z]);
    closestPointOnCylinder = hitPosZ + cylR * (hitPos - hitPosZ).unit(); 
  }
  
  return closestPointOnCylinder;
 
}
 
 
 
CylinderIntersections FastCaloSimCaloExtrapolation::getCylinderIntersections(float cylR, float cylZ, Amg::Vector3D& hitPos1, Amg::Vector3D& hitPos2) const{
  //calculates intersection of infinite line with cylinder --> can have 0 or 2 intersections
  CylinderIntersections intersections;
 
  //look for intersections with the cover of the cylinder
  unsigned int nCoverIntersections = cylinderLineIntersection(cylR, cylZ, hitPos1, hitPos2, intersections.first, intersections.second);
  if(nCoverIntersections == 2){
    ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found two cylinder intersections through cylinder cover.");
    intersections.number = 2;
    return intersections;
  }
  else if (nCoverIntersections == 1){
    
    Amg::Vector3D positiveEndcapIntersection, negativeEndcapIntersection;
    bool IsPositiveEndcapIntersection = cylinderEndcapIntersection(cylR, cylZ, true, hitPos1, hitPos2, positiveEndcapIntersection);
    bool IsNegativeEndcapIntersection = cylinderEndcapIntersection(cylR, cylZ, false, hitPos1, hitPos2, negativeEndcapIntersection);
 
    if(IsPositiveEndcapIntersection && IsNegativeEndcapIntersection){
      //if we have a cover intersection we only expect one additional endcap intersection
      //both endcap intersections can be valid in case the intersection is at the edge of the cylinder cover and endcap
      //in that case take the endcap intersection which is further away from the cylinder cover intersection to prevent taking the same intersection twice
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found intersection through cylinder cover and both endcaps. Intersection seems to be at edge of cover and endcap.");
      intersections.second = (positiveEndcapIntersection - intersections.first).norm() > (negativeEndcapIntersection - intersections.first).norm() ? positiveEndcapIntersection : negativeEndcapIntersection;
      intersections.number = 2;
    }
    else if(IsPositiveEndcapIntersection) {
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found intersection through cylinder cover and positive endcap.");
      intersections.second = positiveEndcapIntersection;
      intersections.number = 2;
    }
    else if(IsNegativeEndcapIntersection) {
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found intersection through cylinder cover and negative endcap.");
      intersections.second = negativeEndcapIntersection;
      intersections.number = 2;
    }
    else{
      //line is tangential to cylinder cover
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found single intersection through cylinder cover.");
      intersections.number = 1;
    }
 
  }
  else{
    //no cylinder cover intersections 
    Amg::Vector3D positiveEndcapIntersection, negativeEndcapIntersection;
    bool IsPositiveEndcapIntersection = cylinderEndcapIntersection(cylR, cylZ, true, hitPos1, hitPos2, positiveEndcapIntersection);
    bool IsNegativeEndcapIntersection = cylinderEndcapIntersection(cylR, cylZ, false, hitPos1, hitPos2, negativeEndcapIntersection);
 
    if(IsPositiveEndcapIntersection && IsNegativeEndcapIntersection){
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found intersections through both endcaps.");
      intersections.first = positiveEndcapIntersection;
      intersections.second = negativeEndcapIntersection;
      intersections.number = 2;
    }
    else if(IsPositiveEndcapIntersection) {
      //dont expect this to ever happen
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found single intersection through positive endcap. This should not happen.");
      intersections.first = positiveEndcapIntersection;
      intersections.number = 1;
    }
    else if(IsNegativeEndcapIntersection) {
      //dont expect this to ever happen
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found single intersection through negative endcap. This should not happen.");
      intersections.first = negativeEndcapIntersection;
      intersections.number = 1;
    }
    else{
      ATH_MSG_DEBUG("[getCylinderIntersections(R="<<cylR<<",Z="<<cylZ<<")] Found no cylinder intersections.");
      //no intersections at all
      intersections.number = 0;
 
    }
  }
 
  return intersections;
 
 
}
 
 
//calculates the intersection between the line defined by pointA and pointB and the cylinder cover definded by cylR and cylZ
int FastCaloSimCaloExtrapolation::cylinderLineIntersection(float cylR, float cylZ, Amg::Vector3D& pointA, Amg::Vector3D& pointB, Amg::Vector3D& intersectA, Amg::Vector3D& intersectB) const{
 
  //projections of points spanning the line onto the xy plane
  Amg::Vector3D projPointA(pointA[Amg::x], pointA[Amg::y], 0);
  Amg::Vector3D projPointB(pointB[Amg::x], pointB[Amg::y], 0);
  Amg::Vector3D projDiff = projPointA - projPointB;
 
  //calculate distance from (0,0,0) to line spanned by projPointA and projPointB
  double projDiffNorm2 = projDiff.dot(projDiff);
  double t  = projPointA.dot(projDiff) / projDiffNorm2;
  double d2 = projPointA.dot(projPointA) - t*t*projDiffNorm2;
 
  if(d2 < 0){
    ATH_MSG_COND("[cylinderLineIntersection] Got negative distance (d2="<<d2<<"). Forcing to 0.", d2 > -0.001);
    d2 = 0;
  }
 
  //if distance larger than cylinder radius then there are no intersection and we are done
  if(d2 > cylR*cylR) return 0;
 
  double k = std::sqrt((cylR*cylR - d2)/projDiffNorm2);
 
  intersectA = pointA + (t+k)*(pointB - pointA);
  intersectB = pointA + (t-k)*(pointB - pointA);
 
  //check if intersection is outside z bounds
  bool IntAisValid = (intersectA[Amg::z] <= cylZ && intersectA[Amg::z] >= -cylZ);
  bool IntBisValid = (intersectB[Amg::z] <= cylZ && intersectB[Amg::z] >= -cylZ);
 
  if(IntAisValid && IntBisValid) return 2;
  else if(IntAisValid) return 1;
  else if(IntBisValid){
    intersectA = intersectB;
    return 1;
  }
 
 
  return 0;
 
}
 
 
bool FastCaloSimCaloExtrapolation::cylinderEndcapIntersection(float cylR, float cylZ, bool positiveEndcap, Amg::Vector3D& pointA, Amg::Vector3D& pointB, Amg::Vector3D& intersection) {
 
  //normal and point on endcap defines the plane
  Amg::Vector3D pointOnEndcap;
  Amg::Vector3D normal(0, 0, 1);
  positiveEndcap ? pointOnEndcap = {0, 0, cylZ} : pointOnEndcap = {0, 0, -cylZ};
  Amg::Vector3D hitDir = (pointB - pointA);
 
  double denom = normal.dot(hitDir);
  if (std::abs(denom) > 1e-6) { 
    double t = normal.dot(pointOnEndcap - pointB)/denom;
    //compute intersection regardless of direction (t>0 or t<0) 
    intersection = pointB + t*hitDir;
    Amg::Vector3D v = intersection - pointOnEndcap; 
 
    //check if intersection is within cylR bounds 
    return std::sqrt(v.dot(v)) <= cylR;
 
  }
 
  return false;
 
  }
 
int FastCaloSimCaloExtrapolation::whichIntersection(float cylR, float cylZ, Amg::Vector3D& hitPos1, Amg::Vector3D& hitPos2, Amg::Vector3D& intersectionA, Amg::Vector3D intersectionB) const{
 
  //check if the hit positions are outside or inside the cylinder surface
  HITPOSITION cylPosHit1 = whereOnCylinder(cylR, cylZ, hitPos1);
  HITPOSITION cylPosHit2 = whereOnCylinder(cylR, cylZ, hitPos2);
  
  if((cylPosHit1 == INSIDE) ^ (cylPosHit2 == INSIDE)){
    /* CASE A: one hit position inside and one outside of the cylinder (travel through surface), 
    one intersection is on cylinder, take intersection closest to line segment */
    ATH_MSG_DEBUG("[whichIntersection] Travel through surface.");
    return getPointLineSegmentDistance(intersectionA, hitPos1, hitPos2) > getPointLineSegmentDistance(intersectionB, hitPos1, hitPos2);
  }
  else if(cylPosHit1 == INSIDE && cylPosHit2 == INSIDE){
    /* CASE B: both hit position inside, take intersection which points towards travel direction of particle */
    Amg::Vector3D directionA = intersectionA - hitPos2;
    Amg::Vector3D directionB = intersectionB - hitPos2;
    Amg::Vector3D hitDir = hitPos2 - hitPos1;
    ATH_MSG_DEBUG("[whichIntersection] Both hit positions inside.");
    return directionA.dot(hitDir) < directionB.dot(hitDir);
  }
  else{
    // /* CASE C: both hit position outside and the intersections lay on the segment, take intersection closest to second hit position */
    // /* CASE D: both hit positions are outside and the intersections are not on the line segment, take intersection closest to one of the hit positions */
    double distHitPosIntersectA = (hitPos2 - intersectionA).norm();
    double distHitPosIntersectB = (hitPos2 - intersectionB).norm();
    ATH_MSG_DEBUG("[whichIntersection] Both hit positions outside.");
    return distHitPosIntersectA > distHitPosIntersectB;
  }
} 
 
#if defined(FLATTEN)
// We compile this package with optimization, even in debug builds; otherwise,
// the heavy use of Eigen makes it too slow.  However, from here we may call
// to out-of-line Eigen code that is linked from other DSOs; in that case,
// it would not be optimized.  Avoid this by forcing all Eigen code
// to be inlined here if possible.
ATH_FLATTEN
#endif
double FastCaloSimCaloExtrapolation::getPointLineSegmentDistance(Amg::Vector3D& point, Amg::Vector3D& hitPos1, Amg::Vector3D& hitPos2) {
 
  Amg::Vector3D hitDir = hitPos2 - hitPos1;
  Amg::Vector3D w = point - hitPos1;
 
  double c1 = w.dot(hitDir);
  if(c1 <= 0) return Amg::distance(point, hitPos1);
  double c2 = hitDir.dot(hitDir);
  if(c2 <= c1) return Amg::distance(point, hitPos2);
  double t = c1/c2;
  Amg::Vector3D vec = hitPos1 + t*hitDir;
  return Amg::distance(point, vec);
 
}
 
enum FastCaloSimCaloExtrapolation::HITPOSITION FastCaloSimCaloExtrapolation::whereOnCylinder(float cylR, float cylZ, Amg::Vector3D& hitPos) {
  //set a 1mm tolerance within which the hit position is considered to be on the cylinder surface
  //setting this higher can lead to extrapolation failures around truth particle eta ~4
  float tolerance = 1;
 
  bool isOnEndcap = hitPos.perp() <= cylR + tolerance && (hitPos[Amg::z] > 0 ? std::abs(hitPos[Amg::z] - cylZ) < tolerance : std::abs(hitPos[Amg::z] + cylZ) < tolerance);
  bool isOnCover  = std::abs(hitPos.perp() - cylR) < tolerance && hitPos[Amg::z] < cylZ && hitPos[Amg::z] > -cylZ;
  
  //check if hit position is on endcap or cover of cylinder
  if(isOnEndcap || isOnCover) return HITPOSITION::ON;
  
  //check if hit position is inside cover
  if(hitPos[Amg::z] < cylZ && hitPos[Amg::z] > -cylZ && hitPos.perp() < cylR) return HITPOSITION::INSIDE;
  
  return HITPOSITION::OUTSIDE;
}
 
bool FastCaloSimCaloExtrapolation::doesTravelThroughSurface(float cylR, float cylZ, Amg::Vector3D& hitPos1, Amg::Vector3D& hitPos2) {
  //travel through surface in case one hit position is outside and the other outside of cylinder surface
  return (whereOnCylinder(cylR, cylZ, hitPos1) == INSIDE) ^ (whereOnCylinder(cylR, cylZ, hitPos2) == INSIDE);
}
 
bool FastCaloSimCaloExtrapolation::isOnSegment(Amg::Vector3D& point, Amg::Vector3D& hitPos1, Amg::Vector3D& hitPos2) {
  return getPointLineSegmentDistance(point, hitPos1, hitPos2) < 0.001;
}
 
bool FastCaloSimCaloExtrapolation::isCaloBarrel(int sample) const
{
  return GetCaloGeometry()->isCaloBarrel(sample);
}
 
double FastCaloSimCaloExtrapolation::deta(int sample, double eta) const
{
  return GetCaloGeometry()->deta(sample, eta);
}
 
void FastCaloSimCaloExtrapolation::minmaxeta(int sample, double eta, double& mineta, double& maxeta) const
{
  GetCaloGeometry()->minmaxeta(sample, eta, mineta, maxeta);
}
 
double FastCaloSimCaloExtrapolation::rmid(int sample, double eta) const
{
  return GetCaloGeometry()->rmid(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::zmid(int sample, double eta) const
{
  return GetCaloGeometry()->zmid(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::rzmid(int sample, double eta) const
{
  return GetCaloGeometry()->rzmid(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::rent(int sample, double eta) const
{
  return GetCaloGeometry()->rent(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::zent(int sample, double eta) const
{
  return GetCaloGeometry()->zent(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::rzent(int sample, double eta) const
{
  return GetCaloGeometry()->rzent(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::rext(int sample, double eta) const
{
  return GetCaloGeometry()->rext(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::zext(int sample, double eta) const
{
  return GetCaloGeometry()->zext(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::rzext(int sample, double eta) const
{
  return GetCaloGeometry()->rzext(sample, eta);
}
 
double FastCaloSimCaloExtrapolation::rpos(int sample, double eta, int subpos) const
{
  return GetCaloGeometry()->rpos(sample, eta, subpos);
}
 
double FastCaloSimCaloExtrapolation::zpos(int sample, double eta, int subpos) const
{
  return GetCaloGeometry()->zpos(sample, eta, subpos);
}
 
double FastCaloSimCaloExtrapolation::rzpos(int sample, double eta, int subpos) const
{
  return GetCaloGeometry()->rzpos(sample, eta, subpos);
}