FastCaloSimCaloExtrapolation Class Reference

#include <FastCaloSimCaloExtrapolation.h>

Inheritance diagram for FastCaloSimCaloExtrapolation:

Collaboration diagram for FastCaloSimCaloExtrapolation:

Public Types
enum	SUBPOS { SUBPOS_MID = TFCSExtrapolationState::SUBPOS_MID , SUBPOS_ENT = TFCSExtrapolationState::SUBPOS_ENT , SUBPOS_EXT = TFCSExtrapolationState::SUBPOS_EXT }
enum	HITPOSITION { INSIDE , OUTSIDE , ON }

Public Member Functions
	FastCaloSimCaloExtrapolation (const std::string &t, const std::string &n, const IInterface *p)
	~FastCaloSimCaloExtrapolation ()=default
virtual StatusCode	initialize () override final
virtual StatusCode	finalize () override final
virtual void	extrapolate (TFCSExtrapolationState &result, const TFCSTruthState *truth, const std::vector< G4FieldTrack > &caloSteps) const override final
virtual void	extrapolate (TFCSExtrapolationState &result, const TFCSTruthState *truth) const override final

Protected Member Functions
const IFastCaloSimGeometryHelper *	GetCaloGeometry () const
bool	extrapolateToCylinder (const std::vector< G4FieldTrack > &caloSteps, float cylR, float cylZ, Amg::Vector3D &extPos, Amg::Vector3D &momDir) const
	Finds best extrapolation extPos from the caloSteps for a cylinder defined by radius cylR and half-length cylZ as well as corresponding momentum direction.
void	extrapolateToID (TFCSExtrapolationState &result, const std::vector< G4FieldTrack > &caloSteps, const TFCSTruthState *truth) const
	Extrapolates to ID using three uniquely defined cylinder surfaces.
void	extrapolateToLayers (TFCSExtrapolationState &result, const std::vector< G4FieldTrack > &caloSteps, const TFCSTruthState *truth) const
	Extrapolates to all other layers of the calorimeter.
void	findPCA (float cylR, float cylZ, Amg::Vector3D &hitPos1, Amg::Vector3D &hitPos2, Amg::Vector3D &PCA) const
	Finds Point of Closest Approach (PCA) on the cylinder defined by radius cylR and half-length cylZ of a line segment spanned by two hit positions to a cylinder.
void	getIterativePCA (float cylR, float cylZ, Amg::Vector3D &BoundA, Amg::Vector3D &BoundB, Amg::Vector3D &PCA) const
	Finds PCA iteratively given two bounds A and B on a line segment, used for (rare) cases with no easy analytical solutions.
bool	extrapolateWithIntersection (const std::vector< G4FieldTrack > &caloSteps, float cylR, float cylZ, Amg::Vector3D &extPos, Amg::Vector3D &momDir) const
	Extrapolates position on cylinder by finding intersections of subsequent hit positions, intersection is considered if we detect a travel through the surface with the line segment or we find a forward intersection (in the travel direction of the particle) which lies on the line segment, returns false if no such postion is found.
bool	extrapolateWithPCA (const std::vector< G4FieldTrack > &caloSteps, float cylR, float cylZ, Amg::Vector3D &extPos, Amg::Vector3D &momDir) const
	Extrapolates to the cylinder using the PCA to the polygon spanned by the individual line segments from the caloSteps.
int	whichIntersection (float cylR, float cylZ, Amg::Vector3D &hitPos1, Amg::Vector3D &hitPos2, Amg::Vector3D &intersectionA, Amg::Vector3D intersectionB) const
	Returns ID of more sensible intersection between line segment spanned by hitPos1 and hitPos2 and cylinder.
int	circleLineIntersection2D (float circR, Amg::Vector3D &pointA, Amg::Vector3D &pointB, Amg::Vector3D &intersectA, Amg::Vector3D &intersectB) const
	Analytically computes 2D intersections between circle of radius circR and (infinite) line spanned by pointA nad pointB.
int	cylinderLineIntersection (float cylR, float cylZ, Amg::Vector3D &pointA, Amg::Vector3D &pointB, Amg::Vector3D &intersectA, Amg::Vector3D &intersectB) const
	Analytically computes the intersection between the (infinite) line defined by pointA and pointB and the cylinder cover (without endcaps)
CylinderIntersections	getCylinderIntersections (float cylR, float cylZ, Amg::Vector3D &hitPos1, Amg::Vector3D &hitPos2) const
	Analytically computes the intersection between the (infinite) line spanned by hitPos1 and hitPos2 with a cylinder.
