fastsim::HelixTrajectory Class Reference

Mathematical representation of a helix. More...

#include <HelixTrajectory.h>

Inheritance diagram for fastsim::HelixTrajectory:

Public Member Functions
bool	crosses (const BarrelSimplifiedGeometry &layer) const override
	Check if an intersection of the trajectory with a barrel layer exists. More...
double	getRadParticle (double phi) const
	Return distance of particle from center of the detector if it was at given angle phi of the helix. More...
	HelixTrajectory (const Particle &particle, double magneticFieldZ)
	Constructor. More...
void	move (double deltaTimeC) override
	Move the particle along the helix trajectory for a given time. More...
double	nextCrossingTimeC (const BarrelSimplifiedGeometry &layer, bool onLayer=false) const override
	Return delta time (t*c) of the next intersection of trajectory and barrel layer. More...
Public Member Functions inherited from fastsim::Trajectory
const math::XYZTLorentzVector &	getMomentum ()
	Simple getter: return momentum of the particle that was used to create trajectory. More...
const math::XYZTLorentzVector &	getPosition ()
	Simple getter: return position of the particle that was used to create trajectory. More...
double	nextCrossingTimeC (const SimplifiedGeometry &layer, bool onLayer=false) const
	Return delta time (t*c) of the next intersection of trajectory and generic layer. More...
double	nextCrossingTimeC (const ForwardSimplifiedGeometry &layer, bool onLayer=false) const
	Return delta time (t*c) of the next intersection of trajectory and forward layer. More...
virtual	~Trajectory ()

Private Attributes
const double	centerR_
	Distance of the center of the helix from the center of the tracker. More...
const double	centerX_
	X-coordinate of the center of the helix. More...
const double	centerY_
	Y-coordinate of the center of the helix. More...
const double	maxR_
	The maximum distance of the helix from the center of the tracker. More...
const double	minR_
	The minimal distance of the helix from the center of the tracker. More...
const double	phi_
	Ranges from 0 to 2PI: 0 corresponds to the positive X direction, phi increases counterclockwise. More...
const double	phiSpeed_
	The angular speed of the particle on the helix trajectory. More...
const double	radius_
	The radius of the helix. More...

Additional Inherited Members
Static Public Member Functions inherited from fastsim::Trajectory
static std::unique_ptr< Trajectory >	createTrajectory (const fastsim::Particle &particle, const double magneticFieldZ)
	Calls constructor of derived classes. More...
Protected Member Functions inherited from fastsim::Trajectory
	Trajectory (const fastsim::Particle &particle)
	Constructor of base class. More...
Protected Attributes inherited from fastsim::Trajectory
math::XYZTLorentzVector	momentum_
	momentum of the particle that was used to create trajectory More...
math::XYZTLorentzVector	position_
	position of the particle that was used to create trajectory More...

Detailed Description

Mathematical representation of a helix.

Reflects trajectory of a charged particle in a magnetic field. The trajectory is defined by cylindrical coordinates (see definition of variables for more information).

Definition at line 18 of file HelixTrajectory.h.

Constructor & Destructor Documentation

fastsim::HelixTrajectory::HelixTrajectory	(	const Particle &	particle,
		double	magneticFieldZ
	)

Constructor.

The magnetic field is (to good approximation) constant between two tracker layers (and only in Z-direction).

Parameters

particle	A (usually charged) particle.
magneticFieldZ	The magnetic field.

