ExtrapolationValidation.cxx

Go to the documentation of this file.

/*
   Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
 */

// ExtrapolationValidation.cxx, (c) ATLAS Detector software
 
// Tracking
#include <cmath>
 
#include "TrkExAlgs/ExtrapolationValidation.h"
#include "TrkSurfaces/CylinderSurface.h"
#include "TrkSurfaces/DiscSurface.h"
#include "TrkSurfaces/PlaneSurface.h"
#include "TrkEventPrimitives/ParticleHypothesis.h"
#include "TrkGeometry/TrackingGeometry.h"
#include "TrkGeometry/TrackingVolume.h"
#include "TrkVolumes/CylinderVolumeBounds.h"
#include "TrkTrack/TrackStateOnSurface.h"
// Validation mode - TTree includes
#include "TTree.h"
#include "GaudiKernel/ITHistSvc.h"
 
//================ Constructor =================================================
 
Trk::ExtrapolationValidation::ExtrapolationValidation(const std::string& name, ISvcLocator* pSvcLocator)
  :
  AthAlgorithm(name, pSvcLocator) {}
 
//================ Destructor =================================================
 
Trk::ExtrapolationValidation::~ExtrapolationValidation() {
  // clear random number generators
  delete m_gaussDist;
  delete m_flatDist;
}
 
//================ Initialisation =================================================
 
StatusCode Trk::ExtrapolationValidation::initialize() {
  // Code entered here will be executed once at program start.
  ATH_MSG_INFO(" initialize()");
 
  // Get Extrapolator from ToolService
  ATH_CHECK(m_extrapolator.retrieve());
 
  // create the new Tree
  m_validationTree = new TTree(m_validationTreeName.value().c_str(),
                               m_validationTreeDescription.value().c_str());
 
  // the branches for the parameters
  m_validationTree->Branch("Parameters", &m_parameters, "params/I");
  m_validationTree->Branch("ParametersLoc1", m_parameterLoc1, "paramLoc1[params]/F");
  m_validationTree->Branch("ParametersLoc2", m_parameterLoc2, "paramLoc2[params]/F");
  m_validationTree->Branch("ParametersPhi", m_parameterPhi, "paramPhi[params]/F");
  m_validationTree->Branch("ParametersTheta", m_parameterTheta, "paramTheta[params]/F");
  m_validationTree->Branch("ParametersEta", m_parameterEta, "paramEta[params]/F");
  m_validationTree->Branch("ParametersQoverP", m_parameterQoverP, "paramQoverP[params]/F");
  // for the covariance diagonals
  m_validationTree->Branch("CovarianceLoc1", m_covarianceLoc1, "covLoc1[params]/F");
  m_validationTree->Branch("CovarianceLoc2", m_covarianceLoc2, "covLoc2[params]/F");
  m_validationTree->Branch("CovariancePhi", m_covariancePhi, "covPhi[params]/F");
  m_validationTree->Branch("CovarianceTheta", m_covarianceTheta, "covTheta[params]/F");
  m_validationTree->Branch("CovarianceQoverP", m_covarianceQoverP, "covQoverP[params]/F");
  m_validationTree->Branch("CovarianceDeterminant", m_covarianceDeterminant, "covDet[params]/F");
  // the start Momentum
  m_validationTree->Branch("StartMomentum", &m_startP, "startP/F");
  // for the start surface
  m_validationTree->Branch("StartSurfaceX", &m_startX, "startX/F");
  m_validationTree->Branch("StartSurfaceY", &m_startY, "startY/F");
  m_validationTree->Branch("StartSurfaceR", &m_startR, "startR/F");
  m_validationTree->Branch("StartSurfaceZ", &m_startZ, "startZ/F");
  // the estimation of the parameters
  m_validationTree->Branch("EstimationSurfaceX", &m_estimationX, "estimateX/F");
  m_validationTree->Branch("EstimationSurfaceY", &m_estimationY, "estimateY/F");
  m_validationTree->Branch("EstimationSurfaceR", &m_estimationR, "estimateR/F");
  m_validationTree->Branch("EstimationSurfaceZ", &m_estimationZ, "estimateZ/F");
  // for the surface type
  m_validationTree->Branch("DestinationSurfaceType", &m_destinationSurfaceType, "surfaceType/I");
  m_validationTree->Branch("DestinationSurfaceX", &m_destinationX, "surfaceX/F");
  m_validationTree->Branch("DestinationSurfaceY", &m_destinationY, "surfaceY/F");
  m_validationTree->Branch("DestinationSurfaceR", &m_destinationR, "surfaceR/F");
  m_validationTree->Branch("DestinationSurfaceZ", &m_destinationZ, "surfaceZ/F");
 
  // now register the Tree
  SmartIF<ITHistSvc> tHistSvc {
    service("THistSvc")
  };
  if (!tHistSvc) {
    ATH_MSG_ERROR("initialize() Could not find Hist Service -> Switching ValidationMode Off !");
    delete m_validationTree;
    m_validationTree = nullptr;
  }
  if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
    ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !");
    delete m_validationTree;
    m_validationTree = nullptr;
  }
 
  // intialize the random number generators
  m_gaussDist = new Rndm::Numbers(randSvc(), Rndm::Gauss(0., 1.));
  m_flatDist = new Rndm::Numbers(randSvc(), Rndm::Flat(0., 1.));
 
  ATH_MSG_INFO("initialize() successful");
  return StatusCode::SUCCESS;
}
 