void	minmaxeta (int sample, double eta, double &mineta, double &maxeta) const
bool	isCaloBarrel (int sample) const
double	deta (int sample, double eta) const
double	rzmid (int sample, double eta) const
double	rzent (int sample, double eta) const
double	rzext (int sample, double eta) const
double	rmid (int sample, double eta) const
double	rent (int sample, double eta) const
double	rext (int sample, double eta) const
double	zmid (int sample, double eta) const
double	zent (int sample, double eta) const
double	zext (int sample, double eta) const
double	rpos (int sample, double eta, int subpos=CaloSubPos::SUBPOS_MID) const
double	zpos (int sample, double eta, int subpos=CaloSubPos::SUBPOS_MID) const
double	rzpos (int sample, double eta, int subpos=CaloSubPos::SUBPOS_MID) const

Static Protected Member Functions
static double	getPointLineSegmentDistance (Amg::Vector3D &point, Amg::Vector3D &hitPos1, Amg::Vector3D &hitPos2)
	Computes the distance between a point and the line segment spanned by hitPos1 and hitPos2.
static bool	isOnSegment (Amg::Vector3D &point, Amg::Vector3D &hitPos1, Amg::Vector3D &hitPos2)
	Returns true if point lies on the line segment spanned by hitPos1 and hitPos2, otherwise returns false.
static bool	cylinderEndcapIntersection (float cylR, float cylZ, bool positiveEndcap, Amg::Vector3D &pointA, Amg::Vector3D &pointB, Amg::Vector3D &intersection)
	Computes intersection between the (infinite) line spanned by pointA and pointB with the positive (negative) endcap of a cylinder, returns true if intersection is found.
static bool	doesTravelThroughSurface (float cylR, float cylZ, Amg::Vector3D &hitPos1, Amg::Vector3D &hitPos2)
	Returns true if the line segment spanned by hitPos1 and hitPos2 crosses the cylinder surface, false otherwise.
static enum HITPOSITION	whereOnCylinder (float cylR, float cylZ, Amg::Vector3D &hitPos)
	Checks if position of hitPos is inside, outside or on the cylinder bounds.
static Amg::Vector3D	projectOnCylinder (float cylR, float cylZ, Amg::Vector3D &hitPos)
	Projects position hitPos onto the cylinder surface and returns projected position.

Protected Attributes
HepPDT::ParticleDataTable *	m_particleDataTable {nullptr}
FloatArrayProperty	m_CaloBoundaryR {this, "CaloBoundaryR", {1148.0,120.0,41.0}}
FloatArrayProperty	m_CaloBoundaryZ {this, "CaloBoundaryZ", {3550.0,4587.0,4587.0}}
PublicToolHandle< IFastCaloSimCaloTransportation >	m_CaloTransportation {this, "CaloTransportation", "FastCaloSimCaloTransportation"}
PublicToolHandle< IFastCaloSimGeometryHelper >	m_CaloGeometryHelper {this, "CaloGeometryHelper", "FastCaloSimGeometryHelper"}

Detailed Description

Definition at line 31 of file FastCaloSimCaloExtrapolation.h.

Member Enumeration Documentation

HITPOSITION

enum FastCaloSimCaloExtrapolation::HITPOSITION

Enumerator
INSIDE
OUTSIDE
ON

Definition at line 48 of file FastCaloSimCaloExtrapolation.h.

                  {
    INSIDE,  //hit position is inside cylinder bounds
    OUTSIDE, //hit position is outside cylinder bounds
    ON       //hit position is on cylinder bounds
  };

SUBPOS

enum FastCaloSimCaloExtrapolation::SUBPOS

Enumerator
SUBPOS_MID
SUBPOS_ENT
SUBPOS_EXT

Definition at line 42 of file FastCaloSimCaloExtrapolation.h.

              { 
    SUBPOS_MID = TFCSExtrapolationState::SUBPOS_MID, //MID=middle of calo layer
    SUBPOS_ENT = TFCSExtrapolationState::SUBPOS_ENT, //ENT=entrance of calo layer
    SUBPOS_EXT = TFCSExtrapolationState::SUBPOS_EXT  //EXT=exit of calo layer
  }; 

Constructor & Destructor Documentation

FastCaloSimCaloExtrapolation()

FastCaloSimCaloExtrapolation::FastCaloSimCaloExtrapolation	(	const std::string &	t,
		const std::string &	n,
		const IInterface *	p )

Definition at line 47 of file FastCaloSimCaloExtrapolation.cxx.

  : base_class(t,n,p)
{
}

~FastCaloSimCaloExtrapolation()

