-
Notifications
You must be signed in to change notification settings - Fork 7
/
propagate-toz-pseudocode.txt
101 lines (87 loc) · 4.17 KB
/
propagate-toz-pseudocode.txt
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
Input track: 6 parameters, symmetric 6x6 covariance matrix (so 21 free parameters)
inputtrk={
par = {-12.806846618652344, -7.723824977874756, 38.13014221191406,0.23732035065189902, -2.613372802734375, 0.35594117641448975},
cov = {6.290299552347278e-07,4.1375109560704004e-08,7.526661534029699e-07,2.0973730840978533e-07,1.5431574240665213e-07,9.626245400795597e-08,-2.804026640189443e-06,
6.219111130687595e-06,2.649119409845118e-07,0.00253512163402557,-2.419662877381737e-07,4.3124190760040646e-07,3.1068903991780678e-09,0.000923913115050627,
0.00040678296006807003,-7.755406890332818e-07,1.68539375883925e-06,6.676875566525437e-08,0.0008420574605423793,7.356584799406111e-05,0.0002306247719158348},
q = 1,
hitIdx = {1, 0, 17, 16, 36, 35, 33, 34, 59, 58, 70, 85, 101, 102, 116, 117, 132, 133, 152, 169, 187, 202}
};
Input hit: 3 parameters, symmetric 3x3 covariance matrix (so 6 free parameters)
inputhit = {
par = {-20.7824649810791, -12.24150276184082, 57.8067626953125},
cov = {2.545517190810642e-06,-2.6680759219743777e-06,2.8030024168401724e-06,0.00014160551654640585,0.00012282167153898627,11.385087966918945}
};
Track track, outtrack;
Hit hit;
Randomly smear parameters of the track:
//par
for (size_t ip=0;ip<6;++ip) {
track.par[ip] = (1+smear*randn(0,1))*inputtrk.par[ip];
}
//cov
for (size_t ip=0;ip<21;++ip) {
track.cov[ip] = (1+smear*randn(0,1))*inputtrk.cov[ip];
}
//q
track.q = inputtrk.q-2*ceil(-0.5 + (float)rand() / RAND_MAX);//fixme check
Randomly smear parameters of the hit:
//par
for (size_t ip=0;ip<3;++ip) {
hit.par[ip] = (1+smear*randn(0,1))*inputhit.pos[ip];
}
//cov
for (size_t ip=0;ip<6;++ip) {
hit.cov[ip] = (1+smear*randn(0,1))*inputhit.cov[ip];
}
Propagate a single track (repeat for nevts and bsize*nb tracks per event)
**(this is the region to be timed)**
hitPos = hit.par
inPar = track.par
inChg = track.q
inErr = track.cov
outPar = outtrack.par
outErr = outtrack.cov
int x_id = 0;
int y_id = 1;
int z_id = 2;
int ipt_id = 3;
int phi_id = 4;
int theta_id = 5;
const float zout = hitPos[z_id];
const float k = inChg*100/3.8;
const float deltaZ = zout - inPar[z_id];
const float pt = 1./inPar[ipt_id];
const float cosP = cosf(inPar[phi_id]);
const float sinP = sinf(inPar[phi_id]);
const float cosT = cosf(inPar[theta_id]);
const float sinT = sinf(inPar[theta_id]);
const float pxin = cosP*pt;
const float pyin = sinP*pt;
const float alpha = deltaZ*sinT*inPar[ipt_id]/(cosT*k);
const float sina = sinf(alpha); // this can be approximated;
const float cosa = cosf(alpha); // this can be approximated;
outPar[x_id] = inPar[x_id] + k*(pxin*sina - pyin*(1.-cosa));
outPar[y_id] = inPar[y_id] + k*(pyin*sina + pxin*(1.-cosa));
outPar[z_id] = zout;
outPar[ipt_id] = inPar[ipt_id];
outPar[phi_id] = inPar[phi_id]+alpha;
outPar[theta_id]= inPar[theta_id];
6x6 matrix errorProp, temp;
for (size_t i=0;i<6;++i) errorProp[i][i] = 1.;
errorProp[0][2] = cosP*sinT*(sinP*cosa*sCosPsina-cosa)/cosT;
errorProp[0][3] = cosP*sinT*deltaZ*cosa*(1.-sinP*sCosPsina)/(cosT*inPar[ipt_id])-k*(cosP*sina-sinP*(1.-cCosPsina))/(inPar[ipt_id]*inPar[ipt_id]);
errorProp[0][4] = (k/inPar[ipt_id])*(-sinP*sina+sinP*sinP*sina*sCosPsina-cosP*(1.-cCosPsina));
errorProp[0][5] = cosP*deltaZ*cosa*(1.-sinP*sCosPsina)/(cosT*cosT);
errorProp[1][2] = cosa*sinT*(cosP*cosP*sCosPsina-sinP)/cosT;
errorProp[1][3] = sinT*deltaZ*cosa*(cosP*cosP*sCosPsina+sinP)/(cosT*inPar[ipt_id])-k*(sinP*sina+cosP*(1.-cCosPsina))/(inPar[ipt_id]*inPar[ipt_id]);
errorProp[1][4] = (k/inPar[ipt_id])*(-sinP*(1.-cCosPsina)-sinP*cosP*sina*sCosPsina+cosP*sina);
errorProp[1][5] = deltaZ*cosa*(cosP*cosP*sCosPsina+sinP)/(cosT*cosT);
errorProp[4][2] = -inPar[ipt_id]*sinT/(cosT*k);
errorProp[4][3] = sinT*deltaZ/(cosT*k);
errorProp[4][5] = inPar[ipt_id]*deltaZ/(cosT*cosT*k);
//Matrix multiplications of covariance matrices
temp = errorProp*inErr; //MultHelixPropEndcap
outErr = transpose(errorProp)*temp; //MultHelixPropTranspEndcap
** End timed section **
Final print statements to report average final values of all tracks