Definition at line 15 of file HelixTrajectory.cc.

    : Trajectory(particle)
      // exact: r = gamma*beta*m_0*c / (q*e*B) = p_T / (q * e * B)
      // momentum in units of GeV/c: r = p_T * 10^9 / (c * q * B)
      // in cmssw units: r = p_T / (c * 10^-4 * q * B)
      ,
      radius_(
          std::abs(momentum_.Pt() / (fastsim::Constants::speedOfLight * 1e-4 * particle.charge() * magneticFieldZ))),
      phi_(std::atan(momentum_.Py() / momentum_.Px()) +
           (momentum_.Px() * particle.charge() < 0 ? 3. * M_PI / 2. : M_PI / 2.))
      // maybe consider (for -pi/2<x<pi/2)
      // cos(atan(x)) = 1 / sqrt(x^2+1)
      // -> cos(atan(x) + pi/2)  = - x / sqrt(x^2+1)
      // -> cos(atan(x) +3*pi/2) = + x / sqrt(x^2+1)
      // sin(atan(x)) = x / sqrt(x^2+1)
      // -> sin(atan(x) + pi/2)  = + 1 / sqrt(x^2+1)
      // -> sin(atan(x) +3*pi/2) = - 1 / sqrt(x^2+1)
      ,
      centerX_(position_.X() -
               radius_ * (momentum_.Py() / momentum_.Px()) /
                   std::sqrt((momentum_.Py() / momentum_.Px()) * (momentum_.Py() / momentum_.Px()) + 1) *
                   (momentum_.Px() * particle.charge() < 0 ? 1. : -1.)),
      centerY_(position_.Y() -
               radius_ * 1 / std::sqrt((momentum_.Py() / momentum_.Px()) * (momentum_.Py() / momentum_.Px()) + 1) *
                   (momentum_.Px() * particle.charge() < 0 ? -1. : 1.))
      //, centerX_(position_.X() - radius_*std::cos(phi_))
      //, centerY_(position_.Y() - radius_*std::sin(phi_))
      ,
      centerR_(std::sqrt(centerX_ * centerX_ + centerY_ * centerY_)),
      minR_(std::abs(centerR_ - radius_)),
      maxR_(centerR_ + radius_)
      // omega = q * e * B / (gamma * m) = q * e *B / (E / c^2) = q * e * B * c^2 / E
      // omega: negative for negative q -> seems to be what we want.
      // energy in units of GeV: omega = q * B * c^2 / (E * 10^9)
      // in cmssw units: omega[1/ns] = q * B * c^2 * 10^-4 / E
      ,
      phiSpeed_(-particle.charge() * magneticFieldZ * fastsim::Constants::speedOfLight *
                fastsim::Constants::speedOfLight * 1e-4 / momentum_.E()) {
  ;
}

Member Function Documentation

bool fastsim::HelixTrajectory::crosses ( const BarrelSimplifiedGeometry &  layer ) const

overridevirtual

Check if an intersection of the trajectory with a barrel layer exists.

Parameters

layer

A barrel layer.

Implements fastsim::Trajectory.

Definition at line 56 of file HelixTrajectory.cc.

References fastsim::BarrelSimplifiedGeometry::getRadius(), maxR_, and minR_.

Referenced by nextCrossingTimeC().

                                                                                {
  return (minR_ < layer.getRadius() && maxR_ > layer.getRadius());
}

double fastsim::HelixTrajectory::getRadParticle

(

double

phi

)

const

Return distance of particle from center of the detector if it was at given angle phi of the helix.

Parameters

phi	angle of the helix

Definition at line 209 of file HelixTrajectory.cc.

References centerX_, centerY_, funct::cos(), radius_, funct::sin(), and mathSSE::sqrt().

Referenced by nextCrossingTimeC().

                                                              {
  return sqrt((centerX_ + radius_ * std::cos(phi)) * (centerX_ + radius_ * std::cos(phi)) +
              (centerY_ + radius_ * std::sin(phi)) * (centerY_ + radius_ * std::sin(phi)));
}

void fastsim::HelixTrajectory::move ( double  deltaTimeC )

overridevirtual

Move the particle along the helix trajectory for a given time.

Parameters

deltaTimeC

Time in units of t*c..

Implements fastsim::Trajectory.

Definition at line 193 of file HelixTrajectory.cc.

References centerX_, centerY_, funct::cos(), SiPixelRawToDigiRegional_cfi::deltaPhi, fastsim::Trajectory::momentum_, phi_, phiSpeed_, fastsim::Trajectory::position_, radius_, funct::sin(), and fastsim::Constants::speedOfLight.

Referenced by Vispa.Gui.WidgetContainer.WidgetContainer::autosize(), Vispa.Gui.VispaWidget.VispaWidget::dragWidget(), Vispa.Gui.VispaWidget.VispaWidget::setZoom(), and Vispa.Gui.PortConnection.PointToPointConnection::updateConnection().

                                                   {
  double deltaT = deltaTimeC / fastsim::Constants::speedOfLight;
  double deltaPhi = phiSpeed_ * deltaT;
  position_.SetXYZT(centerX_ + radius_ * std::cos(phi_ + deltaPhi),
                    centerY_ + radius_ * std::sin(phi_ + deltaPhi),
                    position_.Z() + momentum_.Z() / momentum_.E() * deltaTimeC,
                    position_.T() + deltaT);
  // Rotation defined by
  // x' = x cos θ - y sin θ
  // y' = x sin θ + y cos θ
  momentum_.SetXYZT(momentum_.X() * std::cos(deltaPhi) - momentum_.Y() * std::sin(deltaPhi),
                    momentum_.X() * std::sin(deltaPhi) + momentum_.Y() * std::cos(deltaPhi),
                    momentum_.Z(),
                    momentum_.E());
}

