]>
Commit | Line | Data |
---|---|---|
5a5a1232 | 1 | #ifndef REVE_MCHelixLine_H |
2 | #define REVE_MCHelixLine_H | |
3 | ||
4 | #include <Reve/Track.h> | |
5 | #include <vector> | |
6 | ||
7 | namespace Reve { | |
8 | ||
bbfaa1c4 | 9 | struct MCVertex |
10 | { | |
11 | Float_t x,y,z,t; | |
12 | }; | |
5a5a1232 | 13 | |
14 | ||
bbfaa1c4 | 15 | struct MCStruct |
16 | { | |
17 | MCStruct(const MCStruct&); // Not implemented | |
18 | MCStruct& operator=(const MCStruct&); // Not implemented | |
19 | ||
20 | TrackRnrStyle* fRnrMod; | |
21 | std::vector<MCVertex>* fPoints; | |
22 | MCVertex fV; | |
23 | Float_t fVelocity; // size of particle velocity | |
5a5a1232 | 24 | |
bbfaa1c4 | 25 | MCStruct(TrackRnrStyle* rs, MCVertex* v0 , Float_t vel, std::vector<MCVertex>* tpv) : |
26 | fRnrMod (rs), | |
27 | fPoints (tpv), | |
28 | fV (*v0), | |
29 | fVelocity (vel) | |
b97c26da | 30 | { |
bbfaa1c4 | 31 | fPoints->push_back(fV); |
32 | } | |
33 | virtual ~MCStruct() {} | |
34 | }; | |
5a5a1232 | 35 | |
bbfaa1c4 | 36 | /**************************************************************************/ |
37 | // HELIX | |
38 | /**************************************************************************/ | |
5a5a1232 | 39 | |
bbfaa1c4 | 40 | struct MCHelix : public MCStruct |
41 | { | |
42 | // constant | |
43 | Float_t fA; // contains charge and magnetic field data | |
5a5a1232 | 44 | |
bbfaa1c4 | 45 | //parameters dependend pT and pZ size, set in init function |
46 | Float_t fLam; // momentum ratio pT/pZ | |
47 | Float_t fR; // a/pT | |
48 | Float_t fPhiStep; // step size in xy projection, dependent of RnrMode and momentum | |
49 | Float_t fTimeStep; | |
50 | ||
51 | Int_t fN; // step number in helix; | |
52 | Int_t NMax; // max number of points in helix | |
53 | Float_t x_off, y_off; // offset for fitting daughters | |
54 | Float_t sin, cos; | |
55 | Bool_t crosR; | |
5a5a1232 | 56 | |
bbfaa1c4 | 57 | MCHelix(TrackRnrStyle* rs, MCVertex* v0, Float_t vel, |
58 | std::vector<MCVertex>* tpv, Float_t a) : | |
59 | MCStruct(rs, v0 , vel, tpv), | |
60 | fA (a), | |
61 | fLam (0), fR (0), fPhiStep (0), fTimeStep (0), | |
62 | fN (0), NMax (0), | |
63 | x_off(0), y_off(0), | |
64 | sin (0), cos (0), crosR (0) | |
65 | {} | |
5a5a1232 | 66 | |
bbfaa1c4 | 67 | void Init(Float_t pT, Float_t pZ) |
68 | { | |
69 | fN = 0; | |
70 | crosR = false; | |
71 | x_off = 0; y_off = 0; | |
72 | fLam = pZ/pT; | |
73 | fR = pT/fA; | |
5a5a1232 | 74 | |
bbfaa1c4 | 75 | fPhiStep = fRnrMod->fMinAng *TMath::Pi()/180; |
76 | if(fRnrMod->fDelta < TMath::Abs(fR)){ | |
77 | Float_t ang = 2*TMath::ACos(1 - fRnrMod->fDelta/TMath::Abs(fR)); | |
78 | if (ang < fPhiStep) fPhiStep = ang; | |
5a5a1232 | 79 | } |
bbfaa1c4 | 80 | if(fA<0) fPhiStep = -fPhiStep; |
5a5a1232 | 81 | |
bbfaa1c4 | 82 | // printf("MCHelix::init (%f/%f) labda %f time step %e phi step %f \n", pT, pZ,fLam, fTimeStep,fPhiStep); |
83 | fTimeStep = TMath::Abs(fR*fPhiStep)*TMath::Sqrt(1+fLam*fLam)/fVelocity; | |
84 | fTimeStep *= 0.01; //cm->m | |
85 | ||
86 | sin = TMath::Sin(fPhiStep); | |
87 | cos = TMath::Cos(fPhiStep); | |
88 | } | |
89 | ||
90 | void SetBounds() | |
91 | { | |
92 | // check steps for max orbits | |
93 | NMax = Int_t(fRnrMod->fMaxOrbs*TMath::TwoPi()/TMath::Abs(fPhiStep)); | |
94 | // check steps for Z boundaries | |
95 | Float_t nz; | |
96 | if(fLam > 0) { | |
97 | nz = (fRnrMod->fMaxZ - fV.z)/(fLam*TMath::Abs(fR*fPhiStep)); | |
98 | } else { | |
99 | nz = (-fRnrMod->fMaxZ - fV.z)/(fLam*TMath::Abs(fR*fPhiStep)); | |
100 | } | |
101 | // printf("steps in helix line %d nz %f vz %f\n", NMax, nz, fV.z); | |
102 | if (nz < NMax) NMax = Int_t(nz); | |
5a5a1232 | 103 | |
bbfaa1c4 | 104 | // check steps if circles intersect |
105 | if(TMath::Sqrt(fV.x*fV.x+fV.y*fV.y) < fRnrMod->fMaxR + TMath::Abs(fR)) { | |
106 | crosR = true; | |
5a5a1232 | 107 | } |
bbfaa1c4 | 108 | // printf("end steps in helix line %d \n", NMax); |
109 | } | |
5a5a1232 | 110 | |
111 | ||
bbfaa1c4 | 112 | void Step(Float_t &px, Float_t &py, Float_t &/*pz*/) |
113 | { | |
114 | fV.t += fTimeStep; | |
115 | fV.x += (px*sin - py*(1 - cos))/fA + x_off; | |
116 | fV.y += (py*sin + px*(1 - cos))/fA + y_off; | |
117 | fV.z += fLam*TMath::Abs(fR*fPhiStep); | |
118 | fPoints->push_back(fV); | |
119 | Float_t px_t = px*cos - py*sin ; | |
120 | Float_t py_t = py*cos + px*sin ; | |
121 | px = px_t; | |
122 | py = py_t; | |
b855ed81 | 123 | ++fN; |
bbfaa1c4 | 124 | } |
5a5a1232 | 125 | |
126 | ||
bbfaa1c4 | 127 | Bool_t LoopToVertex(Float_t &px, Float_t &py, Float_t &pz, |
f3e27855 | 128 | Float_t ex, Float_t ey, Float_t ez) |
bbfaa1c4 | 129 | { |
130 | Float_t p0x = px, p0y = py; | |
131 | Float_t zs = fLam*TMath::Abs(fR*fPhiStep); | |
132 | Float_t fnsteps = (ez - fV.z)/zs; | |
133 | Int_t nsteps = Int_t((ez - fV.z)/zs); | |
134 | Float_t sinf = TMath::Sin(fnsteps*fPhiStep); | |
135 | Float_t cosf = TMath::Cos(fnsteps*fPhiStep); | |
5a5a1232 | 136 | |
bbfaa1c4 | 137 | { |
b855ed81 | 138 | if (nsteps > 0) |
139 | { | |
bbfaa1c4 | 140 | Float_t xf = fV.x + (px*sinf - py*(1 - cosf))/fA; |
141 | Float_t yf = fV.y + (py*sinf + px*(1 - cosf))/fA; | |
142 | x_off = (ex - xf)/fnsteps; | |
143 | y_off = (ey - yf)/fnsteps; | |
144 | Float_t xforw, yforw, zforw; | |
b855ed81 | 145 | for (Int_t l=0; l<nsteps; l++) |
146 | { | |
bbfaa1c4 | 147 | xforw = fV.x + (px*sin - py*(1 - cos))/fA + x_off; |
148 | yforw = fV.y + (py*sin + px*(1 - cos))/fA + y_off; | |
149 | zforw = fV.z + fLam*TMath::Abs(fR*fPhiStep); | |
150 | if ((xforw*xforw+yforw*yforw > fRnrMod->fMaxR*fRnrMod->fMaxR) ||(TMath::Abs(zforw) > fRnrMod->fMaxZ) ) { | |
151 | return false; | |
152 | } | |
b855ed81 | 153 | Step(px,py,pz); |
154 | if (fN > NMax) | |
155 | break; | |
bbfaa1c4 | 156 | } |
5a5a1232 | 157 | |
5a5a1232 | 158 | } |
bbfaa1c4 | 159 | // set time to the end point |
160 | fV.t += TMath::Sqrt((fV.x-ex)*(fV.x-ex)+(fV.y-ey)*(fV.y-ey) +(fV.z-ez)*(fV.z-ez))/fVelocity; | |
161 | fV.x = ex; fV.y = ey; fV.z = ez; | |
162 | fPoints->push_back(fV); | |
163 | } | |
5a5a1232 | 164 | |
bbfaa1c4 | 165 | { // fix momentum in the remaining part |
166 | Float_t cosr = TMath::Cos((fnsteps-nsteps)*fPhiStep); | |
167 | Float_t sinr = TMath::Sin((fnsteps-nsteps)*fPhiStep); | |
168 | Float_t px_t = px*cosr - py*sinr ; | |
169 | Float_t py_t = py*cosr + px*sinr ; | |
170 | px = px_t; | |
171 | py = py_t; | |
172 | } | |
173 | { // calculate direction of faked px,py | |
174 | Float_t pxf = (p0x*cosf - p0y*sinf)/TMath::Abs(fA) + x_off/fPhiStep; | |
175 | Float_t pyf = (p0y*cosf + p0x*sinf)/TMath::Abs(fA) + y_off/fPhiStep; | |
176 | Float_t fac = TMath::Sqrt(p0x*p0x + p0y*p0y)/TMath::Sqrt(pxf*pxf + pyf*pyf); | |
177 | px = fac*pxf; | |
178 | py = fac*pyf; | |
5a5a1232 | 179 | } |
bbfaa1c4 | 180 | return true; |
181 | } | |
5a5a1232 | 182 | |
bbfaa1c4 | 183 | Bool_t LoopToBounds(Float_t &px, Float_t &py, Float_t &pz) |
184 | { | |
185 | // printf("MC helix loop_to_bounds\n"); | |
186 | SetBounds(); | |
b855ed81 | 187 | if (NMax > 0) |
188 | { | |
bbfaa1c4 | 189 | // printf("NMAx MC helix loop_to_bounds\n"); |
bbfaa1c4 | 190 | Float_t xforw,yforw,zforw; |
b855ed81 | 191 | while (fN < NMax) |
192 | { | |
bbfaa1c4 | 193 | xforw = fV.x + (px*sin - py*(1 - cos))/fA + x_off; |
194 | yforw = fV.y + (py*sin + px*(1 - cos))/fA + y_off; | |
195 | zforw = fV.z + fLam*TMath::Abs(fR*fPhiStep); | |
5a5a1232 | 196 | |
bbfaa1c4 | 197 | if ((crosR && (xforw*xforw+yforw*yforw > fRnrMod->fMaxR*fRnrMod->fMaxR)) ||(TMath::Abs(zforw) > fRnrMod->fMaxZ)) { |
198 | return false; | |
199 | } | |
200 | Step(px,py,pz); | |
5a5a1232 | 201 | } |
bbfaa1c4 | 202 | return true; |
5a5a1232 | 203 | } |
bbfaa1c4 | 204 | return false; |
205 | } | |
206 | }; | |
5a5a1232 | 207 | |
bbfaa1c4 | 208 | /**************************************************************************/ |
209 | // LINE | |
210 | /**************************************************************************/ | |
5a5a1232 | 211 | |
bbfaa1c4 | 212 | struct MCLine : public MCStruct |
213 | { | |
214 | MCLine(TrackRnrStyle* rs, MCVertex* v0 ,Float_t vel, std::vector<MCVertex>* tpv): | |
215 | MCStruct(rs, v0 , vel, tpv) | |
216 | {} | |
5a5a1232 | 217 | |
bbfaa1c4 | 218 | Bool_t InBounds(Float_t ex, Float_t ey, Float_t ez) |
219 | { | |
220 | if(TMath::Abs(ez) > fRnrMod->fMaxZ || | |
221 | ex*ex + ey*ey > fRnrMod->fMaxR*fRnrMod->fMaxR) | |
222 | return false; | |
223 | else | |
224 | return true; | |
225 | } | |
5a5a1232 | 226 | |
bbfaa1c4 | 227 | void GotoVertex(Float_t x1, Float_t y1, Float_t z1) |
228 | { | |
bbfaa1c4 | 229 | fV.t += TMath::Sqrt((fV.x-x1)*(fV.x-x1)+(fV.y-y1)*(fV.y-y1)+(fV.z-z1)*(fV.z-z1))/fVelocity; |
230 | fV.x=x1; fV.y=y1; fV.z=z1; | |
231 | fPoints->push_back(fV); | |
232 | } | |
5a5a1232 | 233 | |
234 | ||
bbfaa1c4 | 235 | void GotoBounds( Float_t px, Float_t py, Float_t pz) |
236 | { | |
237 | Float_t tZ = 0,Tb = 0; | |
238 | // time where particle intersect +/- fMaxZ | |
239 | if (pz > 0) { | |
240 | tZ = (fRnrMod->fMaxZ - fV.z)/pz; | |
241 | } | |
242 | else if (pz < 0 ) { | |
243 | tZ = (-1)*(fRnrMod->fMaxZ + fV.z)/pz; | |
244 | } | |
245 | // time where particle intersects cylinder | |
246 | Float_t tR=0; | |
247 | Double_t a = px*px + py*py; | |
248 | Double_t b = 2*(fV.x*px + fV.y*py); | |
249 | Double_t c = fV.x*fV.x + fV.y*fV.y - fRnrMod->fMaxR*fRnrMod->fMaxR; | |
250 | Double_t D = b*b - 4*a*c; | |
251 | if(D >= 0) { | |
252 | Double_t D_sqrt=TMath::Sqrt(D); | |
253 | tR = ( -b - D_sqrt )/(2*a); | |
254 | if( tR < 0) { | |
255 | tR = ( -b + D_sqrt )/(2*a); | |
5a5a1232 | 256 | } |
257 | ||
bbfaa1c4 | 258 | // compare the two times |
259 | Tb = tR < tZ ? tR : tZ; | |
260 | } else { | |
261 | Tb = tZ; | |
5a5a1232 | 262 | } |
bbfaa1c4 | 263 | |
264 | GotoVertex(fV.x+px*Tb, fV.y+py*Tb, fV.z+ pz*Tb); | |
265 | } | |
266 | }; // struct Line | |
5a5a1232 | 267 | |
268 | ||
269 | } // namespace Reve | |
270 | ||
271 | #endif |