FastCaloSimCaloExtrapolation::~FastCaloSimCaloExtrapolation ( )

default

Member Function Documentation

circleLineIntersection2D()

int FastCaloSimCaloExtrapolation::circleLineIntersection2D ( float circR, Amg::Vector3D & pointA, Amg::Vector3D & pointB, Amg::Vector3D & intersectA, Amg::Vector3D & intersectB ) const

protected

Analytically computes 2D intersections between circle of radius circR and (infinite) line spanned by pointA nad pointB.

Definition at line 565 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                       {
    //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;
    }
 
 
}

cylinderEndcapIntersection()

bool FastCaloSimCaloExtrapolation::cylinderEndcapIntersection ( float cylR, float cylZ, bool positiveEndcap, Amg::Vector3D & pointA, Amg::Vector3D & pointB, Amg::Vector3D & intersection )

staticprotected

Computes intersection between the (infinite) line spanned by pointA and pointB with the positive (negative) endcap of a cylinder, returns true if intersection is found.

Definition at line 772 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                              {
 
  //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;
 
  }

cylinderLineIntersection()

int FastCaloSimCaloExtrapolation::cylinderLineIntersection ( float cylR, float cylZ, Amg::Vector3D & pointA, Amg::Vector3D & pointB, Amg::Vector3D & intersectA, Amg::Vector3D & intersectB ) const

protected

Analytically computes the intersection between the (infinite) line defined by pointA and pointB and the cylinder cover (without endcaps)

Definition at line 730 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                                  {
 
  //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;
 
}

deta()

double FastCaloSimCaloExtrapolation::deta ( int sample, double eta ) const

protected

Definition at line 880 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->deta(sample, eta);
}

doesTravelThroughSurface()

bool FastCaloSimCaloExtrapolation::doesTravelThroughSurface ( float cylR, float cylZ, Amg::Vector3D & hitPos1, Amg::Vector3D & hitPos2 )

staticprotected

Returns true if the line segment spanned by hitPos1 and hitPos2 crosses the cylinder surface, false otherwise.

Definition at line 866 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                              {
  //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);
}

extrapolate() [1/2]

void FastCaloSimCaloExtrapolation::extrapolate ( TFCSExtrapolationState & result, const TFCSTruthState * truth ) const

finaloverridevirtual

Definition at line 84 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                               { 
  
  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);
}

extrapolate() [2/2]

void FastCaloSimCaloExtrapolation::extrapolate ( TFCSExtrapolationState & result, const TFCSTruthState * truth, const std::vector< G4FieldTrack > & caloSteps ) const

finaloverridevirtual

Definition at line 72 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                         { 
  
  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");
 
}

extrapolateToCylinder()

bool FastCaloSimCaloExtrapolation::extrapolateToCylinder ( const std::vector< G4FieldTrack > & caloSteps, float cylR, float cylZ, Amg::Vector3D & extPos, Amg::Vector3D & momDir ) const

protected

Finds best extrapolation extPos from the caloSteps for a cylinder defined by radius cylR and half-length cylZ as well as corresponding momentum direction.

Definition at line 237 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                        {
 
  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;
 
}

extrapolateToID()

void FastCaloSimCaloExtrapolation::extrapolateToID ( TFCSExtrapolationState & result, const std::vector< G4FieldTrack > & caloSteps, const TFCSTruthState * truth ) const

protected

Extrapolates to ID using three uniquely defined cylinder surfaces.

Definition at line 94 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                             {
 
  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()");
 
} 

extrapolateToLayers()

void FastCaloSimCaloExtrapolation::extrapolateToLayers ( TFCSExtrapolationState & result, const std::vector< G4FieldTrack > & caloSteps, const TFCSTruthState * truth ) const

protected

Extrapolates to all other layers of the calorimeter.

Definition at line 166 of file FastCaloSimCaloExtrapolation.cxx.

{
  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()");
}

extrapolateWithIntersection()

bool FastCaloSimCaloExtrapolation::extrapolateWithIntersection ( const std::vector< G4FieldTrack > & caloSteps, float cylR, float cylZ, Amg::Vector3D & extPos, Amg::Vector3D & momDir ) const

protected

Extrapolates position on cylinder by finding intersections of subsequent hit positions, intersection is considered if we detect a travel through the surface with the line segment or we find a forward intersection (in the travel direction of the particle) which lies on the line segment, returns false if no such postion is found.

Definition at line 264 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                              {
 
  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;
}

extrapolateWithPCA()

bool FastCaloSimCaloExtrapolation::extrapolateWithPCA ( const std::vector< G4FieldTrack > & caloSteps, float cylR, float cylZ, Amg::Vector3D & extPos, Amg::Vector3D & momDir ) const