double fastsim::HelixTrajectory::nextCrossingTimeC ( const BarrelSimplifiedGeometry &  layer, bool  onLayer = false  ) const

overridevirtual

Return delta time (t*c) of the next intersection of trajectory and barrel layer.

This function solves the quadratic equation (basically intersection of two circles with a given radius) in order to calculate the moment in time when the particle's trajectory intersects with a given barrel layer.

Parameters

layer	A barrel layer.
onLayer	Specify if the particle already is on the layer (in that case the second solution has to be picked).

Returns

t*c [ns * cm/ns] of next intersection (-1 if there is none).

Implements fastsim::Trajectory.

Definition at line 60 of file HelixTrajectory.cc.

References a, funct::abs(), b, HltBtagPostValidation_cff::c, centerX_, centerY_, funct::cos(), crosses(), dumpMFGeometry_cfg::delta, MillePedeFileConverter_cfg::e, F(), callgraph::G, fastsim::BarrelSimplifiedGeometry::getRadius(), getRadParticle(), M_PI, min(), fastsim::Trajectory::momentum_, fastsim::StraightTrajectory::nextCrossingTimeC(), phi_, phiSpeed_, fastsim::Trajectory::position_, radius_, funct::sin(), fastsim::Constants::speedOfLight, mathSSE::sqrt(), RandomServiceHelper::t1, and RandomServiceHelper::t2.

                                                                                                          {
  if (!crosses(layer))
    return -1;

  // solve the following equation for sin(phi)
  // (x^2 + y^2 = R_L^2)     (1)      the layer
  // x = x_c + R_H*cos(phi)  (2)      the helix in the xy plane
  // y = y_c + R_H*sin(phi)  (3)      the helix in the xy plane
  // with
  // R_L: the radius of the layer
  // x_c,y_c the center of the helix in xy plane
  // R_H, the radius of the helix
  // phi, the phase of the helix
  //
  // substitute (2) and (3) in (1)
  // =>
  //   x_c^2 + 2*x_c*R_H*cos(phi) + R_H^2*cos^2(phi)
  // + y_c^2 + 2*y_c*R_H*sin(phi) + R_H^2*sin^2(phi)
  // = R_L^2
  // =>
  // (x_c^2 + y_c^2 + R_H^2 - R_L^2) + (2*y_c*R_H)*sin(phi) = -(2*x_c*R_H)*cos(phi)
  //
  // rewrite
  //               E                 +       F    *sin(phi) =      G     *cos(phi)
  // =>
  // E^2 + 2*E*F*sin(phi) + F^2*sin^2(phi) = G^2*(1-sin^2(phi))
  // rearrange
  // sin^2(phi)*(F^2 + G^2) + sin(phi)*(2*E*F) + (E^2 - G^2) = 0
  //
  // rewrite
  // sin^2(phi)*     a      + sin(phi)*   b    +      c      = 0
  // => sin(phi) = (-b +/- sqrt(b^2 - 4*ac)) / (2*a)
  // with
  // a = F^2 + G^2
  // b = 2*E*F
  // c = E^2 - G^2

  double E = centerX_ * centerX_ + centerY_ * centerY_ + radius_ * radius_ - layer.getRadius() * layer.getRadius();
  double F = 2 * centerY_ * radius_;
  double G = 2 * centerX_ * radius_;

  double a = F * F + G * G;
  double b = 2 * E * F;
  double c = E * E - G * G;

  double delta = b * b - 4 * a * c;

  // case of no solution
  if (delta < 0) {
    // Should not be reached: Full Propagation does always have a solution "if(crosses(layer)) == -1"
    // Even if particle is outside all layers -> can turn around in magnetic field
    return -1;
  }

  // Uses a numerically more stable procedure:
  // https://people.csail.mit.edu/bkph/articles/Quadratics.pdf
  double sqrtDelta = sqrt(delta);
  double phi1 = 0, phi2 = 0;
  if (b < 0) {
    phi1 = std::asin((2. * c) / (-b + sqrtDelta));
    phi2 = std::asin((-b + sqrtDelta) / (2. * a));
  } else {
    phi1 = std::asin((-b - sqrtDelta) / (2. * a));
    phi2 = std::asin((2. * c) / (-b - sqrtDelta));
  }

  // asin is ambiguous, make sure to have the right solution
  if (std::abs(layer.getRadius() - getRadParticle(phi1)) > 1.0e-2) {
    phi1 = -phi1 + M_PI;
  }
  if (std::abs(layer.getRadius() - getRadParticle(phi2)) > 1.0e-2) {
    phi2 = -phi2 + M_PI;
  }

  // another ambiguity
  if (phi1 < 0) {
    phi1 += 2. * M_PI;
  }
  if (phi2 < 0) {
    phi2 += 2. * M_PI;
  }

  // find the corresponding times when the intersection occurs
  // make sure they are positive
  double t1 = (phi1 - phi_) / phiSpeed_;
  while (t1 < 0) {
    t1 += 2 * M_PI / std::abs(phiSpeed_);
  }
  double t2 = (phi2 - phi_) / phiSpeed_;
  while (t2 < 0) {
    t2 += 2 * M_PI / std::abs(phiSpeed_);
  }

  // Check if propagation successful (numerical reasons): both solutions (phi1, phi2) have to be on the layer (same radius)
  // Can happen due to numerical instabilities of geometrical function (if momentum is almost parallel to x/y axis)
  // Get crossingTimeC from StraightTrajectory as good approximation
  if (std::abs(layer.getRadius() - getRadParticle(phi1)) > 1.0e-2 ||
      std::abs(layer.getRadius() - getRadParticle(phi2)) > 1.0e-2) {
    StraightTrajectory traj(*this);
    return traj.nextCrossingTimeC(layer, onLayer);
  }

  // if the particle is already on the layer, we need to make sure the 2nd solution is picked.
  // happens if particle turns around in the magnetic field instead of hitting the next layer
  if (onLayer) {
    bool particleMovesInwards = momentum_.X() * position_.X() + momentum_.Y() * position_.Y() < 0;

    double posX1 = centerX_ + radius_ * std::cos(phi1);
    double posY1 = centerY_ + radius_ * std::sin(phi1);
    double momX1 = momentum_.X() * std::cos(phi1 - phi_) - momentum_.Y() * std::sin(phi1 - phi_);
    double momY1 = momentum_.X() * std::sin(phi1 - phi_) + momentum_.Y() * std::cos(phi1 - phi_);
    bool particleMovesInwards1 = momX1 * posX1 + momY1 * posY1 < 0;

    double posX2 = centerX_ + radius_ * std::cos(phi2);
    double posY2 = centerY_ + radius_ * std::sin(phi2);
    double momX2 = momentum_.X() * std::cos(phi2 - phi_) - momentum_.Y() * std::sin(phi2 - phi_);
    double momY2 = momentum_.X() * std::sin(phi2 - phi_) + momentum_.Y() * std::cos(phi2 - phi_);
    bool particleMovesInwards2 = momX2 * posX2 + momY2 * posY2 < 0;

    if (particleMovesInwards1 != particleMovesInwards) {
      return t1 * fastsim::Constants::speedOfLight;
    } else if (particleMovesInwards2 != particleMovesInwards) {
      return t2 * fastsim::Constants::speedOfLight;
    }
    // try to catch numerical issues again..
    else {
      return -1;
    }
  }

  return std::min(t1, t2) * fastsim::Constants::speedOfLight;
}

