|
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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
#include "DQM/SiStripCommissioningSources/plugins/tracking/SiStripFineDelayTOF.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <DataFormats/TrackReco/interface/Track.h>
#include <iostream>
#include "TVector3.h"
#include "TMath.h"
//#define TIF_COSMIC_SETUP
double SiStripFineDelayTOF::timeOfFlight(
bool cosmics, bool field, double* trackParameters, double* hit, double* phit, bool onDisk) {
// case 1: cosmics with no field.
if (cosmics && !field) {
return timeOfFlightCosmic(hit, phit);
}
// case 2: cosmics with field.
else if (cosmics && field) {
return timeOfFlightCosmicB(trackParameters, hit, phit, onDisk);
}
// case 3: beam with no field.
else if (!cosmics && !field) {
return timeOfFlightBeam(hit, phit);
}
// case 4: beam with field.
else {
return timeOfFlightBeamB(trackParameters, hit, phit, onDisk);
}
}
double SiStripFineDelayTOF::timeOfFlightCosmic(double* hit, double* phit) {
// constants
const double c = 30; // cm/ns
#ifndef TIF_COSMIC_SETUP
const double Rmu = 385; // cm
const double zmu = 560; // cm
// estimate the time for crossing the barrel
TVector3 r0(hit[0], hit[1], hit[2]);
TVector3 pr(phit[0], phit[1], phit[2]);
TVector2 r0_xy = r0.XYvector();
TVector2 pr_xy = pr.Unit().XYvector();
double t_barrel =
((r0_xy * pr_xy) + sqrt((r0_xy * pr_xy) * (r0_xy * pr_xy) - r0_xy.Mod2() + Rmu * Rmu)) / c * pr.Mag() / pr.Perp();
// estimate the time for crossing endcaps
double t_endcap = fabs(((phit[2] / fabs(phit[2]) * zmu) + hit[2]) / c * pr.Mag() / pr.Pz());
// take the minimum
return t_barrel < t_endcap ? t_barrel : t_endcap;
#else
const double y_trigger = 100; //cm
double p = sqrt(phit[0] * phit[0] + phit[1] * phit[1] + phit[2] * phit[2]);
// LogDebug("TOF") << "momentum:" << phit[0] << " " << phit[1] << " " << phit[2];
// LogDebug("TOF") << "p/py=" << p/phit[1] << " Y0,Y,dY = " << y_trigger << " " << hit[1] << " " << (y_trigger-hit[1]);
// LogDebug("TOF") << "d, t=d/c : " << ((y_trigger-hit[1])*(p/phit[1])) << " " << ((y_trigger-hit[1])*(p/phit[1])/c);
return fabs((y_trigger - hit[1]) * (p / phit[1]) / c);
#endif
}
double SiStripFineDelayTOF::timeOfFlightCosmicB(double* trackParameters, double* hit, double* phit, bool onDisk) {
// constants
const double Rmu = 385; // cm
const double zmu = 560; // cm
const double c = 30; // cm/ns
// track parameters
double& kappa = trackParameters[0];
double& theta = trackParameters[1];
double& phi_0 = trackParameters[2];
double& d_0 = trackParameters[3];
double& z_0 = trackParameters[4];
double invkappa = 1 / kappa;
// computes the value of the track parameter that correspond to the hit, relative to phi_0
//double phi = kappa*tan(theta)*(hit[2]-z_0);
double phi = getPhi(trackParameters, hit, onDisk) - phi_0;
// computes the value of the track parameter that correspond to the muon system, relative to phi_0
// phi_mu = phi_mu0 - phi_0
double phi_mu_b = (kappa > 0 ? -1 : 1) * acos((Rmu * Rmu - d_0 * d_0 - 2 * invkappa * invkappa + 2 * invkappa * d_0) /
(2 * invkappa * d_0 - 2 * invkappa * invkappa));
double phi_mu_e = kappa * tan(theta) * ((phit[2] < 0 ? 1 : -1) * zmu - z_0);
// estimate the time for crossing the barrel
double t_barrel = fabs(invkappa / (c * sin(theta)) * (phi - phi_mu_b));
// estimate the time for crossing endcaps
double t_endcap = fabs(invkappa / (c * sin(theta)) * (phi - phi_mu_e));
// take the minimum
return t_barrel < t_endcap ? t_barrel : t_endcap;
}
double SiStripFineDelayTOF::timeOfFlightBeam(double* hit, double*) {
// constants
const double c = 30; // cm/ns
TVector3 r0(hit[0], hit[1], hit[2]);
return r0.Mag() / c;
}
double SiStripFineDelayTOF::timeOfFlightBeamB(double* trackParameters, double* hit, double* phit, bool onDisk) {
// constants
const double c = 30; // cm/ns
// track parameters
double& theta = trackParameters[1];
// returns the time of flight from the origin
return onDisk ? fabs(hit[2] / (c * cos(theta))) : fabs(sqrt(hit[0] * hit[0] + hit[1] * hit[1]) / (c * sin(theta)));
}
double SiStripFineDelayTOF::x(double* trackParameters, double phi) {
return trackParameters[3] * sin(trackParameters[2]) + (1 / trackParameters[0]) * (sin(phi) - sin(trackParameters[2]));
}
double SiStripFineDelayTOF::y(double* trackParameters, double phi) {
return -trackParameters[3] * cos(trackParameters[2]) -
(1 / trackParameters[0]) * (cos(phi) - cos(trackParameters[2]));
}
double SiStripFineDelayTOF::z(double* trackParameters, double phi) {
return trackParameters[4] + (phi - trackParameters[2]) / (trackParameters[0] * tan(trackParameters[1]));
}
double SiStripFineDelayTOF::getPhi(double* trackParameters, double* hit, bool onDisk) {
if (onDisk) //use z coordinate to find phi
return trackParameters[2] + trackParameters[0] * tan(trackParameters[1]) * (hit[2] - trackParameters[4]);
else // use x,y coordinate to find phi
{
double phi = 0;
if (trackParameters[2] > -2.35 && trackParameters[2] < -0.78) {
// first use y
phi = acos(-trackParameters[0] *
(hit[1] - (trackParameters[3] - 1 / trackParameters[0]) * (-cos(trackParameters[2]))));
// use x to resolve the ambiguity
if (fabs(x(trackParameters, phi) - hit[0]) > fabs(x(trackParameters, -phi) - hit[0]))
phi *= -1.;
} else {
// first use x
phi =
asin(trackParameters[0] * (hit[0] - (trackParameters[3] - 1 / trackParameters[0]) * sin(trackParameters[2])));
// use y to resolve the ambiguity
if ((fabs(y(trackParameters, phi) - hit[1]) > fabs(y(trackParameters, TMath::Pi() - phi) - hit[1])))
phi = phi > 0 ? TMath::Pi() - phi : -TMath::Pi() - phi;
}
return phi;
}
return 0.;
}
void SiStripFineDelayTOF::trackParameters(const reco::Track& tk, double* trackParameters) {
math::XYZPoint position = tk.outerPosition();
const math::XYZVector& momentum = tk.outerMomentum();
LogDebug("trackParameters") << "outer position: " << position.x() << " " << position.y() << " " << position.z();
LogDebug("trackParameters") << "outer momentum: " << momentum.x() << " " << momentum.y() << " " << momentum.z();
math::XYZVector fieldDirection(0, 0, 1);
// computes the center of curvature
math::XYZVector radius = momentum.Cross(fieldDirection).Unit();
position -= radius / trackParameters[0];
LogDebug("trackParameters") << "center of curvature: " << position.x() << " " << position.y() << " " << position.z();
// the transverse IP
trackParameters[3] = position.rho() - fabs(1 / trackParameters[0]);
if (trackParameters[0] > 0)
trackParameters[3] *= -1.;
// phi_0
double phi_out = trackParameters[2];
double phi_0 = position.phi() - TMath::Pi() / 2;
if (trackParameters[0] < 0)
phi_0 -= TMath::Pi();
if (phi_0 < -2 * TMath::Pi())
phi_0 += 2 * TMath::Pi();
if (phi_0 > 2 * TMath::Pi())
phi_0 -= 2 * TMath::Pi();
if (phi_0 > 0)
phi_0 -= 2 * TMath::Pi();
LogDebug("trackParameters") << "phi_0: " << phi_0;
trackParameters[2] = phi_0;
// z_0
trackParameters[4] = tk.outerPosition().z() - (phi_out - phi_0) / TMath::Tan(trackParameters[1]) / trackParameters[0];
LogDebug("trackParameters") << "z_0: " << tk.outerPosition().z() << " " << trackParameters[4] << " "
<< tk.innerPosition().z();
}
|