protected

Extrapolates to the cylinder using the PCA to the polygon spanned by the individual line segments from the caloSteps.

Definition at line 342 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                     {
 
  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;
} 

finalize()

StatusCode FastCaloSimCaloExtrapolation::finalize ( )

finaloverridevirtual

Definition at line 66 of file FastCaloSimCaloExtrapolation.cxx.

                                                 {
  ATH_MSG_INFO( "Finalizing FastCaloSimCaloExtrapolation" );
  return StatusCode::SUCCESS;
}

findPCA()

void FastCaloSimCaloExtrapolation::findPCA ( float cylR, float cylZ, Amg::Vector3D & hitPos1, Amg::Vector3D & hitPos2, Amg::Vector3D & PCA ) const

protected

Finds Point of Closest Approach (PCA) on the cylinder defined by radius cylR and half-length cylZ of a line segment spanned by two hit positions to a cylinder.

Definition at line 382 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                    {
  //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);
  
}

GetCaloGeometry()

const IFastCaloSimGeometryHelper * FastCaloSimCaloExtrapolation::GetCaloGeometry ( ) const

inlineprotected

Definition at line 59 of file FastCaloSimCaloExtrapolation.h.

{return &(*m_CaloGeometryHelper);};

getCylinderIntersections()

CylinderIntersections FastCaloSimCaloExtrapolation::getCylinderIntersections ( float cylR, float cylZ, Amg::Vector3D & hitPos1, Amg::Vector3D & hitPos2 ) const

protected

Analytically computes the intersection between the (infinite) line spanned by hitPos1 and hitPos2 with a cylinder.

Definition at line 649 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                    {
  //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;
 
 
}

getIterativePCA()

void FastCaloSimCaloExtrapolation::getIterativePCA ( float cylR, float cylZ, Amg::Vector3D & BoundA, Amg::Vector3D & BoundB, Amg::Vector3D & PCA ) const

protected

Finds PCA iteratively given two bounds A and B on a line segment, used for (rare) cases with no easy analytical solutions.

Definition at line 488 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                          {
 
    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);
 
}

getPointLineSegmentDistance()

double FastCaloSimCaloExtrapolation::getPointLineSegmentDistance ( Amg::Vector3D & point, Amg::Vector3D & hitPos1, Amg::Vector3D & hitPos2 )

staticprotected

Computes the distance between a point and the line segment spanned by hitPos1 and hitPos2.

Definition at line 834 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                               {
 
  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);
 
}

initialize()

StatusCode FastCaloSimCaloExtrapolation::initialize ( )

finaloverridevirtual

Definition at line 52 of file FastCaloSimCaloExtrapolation.cxx.

{
  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;
 
 
}

isCaloBarrel()

bool FastCaloSimCaloExtrapolation::isCaloBarrel ( int sample ) const

protected

Definition at line 875 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->isCaloBarrel(sample);
}

isOnSegment()

bool FastCaloSimCaloExtrapolation::isOnSegment ( Amg::Vector3D & point, Amg::Vector3D & hitPos1, Amg::Vector3D & hitPos2 )

staticprotected

Returns true if point lies on the line segment spanned by hitPos1 and hitPos2, otherwise returns false.

Definition at line 871 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                             {
  return getPointLineSegmentDistance(point, hitPos1, hitPos2) < 0.001;
}

minmaxeta()

void FastCaloSimCaloExtrapolation::minmaxeta ( int sample, double eta, double & mineta, double & maxeta ) const

protected

Definition at line 885 of file FastCaloSimCaloExtrapolation.cxx.

{
  GetCaloGeometry()->minmaxeta(sample, eta, mineta, maxeta);
}

projectOnCylinder()

Amg::Vector3D FastCaloSimCaloExtrapolation::projectOnCylinder ( float cylR, float cylZ, Amg::Vector3D & hitPos )

staticprotected

Projects position hitPos onto the cylinder surface and returns projected position.

Definition at line 614 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                       {
        
  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;
 
}

rent()

double FastCaloSimCaloExtrapolation::rent ( int sample, double eta ) const

protected

Definition at line 905 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rent(sample, eta);
}

rext()

double FastCaloSimCaloExtrapolation::rext ( int sample, double eta ) const

protected

Definition at line 920 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rext(sample, eta);
}

rmid()

double FastCaloSimCaloExtrapolation::rmid ( int sample, double eta ) const

protected

Definition at line 890 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rmid(sample, eta);
}

rpos()