Member Data Documentation

const double fastsim::HelixTrajectory::centerR_

private

Distance of the center of the helix from the center of the tracker.

Definition at line 62 of file HelixTrajectory.h.

const double fastsim::HelixTrajectory::centerX_

private

X-coordinate of the center of the helix.

Definition at line 60 of file HelixTrajectory.h.

Referenced by getRadParticle(), move(), and nextCrossingTimeC().

const double fastsim::HelixTrajectory::centerY_

private

Y-coordinate of the center of the helix.

Definition at line 61 of file HelixTrajectory.h.

Referenced by getRadParticle(), move(), and nextCrossingTimeC().

const double fastsim::HelixTrajectory::maxR_

private

The maximum distance of the helix from the center of the tracker.

Definition at line 64 of file HelixTrajectory.h.

Referenced by crosses().

const double fastsim::HelixTrajectory::minR_

private

The minimal distance of the helix from the center of the tracker.

Definition at line 63 of file HelixTrajectory.h.

Referenced by crosses().

const double fastsim::HelixTrajectory::phi_

private

Ranges from 0 to 2PI: 0 corresponds to the positive X direction, phi increases counterclockwise.

The angle of the particle alone the helix.

Definition at line 58 of file HelixTrajectory.h.

Referenced by move(), and nextCrossingTimeC().

const double fastsim::HelixTrajectory::phiSpeed_

private

The angular speed of the particle on the helix trajectory.

Definition at line 65 of file HelixTrajectory.h.

Referenced by move(), and nextCrossingTimeC().

const double fastsim::HelixTrajectory::radius_

private

The radius of the helix.

Definition at line 57 of file HelixTrajectory.h.

Referenced by getRadParticle(), move(), and nextCrossingTimeC().