|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
|
#include "FastSimulation/SimplifiedGeometryPropagator/interface/StraightTrajectory.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/Constants.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometry.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/BarrelSimplifiedGeometry.h"
#include "FastSimulation/SimplifiedGeometryPropagator/interface/ForwardSimplifiedGeometry.h"
double fastsim::StraightTrajectory::nextCrossingTimeC(const fastsim::BarrelSimplifiedGeometry& layer,
bool onLayer) const {
//
// solve the equation
// (x^2 + y^2 = R^2), (1)
// x = x_0 + v_x*t, (2)
// y = y_0 + v_y*t (3)
// for t
//
// subsitute (2) and (3) in (1):
// => 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
// => t = (-b +/- sqrt(b^2 - ac) ) / a see https://en.wikipedia.org/wiki/Quadratic_formula, divide equation by 2
// with a = v_x^2 + v_y^2;
// with b = x_0*v_x + y_0*v_y
// with c = x_0^2 + y_0^2 - R^2
//
// substitute: v_x = p_x / E * c ( note: extra * c absorbed in p_x units)
// substitute: v_y = p_y / E * c ( note: extra * c absorbed in p_y units)
// => t*c = (-b +/- sqrt(b^2 - ac) ) / a * E
// with a = p_x^2 + p_y^2
// with b = p_x*x_0 + p_y*y_0
// with c = x_0^2 + y_0^2 - R^2
double a = momentum_.Perp2();
double b = (position_.X() * momentum_.X() + position_.Y() * momentum_.Y());
double c = position_.Perp2() - layer.getRadius() * layer.getRadius();
double delta = b * b - a * c;
if (delta < 0) {
return -1;
}
double sqrtDelta = sqrt(delta);
//
// return the earliest solution > 0,
// return -1 if no positive solution exists
// note: 'a' always positive, sqrtDelta always positive
//
double tc1 = (-b - sqrtDelta) / a * momentum_.E();
double tc2 = (-b + sqrtDelta) / a * momentum_.E();
if (onLayer && tc2 > 0) {
bool particleMovesInwards = momentum_.X() * position_.X() + momentum_.Y() * position_.Y() < 0;
double posX2 = position_.X() + momentum_.X() / momentum_.E() * tc2;
double posY2 = position_.Y() + momentum_.Y() / momentum_.E() * tc2;
bool particleMovesInwards2 = momentum_.X() * posX2 + momentum_.Y() * posY2 < 0;
if (particleMovesInwards != particleMovesInwards2) {
return tc2;
}
return -1;
}
if (tc1 > 0) {
return tc1;
} else if (tc2 > 0) {
return tc2;
}
return -1.;
}
void fastsim::StraightTrajectory::move(double deltaTimeC) {
position_.SetXYZT(position_.X() + momentum_.X() / momentum_.E() * deltaTimeC,
position_.Y() + momentum_.Y() / momentum_.E() * deltaTimeC,
position_.Z() + momentum_.Z() / momentum_.E() * deltaTimeC,
position_.T() + deltaTimeC / fastsim::Constants::speedOfLight);
}
|