double FastCaloSimCaloExtrapolation::rpos ( int sample, double eta, int subpos = CaloSubPos::SUBPOS_MID ) const

protected

Definition at line 935 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rpos(sample, eta, subpos);
}

rzent()

double FastCaloSimCaloExtrapolation::rzent ( int sample, double eta ) const

protected

Definition at line 915 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rzent(sample, eta);
}

rzext()

double FastCaloSimCaloExtrapolation::rzext ( int sample, double eta ) const

protected

Definition at line 930 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rzext(sample, eta);
}

rzmid()

double FastCaloSimCaloExtrapolation::rzmid ( int sample, double eta ) const

protected

Definition at line 900 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rzmid(sample, eta);
}

rzpos()

double FastCaloSimCaloExtrapolation::rzpos ( int sample, double eta, int subpos = CaloSubPos::SUBPOS_MID ) const

protected

Definition at line 945 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->rzpos(sample, eta, subpos);
}

whereOnCylinder()

enum FastCaloSimCaloExtrapolation::HITPOSITION FastCaloSimCaloExtrapolation::whereOnCylinder ( float cylR, float cylZ, Amg::Vector3D & hitPos )

staticprotected

Checks if position of hitPos is inside, outside or on the cylinder bounds.

Definition at line 849 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                      {
  //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;
}

whichIntersection()

int FastCaloSimCaloExtrapolation::whichIntersection ( float cylR, float cylZ, Amg::Vector3D & hitPos1, Amg::Vector3D & hitPos2, Amg::Vector3D & intersectionA, Amg::Vector3D intersectionB ) const

protected

Returns ID of more sensible intersection between line segment spanned by hitPos1 and hitPos2 and cylinder.

Definition at line 796 of file FastCaloSimCaloExtrapolation.cxx.

                                                                                                                                                                                  {
 
  //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;
  }
} 

zent()

double FastCaloSimCaloExtrapolation::zent ( int sample, double eta ) const

protected

Definition at line 910 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->zent(sample, eta);
}

zext()

double FastCaloSimCaloExtrapolation::zext ( int sample, double eta ) const

protected

Definition at line 925 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->zext(sample, eta);
}

zmid()

double FastCaloSimCaloExtrapolation::zmid ( int sample, double eta ) const

protected

Definition at line 895 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->zmid(sample, eta);
}

zpos()

double FastCaloSimCaloExtrapolation::zpos ( int sample, double eta, int subpos = CaloSubPos::SUBPOS_MID ) const

protected

Definition at line 940 of file FastCaloSimCaloExtrapolation.cxx.

{
  return GetCaloGeometry()->zpos(sample, eta, subpos);
}

Member Data Documentation

m_CaloBoundaryR

FloatArrayProperty FastCaloSimCaloExtrapolation::m_CaloBoundaryR {this, "CaloBoundaryR", {1148.0,120.0,41.0}}

protected

Definition at line 122 of file FastCaloSimCaloExtrapolation.h.

{this, "CaloBoundaryR", {1148.0,120.0,41.0}};

m_CaloBoundaryZ

FloatArrayProperty FastCaloSimCaloExtrapolation::m_CaloBoundaryZ {this, "CaloBoundaryZ", {3550.0,4587.0,4587.0}}

protected

Definition at line 123 of file FastCaloSimCaloExtrapolation.h.

{this, "CaloBoundaryZ", {3550.0,4587.0,4587.0}};

m_CaloGeometryHelper

PublicToolHandle<IFastCaloSimGeometryHelper> FastCaloSimCaloExtrapolation::m_CaloGeometryHelper {this, "CaloGeometryHelper", "FastCaloSimGeometryHelper"}

protected

Definition at line 128 of file FastCaloSimCaloExtrapolation.h.

{this, "CaloGeometryHelper", "FastCaloSimGeometryHelper"};

m_CaloTransportation

PublicToolHandle<IFastCaloSimCaloTransportation> FastCaloSimCaloExtrapolation::m_CaloTransportation {this, "CaloTransportation", "FastCaloSimCaloTransportation"}

protected

Definition at line 126 of file FastCaloSimCaloExtrapolation.h.

{this, "CaloTransportation", "FastCaloSimCaloTransportation"};

m_particleDataTable

HepPDT::ParticleDataTable* FastCaloSimCaloExtrapolation::m_particleDataTable {nullptr}

protected

Definition at line 119 of file FastCaloSimCaloExtrapolation.h.

{nullptr};

---

The documentation for this class was generated from the following files:

• FastCaloSimCaloExtrapolation.h
• FastCaloSimCaloExtrapolation.cxx