]>
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; | |
123 | fN++; | |
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 | { |
bbfaa1c4 | 138 | if(nsteps > 0){ |
139 | Float_t xf = fV.x + (px*sinf - py*(1 - cosf))/fA; | |
140 | Float_t yf = fV.y + (py*sinf + px*(1 - cosf))/fA; | |
141 | x_off = (ex - xf)/fnsteps; | |
142 | y_off = (ey - yf)/fnsteps; | |
143 | Float_t xforw, yforw, zforw; | |
144 | for (Int_t l=0; l<nsteps; l++) { | |
145 | xforw = fV.x + (px*sin - py*(1 - cos))/fA + x_off; | |
146 | yforw = fV.y + (py*sin + px*(1 - cos))/fA + y_off; | |
147 | zforw = fV.z + fLam*TMath::Abs(fR*fPhiStep); | |
148 | if ((xforw*xforw+yforw*yforw > fRnrMod->fMaxR*fRnrMod->fMaxR) ||(TMath::Abs(zforw) > fRnrMod->fMaxZ) ) { | |
149 | return false; | |
150 | } | |
151 | Step(px,py,pz); | |
152 | } | |
5a5a1232 | 153 | |
5a5a1232 | 154 | } |
bbfaa1c4 | 155 | // set time to the end point |
156 | fV.t += TMath::Sqrt((fV.x-ex)*(fV.x-ex)+(fV.y-ey)*(fV.y-ey) +(fV.z-ez)*(fV.z-ez))/fVelocity; | |
157 | fV.x = ex; fV.y = ey; fV.z = ez; | |
158 | fPoints->push_back(fV); | |
159 | } | |
5a5a1232 | 160 | |
bbfaa1c4 | 161 | { // fix momentum in the remaining part |
162 | Float_t cosr = TMath::Cos((fnsteps-nsteps)*fPhiStep); | |
163 | Float_t sinr = TMath::Sin((fnsteps-nsteps)*fPhiStep); | |
164 | Float_t px_t = px*cosr - py*sinr ; | |
165 | Float_t py_t = py*cosr + px*sinr ; | |
166 | px = px_t; | |
167 | py = py_t; | |
168 | } | |
169 | { // calculate direction of faked px,py | |
170 | Float_t pxf = (p0x*cosf - p0y*sinf)/TMath::Abs(fA) + x_off/fPhiStep; | |
171 | Float_t pyf = (p0y*cosf + p0x*sinf)/TMath::Abs(fA) + y_off/fPhiStep; | |
172 | Float_t fac = TMath::Sqrt(p0x*p0x + p0y*p0y)/TMath::Sqrt(pxf*pxf + pyf*pyf); | |
173 | px = fac*pxf; | |
174 | py = fac*pyf; | |
5a5a1232 | 175 | } |
bbfaa1c4 | 176 | return true; |
177 | } | |
5a5a1232 | 178 | |
bbfaa1c4 | 179 | Bool_t LoopToBounds(Float_t &px, Float_t &py, Float_t &pz) |
180 | { | |
181 | // printf("MC helix loop_to_bounds\n"); | |
182 | SetBounds(); | |
183 | if(NMax > 0){ | |
184 | // printf("NMAx MC helix loop_to_bounds\n"); | |
bbfaa1c4 | 185 | Float_t xforw,yforw,zforw; |
186 | while(fN < NMax){ | |
187 | xforw = fV.x + (px*sin - py*(1 - cos))/fA + x_off; | |
188 | yforw = fV.y + (py*sin + px*(1 - cos))/fA + y_off; | |
189 | zforw = fV.z + fLam*TMath::Abs(fR*fPhiStep); | |
5a5a1232 | 190 | |
bbfaa1c4 | 191 | if ((crosR && (xforw*xforw+yforw*yforw > fRnrMod->fMaxR*fRnrMod->fMaxR)) ||(TMath::Abs(zforw) > fRnrMod->fMaxZ)) { |
192 | return false; | |
193 | } | |
194 | Step(px,py,pz); | |
5a5a1232 | 195 | } |
bbfaa1c4 | 196 | return true; |
5a5a1232 | 197 | } |
bbfaa1c4 | 198 | return false; |
199 | } | |
200 | }; | |
5a5a1232 | 201 | |
bbfaa1c4 | 202 | /**************************************************************************/ |
203 | // LINE | |
204 | /**************************************************************************/ | |
5a5a1232 | 205 | |
bbfaa1c4 | 206 | struct MCLine : public MCStruct |
207 | { | |
208 | MCLine(TrackRnrStyle* rs, MCVertex* v0 ,Float_t vel, std::vector<MCVertex>* tpv): | |
209 | MCStruct(rs, v0 , vel, tpv) | |
210 | {} | |
5a5a1232 | 211 | |
bbfaa1c4 | 212 | Bool_t InBounds(Float_t ex, Float_t ey, Float_t ez) |
213 | { | |
214 | if(TMath::Abs(ez) > fRnrMod->fMaxZ || | |
215 | ex*ex + ey*ey > fRnrMod->fMaxR*fRnrMod->fMaxR) | |
216 | return false; | |
217 | else | |
218 | return true; | |
219 | } | |
5a5a1232 | 220 | |
bbfaa1c4 | 221 | void GotoVertex(Float_t x1, Float_t y1, Float_t z1) |
222 | { | |
bbfaa1c4 | 223 | fV.t += TMath::Sqrt((fV.x-x1)*(fV.x-x1)+(fV.y-y1)*(fV.y-y1)+(fV.z-z1)*(fV.z-z1))/fVelocity; |
224 | fV.x=x1; fV.y=y1; fV.z=z1; | |
225 | fPoints->push_back(fV); | |
226 | } | |
5a5a1232 | 227 | |
228 | ||
bbfaa1c4 | 229 | void GotoBounds( Float_t px, Float_t py, Float_t pz) |
230 | { | |
231 | Float_t tZ = 0,Tb = 0; | |
232 | // time where particle intersect +/- fMaxZ | |
233 | if (pz > 0) { | |
234 | tZ = (fRnrMod->fMaxZ - fV.z)/pz; | |
235 | } | |
236 | else if (pz < 0 ) { | |
237 | tZ = (-1)*(fRnrMod->fMaxZ + fV.z)/pz; | |
238 | } | |
239 | // time where particle intersects cylinder | |
240 | Float_t tR=0; | |
241 | Double_t a = px*px + py*py; | |
242 | Double_t b = 2*(fV.x*px + fV.y*py); | |
243 | Double_t c = fV.x*fV.x + fV.y*fV.y - fRnrMod->fMaxR*fRnrMod->fMaxR; | |
244 | Double_t D = b*b - 4*a*c; | |
245 | if(D >= 0) { | |
246 | Double_t D_sqrt=TMath::Sqrt(D); | |
247 | tR = ( -b - D_sqrt )/(2*a); | |
248 | if( tR < 0) { | |
249 | tR = ( -b + D_sqrt )/(2*a); | |
5a5a1232 | 250 | } |
251 | ||
bbfaa1c4 | 252 | // compare the two times |
253 | Tb = tR < tZ ? tR : tZ; | |
254 | } else { | |
255 | Tb = tZ; | |
5a5a1232 | 256 | } |
bbfaa1c4 | 257 | |
258 | GotoVertex(fV.x+px*Tb, fV.y+py*Tb, fV.z+ pz*Tb); | |
259 | } | |
260 | }; // struct Line | |
5a5a1232 | 261 | |
262 | ||
263 | } // namespace Reve | |
264 | ||
265 | #endif |