SimplifiedGeometryPropagator/src/StraightTrajectory.cc

0001 #include "FastSimulation/SimplifiedGeometryPropagator/interface/StraightTrajectory.h"
0002 #include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
0003
0004 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
0005 #include "FastSimulation/SimplifiedGeometryPropagator/interface/BarrelSimplifiedGeometry.h"
0006 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ForwardSimplifiedGeometry.h"
0007
0008 double fastsim::StraightTrajectory::nextCrossingTimeC(const fastsim::BarrelSimplifiedGeometry& layer,
0009                                                       bool onLayer) const {
0010   //
0011   // solve the equation
0012   // (x^2 + y^2 = R^2),   (1)
0013   // x = x_0 + v_x*t,     (2)
0014   // y = y_0 + v_y*t      (3)
0015   // for t
0016   //
0017   // subsitute (2) and (3) in (1):
0018   // => t^2*(v_x^2 + v_y^2) + t*( 2*x_0*v_x + 2*y_0*v_y ) + x_0^2 + y_0^2 - R^2 = 0
0019   // => t = (-b +/- sqrt(b^2 - ac) ) / a       see https://en.wikipedia.org/wiki/Quadratic_formula, divide equation by 2
0020   // with a = v_x^2 + v_y^2;
0021   // with b = x_0*v_x + y_0*v_y
0022   // with c = x_0^2 + y_0^2 - R^2
0023   //
0024   // substitute: v_x = p_x / E * c  ( note: extra * c absorbed in p_x units)
0025   // substitute: v_y = p_y / E * c  ( note: extra * c absorbed in p_y units)
0026   // => t*c = (-b +/- sqrt(b^2 - ac) ) / a * E
0027   // with a = p_x^2 + p_y^2
0028   // with b = p_x*x_0 + p_y*y_0
0029   // with c = x_0^2 + y_0^2 - R^2
0030   double a = momentum_.Perp2();
0031   double b = (position_.X() * momentum_.X() + position_.Y() * momentum_.Y());
0032   double c = position_.Perp2() - layer.getRadius() * layer.getRadius();
0033
0034   double delta = b * b - a * c;
0035   if (delta < 0) {
0036     return -1;
0037   }
0038   double sqrtDelta = sqrt(delta);
0039
0040   //
0041   // return the earliest solution > 0,
0042   // return -1 if no positive solution exists
0043   // note: 'a' always positive, sqrtDelta always positive
0044   //
0045   double tc1 = (-b - sqrtDelta) / a * momentum_.E();
0046   double tc2 = (-b + sqrtDelta) / a * momentum_.E();
0047
0048   if (onLayer && tc2 > 0) {
0049     bool particleMovesInwards = momentum_.X() * position_.X() + momentum_.Y() * position_.Y() < 0;
0050
0051     double posX2 = position_.X() + momentum_.X() / momentum_.E() * tc2;
0052     double posY2 = position_.Y() + momentum_.Y() / momentum_.E() * tc2;
0053     bool particleMovesInwards2 = momentum_.X() * posX2 + momentum_.Y() * posY2 < 0;
0054
0055     if (particleMovesInwards != particleMovesInwards2) {
0056       return tc2;
0057     }
0058
0059     return -1;
0060   }
0061
0062   if (tc1 > 0) {
0063     return tc1;
0064   } else if (tc2 > 0) {
0065     return tc2;
0066   }
0067
0068   return -1.;
0069 }
0070
0071 void fastsim::StraightTrajectory::move(double deltaTimeC) {
0072   position_.SetXYZT(position_.X() + momentum_.X() / momentum_.E() * deltaTimeC,
0073                     position_.Y() + momentum_.Y() / momentum_.E() * deltaTimeC,
0074                     position_.Z() + momentum_.Z() / momentum_.E() * deltaTimeC,
0075                     position_.T() + deltaTimeC / fastsim::Constants::speedOfLight);
0076 }