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