//================ Finalisation =================================================
 
StatusCode Trk::ExtrapolationValidation::finalize() {
  // Code entered here will be executed once at the end of the program run.
  ATH_MSG_INFO("================== Output Statistics =========================");
  ATH_MSG_INFO("= Navigation : ");
  ATH_MSG_INFO("=  - breaks fwd : " << static_cast<double>(m_breaksFront) / static_cast<double>(m_triesFront)
                                    << " (" << m_breaksFront << "/" << m_triesFront << ")");
  ATH_MSG_INFO("=  - breaks bwd : " << static_cast<double>(m_breaksBack) / static_cast<double>(m_triesBack)
                                    << " (" << m_breaksBack << "/" << m_triesBack << ")");
  if (m_materialCollectionValidation) {
    ATH_MSG_INFO("= Material collection : ");
    ATH_MSG_INFO("=  - layer collected fwd : " << m_collectedLayerFront);
    ATH_MSG_INFO("=  - layer collected bwd : " << m_collectedLayerBack);
  }
 
  ATH_MSG_INFO("==============================================================");
 
  return StatusCode::SUCCESS;
}
 
//================ Execution ====================================================
 
StatusCode Trk::ExtrapolationValidation::execute(const EventContext& ctx) {
  // get the overall dimensions
  if (!m_highestVolume) {
    // get TrackingGeometry and highest volume
    const Trk::TrackingGeometry* trackingGeometry = m_extrapolator->trackingGeometry();
    m_highestVolume = trackingGeometry ? trackingGeometry->highestTrackingVolume() : nullptr;
    const Trk::CylinderVolumeBounds* cylBounds = m_highestVolume ?
                                                 dynamic_cast<const Trk::CylinderVolumeBounds*>(&(m_highestVolume->
                                                                                                   volumeBounds())) :
                                                 nullptr;
    // bail out
    if (!cylBounds) {
      ATH_MSG_WARNING("No highest TrackingVolume / no VolumeBounds ... pretty useless! ");
      return StatusCode::SUCCESS;
    }
    // get the numbers
    m_maximumR = cylBounds->outerRadius();
    m_maximumZ = cylBounds->halflengthZ();
  }
 
  // intialize the values
  m_parameters = 0;
  m_destinationSurfaceType = 0;
  // -----------> start
  m_startX = 0.;
  m_startY = 0.;
  m_startR = 0.;
  m_startZ = 0.;
  // -----------> estimation
  m_estimationX = 0.;
  m_estimationY = 0.;
  m_estimationR = 0.;
  m_estimationZ = 0.;
  // -----------> destination
  m_destinationX = 0.;
  m_destinationY = 0.;
  m_destinationR = 0.;
  m_destinationZ = 0.;
 
  // the local start parameters
  m_parameterLoc1[m_parameters] = m_sigmaLoc * m_gaussDist->shoot();
  m_parameterLoc2[m_parameters] = m_sigmaLoc * m_gaussDist->shoot();
  // are adopted for planar and straight line surfaces
  m_parameterPhi[m_parameters] = M_PI * m_flatDist->shoot();
  m_parameterPhi[m_parameters] *= (m_flatDist->shoot() > 0.5) ? -1. : 1.;
  m_parameterEta[m_parameters] = m_minEta + m_flatDist->shoot() * (m_maxEta - m_minEta);
  m_parameterTheta[m_parameters] = 2. * std::atan(std::exp(-m_parameterEta[m_parameters]));
 
  m_covarianceLoc1[m_parameters] = std::abs(m_parameterLoc1[m_parameters] * 0.1);
  m_covarianceLoc2[m_parameters] = std::abs(m_parameterLoc2[m_parameters] * 0.1);
  m_covariancePhi[m_parameters] = std::abs(m_parameterPhi[m_parameters] * 0.1);
  m_covarianceTheta[m_parameters] = std::abs(m_parameterTheta[m_parameters] * 0.1);
 
  // this is fine
  double p = m_minP + m_flatDist->shoot() * (m_maxP - m_minP);
  double charge = (m_flatDist->shoot() > 0.5) ? -1. : 1.;
  m_parameterQoverP[m_parameters] = charge / p;
 
  m_covarianceQoverP[m_parameters] = std::abs(m_parameterQoverP[m_parameters] * 0.1);
 
  // for the momentum logging
  m_startP = p;
 
  // start
  m_startR = std::abs(m_sigmaR * m_gaussDist->shoot());
  double surfacePhi = M_PI * m_flatDist->shoot();
  surfacePhi *= (m_flatDist->shoot() > 0.5) ? -1. : 1.;
  m_startX = m_startR * cos(surfacePhi);
  m_startY = m_startR * sin(surfacePhi);
  m_startZ = m_sigmaZ * m_gaussDist->shoot();
 
  // rotate it around Z
  double alphaZ = M_PI * m_flatDist->shoot();
  alphaZ *= (m_flatDist->shoot() > 0.5) ? -1. : 1.;
 
  // create the plane surface
  Trk::PlaneSurface startSurface(createTransform(m_startX,
                                                 m_startY,
                                                 m_startZ,
                                                 m_parameterPhi[m_parameters],
                                                 m_parameterTheta[m_parameters],
                                                 alphaZ),
                                 10e3, 10e3);
 
 
 
  AmgSymMatrix(5) covariance;
  covariance.setZero();
  ATH_MSG_VERBOSE(m_covarianceLoc1[m_parameters]);
  covariance(0, 0) = m_covarianceLoc1[m_parameters];
  ATH_MSG_VERBOSE(m_covarianceLoc2[m_parameters]);
  covariance(1, 1) = m_covarianceLoc2[m_parameters];
  ATH_MSG_VERBOSE(m_covariancePhi[m_parameters]);
  covariance(2, 2) = m_covariancePhi[m_parameters];
  ATH_MSG_VERBOSE(m_covarianceTheta[m_parameters]);
  covariance(3, 3) = m_covarianceTheta[m_parameters];
  ATH_MSG_VERBOSE(m_covarianceQoverP[m_parameters]);
  covariance(4, 4) = m_covarianceQoverP[m_parameters];
  ATH_MSG_VERBOSE("Initial Setting: \n" << covariance);
 
 
  m_covarianceDeterminant[m_parameters] = covariance.determinant();
 
  // the initial perigee with random numbers
  Trk::AtaPlane startParameters(m_parameterLoc1[m_parameters],
                                m_parameterLoc2[m_parameters],
                                m_parameterPhi[m_parameters],
                                m_parameterTheta[m_parameters],
                                m_parameterQoverP[m_parameters],
                                startSurface,
                                std::move(covariance));
 
  ATH_MSG_VERBOSE("Start Parameters : " << startParameters);
  if (startParameters.covariance()) ATH_MSG_VERBOSE("Start Covariance : \n" << *startParameters.covariance());
 
 
  // destination position
  m_estimationR = m_maximumR * m_flatDist->shoot();
 
  // --------------- propagate to find a first intersection ---------------------
  Amg::Transform3D CylTrf;
  CylTrf.setIdentity();
  Trk::CylinderSurface estimationCylinder(CylTrf, m_estimationR, 10e10);
  const Trk::TrackParameters* estimationParameters = m_extrapolator->extrapolateDirectly(ctx,
                                                                                         startParameters,
                                                                                         estimationCylinder,
                                                                                         Trk::alongMomentum,
                                                                                         false).release();
  if (!estimationParameters) {
    ATH_MSG_VERBOSE("Estimation of intersection did not work - skip event !");
    return StatusCode::SUCCESS;
  } else if (m_highestVolume && estimationParameters && !(m_highestVolume->inside(estimationParameters->position()))) {
    ATH_MSG_VERBOSE("Estimation of intersection is outside the known world - skip event !");
    delete estimationParameters;
    return StatusCode::SUCCESS;
  }
 
  ATH_MSG_VERBOSE("Estimation Parameters: " << *estimationParameters);
 
  // record the estimation parameters
  ++m_triesFront;
  ++m_parameters;
  m_parameterLoc1[m_parameters] = estimationParameters->parameters()[Trk::loc1];
  m_parameterLoc2[m_parameters] = estimationParameters->parameters()[Trk::loc2];
  m_parameterPhi[m_parameters] = estimationParameters->parameters()[Trk::phi];
  m_parameterEta[m_parameters] = estimationParameters->eta();
  m_parameterTheta[m_parameters] = estimationParameters->parameters()[Trk::theta];
  m_parameterQoverP[m_parameters] = estimationParameters->parameters()[Trk::qOverP];
  if (estimationParameters->covariance()) {
    m_covarianceLoc1[m_parameters] = (*estimationParameters->covariance())(0, 0);
    m_covarianceLoc2[m_parameters] = (*estimationParameters->covariance())(1, 1);
    m_covariancePhi[m_parameters] = (*estimationParameters->covariance())(2, 2);
    m_covarianceTheta[m_parameters] = (*estimationParameters->covariance())(3, 3);
    m_covarianceQoverP[m_parameters] = (*estimationParameters->covariance())(4, 4);
    m_covarianceDeterminant[m_parameters] = (estimationParameters->covariance())->determinant();
  } else {
    m_covarianceLoc1[m_parameters] = 0;
    m_covarianceLoc2[m_parameters] = 0;
    m_covariancePhi[m_parameters] = 0;
    m_covarianceTheta[m_parameters] = 0;
    m_covarianceQoverP[m_parameters] = 0;
    m_covarianceDeterminant[m_parameters] = 0;
  }
  // the start Momentum
 
  // get the estimated position
  const Amg::Vector3D& estimatedPosition = estimationParameters->position();
 
  m_estimationX = estimatedPosition.x();
  m_estimationY = estimatedPosition.y();
  m_estimationZ = estimatedPosition.z();
 
  // cleanup for memory reasons
  delete estimationParameters;
  estimationParameters = nullptr;
 
  // create the radom surface at the destination point
  Trk::PlaneSurface destinationSurface(createTransform(m_estimationX,
                                                       m_estimationY,
                                                       m_estimationZ,
                                                       m_parameterPhi[m_parameters],
                                                       m_parameterTheta[m_parameters]), 10e5, 10e5);
 
 
  ATH_MSG_VERBOSE("Extrapolation to Destination Surface: " << destinationSurface);
 
  // the destination parameters
  const Trk::TrackParameters* destParameters = nullptr;
  // the standard validation ...
  if (!m_materialCollectionValidation && !m_direct)
    destParameters = m_extrapolator->extrapolate(ctx,
                                                 startParameters,
                                                 destinationSurface,
                                                 Trk::alongMomentum,
                                                 false,
                                                 static_cast<Trk::ParticleHypothesis>(m_particleType.value()),
                                                 Trk::addNoise).release();
  else if (!m_direct) { // material collection validation
    // get the vector of TrackStateOnSurfaces back
    const std::vector<const Trk::TrackStateOnSurface*>*
      collectedMaterial = m_extrapolator->extrapolateM(ctx,
                                                       startParameters,
                                                       destinationSurface,
                                                       Trk::alongMomentum,
                                                       false,
                                                       static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
 
    // get the last one and clone it
    if (collectedMaterial && !collectedMaterial->empty()) {
      // get the last track state on surface & clone the destination parameters
      const Trk::TrackStateOnSurface* destinationState = collectedMaterial->back();
      destParameters = destinationState->trackParameters() ? destinationState->trackParameters()->clone() : nullptr;
      m_collectedLayerFront += collectedMaterial->size();
      // delete the layers / cleanup
      for (const auto* tsos : *collectedMaterial) {
        delete tsos;
      }
    }
  } else {
    destParameters = m_extrapolator->extrapolateDirectly(ctx,
                                                         startParameters,
                                                         destinationSurface,
                                                         Trk::alongMomentum,
                                                         false,
                                                         static_cast<Trk::ParticleHypothesis>(m_particleType.value())).
                      release();
  }
  // ----------------------- check if forward call was successful and continue then
  if (destParameters) {
    // successful tries
    ++m_triesBack;
    // record the destination parameters
    ++m_parameters;
    m_parameterLoc1[m_parameters] = destParameters->parameters()[Trk::loc1];
    m_parameterLoc2[m_parameters] = destParameters->parameters()[Trk::loc2];
    m_parameterPhi[m_parameters] = destParameters->parameters()[Trk::phi];
    m_parameterEta[m_parameters] = destParameters->eta();
    m_parameterTheta[m_parameters] = destParameters->parameters()[Trk::theta];
    m_parameterQoverP[m_parameters] = destParameters->parameters()[Trk::qOverP];
    if (destParameters->covariance()) {
      m_covarianceLoc1[m_parameters] = (*destParameters->covariance())(0, 0);
      m_covarianceLoc2[m_parameters] = (*destParameters->covariance())(1, 1);
      m_covariancePhi[m_parameters] = (*destParameters->covariance())(2, 2);
      m_covarianceTheta[m_parameters] = (*destParameters->covariance())(3, 3);
      m_covarianceQoverP[m_parameters] = (*destParameters->covariance())(4, 4);
      m_covarianceDeterminant[m_parameters] = (destParameters->covariance())->determinant();
    } else {
      m_covarianceLoc1[m_parameters] = 0;
      m_covarianceLoc2[m_parameters] = 0;
      m_covariancePhi[m_parameters] = 0;
      m_covarianceTheta[m_parameters] = 0;
      m_covarianceQoverP[m_parameters] = 0;
      m_covarianceDeterminant[m_parameters] = 0;
    }
    // record the destination parameters
    const Amg::Vector3D& destinationPosition = destParameters->position();
    m_destinationX = destinationPosition.x();
    m_destinationY = destinationPosition.y();
    m_destinationZ = destinationPosition.z();
    m_destinationR = destinationPosition.perp();
 
    // now simply go backwards
    const Trk::TrackParameters* backParameters = nullptr;
    // the standard validation ...
    if (!m_materialCollectionValidation && !m_direct)
      backParameters = m_extrapolator->extrapolate(ctx,
                                                   *destParameters,
                                                   startSurface,
                                                   Trk::oppositeMomentum,
                                                   false,
                                                   static_cast<Trk::ParticleHypothesis>(m_particleType.value()),
                                                   Trk::removeNoise).release();
    else if (!m_direct) { // material collection validation
      // get the vector of TrackStateOnSurfaces back
      const std::vector<const Trk::TrackStateOnSurface*>*
        collectedBackMaterial = m_extrapolator->extrapolateM(ctx,
                                                             *destParameters,
                                                             startSurface,
                                                             Trk::oppositeMomentum,
                                                             false,
                                                             static_cast<Trk::ParticleHypothesis>(m_particleType.value()));
      // get the last one and clone it
      if (collectedBackMaterial && !collectedBackMaterial->empty()) {
        // get the last track state on surface & clone the destination parameters
        const Trk::TrackStateOnSurface* startState = collectedBackMaterial->back();
        // assign the last ones of the call
        backParameters = startState->trackParameters() ? startState->trackParameters()->clone() : nullptr;
        m_collectedLayerBack += collectedBackMaterial->size();
        // delete the layers / cleanup
        for (const auto* tsos : *collectedBackMaterial) {
          delete tsos;
        }
      }
    } else {
      backParameters = m_extrapolator->extrapolateDirectly(ctx,
                                                           *destParameters,
                                                           startSurface,
                                                           Trk::oppositeMomentum,
                                                           false,
                                                           static_cast<Trk::ParticleHypothesis>(m_particleType.value()))
                        .release();
    }
    // ----------------------- check if backward call was successful and continue then
    if (backParameters) {
      ATH_MSG_VERBOSE("Back Parameters : " << *backParameters);
 
      // record the back extrapolated ones
      ++m_parameters;
      m_parameterLoc1[m_parameters] = backParameters->parameters()[Trk::loc1];
      m_parameterLoc2[m_parameters] = backParameters->parameters()[Trk::loc2];
      m_parameterPhi[m_parameters] = backParameters->parameters()[Trk::phi];
      m_parameterEta[m_parameters] = backParameters->eta();
      m_parameterTheta[m_parameters] = backParameters->parameters()[Trk::theta];
      m_parameterQoverP[m_parameters] = backParameters->parameters()[Trk::qOverP];
      if (backParameters->covariance()) {
        m_covarianceLoc1[m_parameters] = (*backParameters->covariance())(0, 0);
        m_covarianceLoc2[m_parameters] = (*backParameters->covariance())(1, 1);
        m_covariancePhi[m_parameters] = (*backParameters->covariance())(2, 2);
        m_covarianceTheta[m_parameters] = (*backParameters->covariance())(3, 3);
        m_covarianceQoverP[m_parameters] = (*backParameters->covariance())(4, 4);
        m_covarianceDeterminant[m_parameters] = (backParameters->covariance())->determinant();
      } else {
        m_covarianceLoc1[m_parameters] = 0;
        m_covarianceLoc2[m_parameters] = 0;
        m_covariancePhi[m_parameters] = 0;
        m_covarianceTheta[m_parameters] = 0;
        m_covarianceQoverP[m_parameters] = 0;
        m_covarianceDeterminant[m_parameters] = 0;
      }
      // memory cleanup
      delete backParameters;
    } else ++m_breaksBack;
    // memory cleanup
    delete destParameters;
  } else ++m_breaksFront;
  // increase ones more
  ++m_parameters;
  // memory cleanup
 
  if (m_validationTree) m_validationTree->Fill();
 
 
  //std::cout<<"Cleaning up..."<<std::endl;
  //delete covariance;
 
  return StatusCode::SUCCESS;
}
 
//============================================================================================
Amg::Transform3D
Trk::ExtrapolationValidation::createTransform(double x, double y, double z, double phi, double theta, double alphaZ) {
  if (phi != 0. && theta != 0.) {
    // create the Start Surface
    Amg::Vector3D surfacePosition(x, y, z);
    // z direction
    Amg::Vector3D surfaceZdirection(cos(phi) * sin(theta),
                                    sin(phi) * sin(theta),
                                    cos(theta));
    // the global z axis
    Amg::Vector3D zAxis(0., 0., 1.);
    // the y direction
    Amg::Vector3D surfaceYdirection(zAxis.cross(surfaceZdirection));
    // the x direction
    Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
 
    double nx = 1. / sqrt(
      surfaceXdirection[0] * surfaceXdirection[0] + surfaceXdirection[1] * surfaceXdirection[1] + surfaceXdirection[2] *
      surfaceXdirection[2]);
    double ny = 1. / sqrt(
      surfaceYdirection[0] * surfaceYdirection[0] + surfaceYdirection[1] * surfaceYdirection[1] + surfaceYdirection[2] *
      surfaceYdirection[2]);
    surfaceXdirection[0] *= nx;
    surfaceXdirection[1] *= nx;
    surfaceXdirection[2] *= nx;
 
    surfaceYdirection[0] *= ny;
    surfaceYdirection[1] *= ny;
    surfaceYdirection[2] *= ny;
    // the rotation
    Amg::RotationMatrix3D surfaceRotation;
    surfaceRotation.col(0) = surfaceXdirection;
    surfaceRotation.col(1) = surfaceYdirection;
    surfaceRotation.col(2) = surfaceZdirection;
    // return it
    if (alphaZ == 0.) return Amg::Transform3D(surfaceRotation, surfacePosition);
 
    Amg::Transform3D nominalTransform(surfaceRotation, surfacePosition);
    return Amg::Transform3D(nominalTransform * Amg::AngleAxis3D(alphaZ, zAxis));
  }
 
  return Amg::Transform3D(Amg::Translation3D(x, y, z));
}