Automatic treatment of the magnetic field value
[u/mrichter/AliRoot.git] / ITS / AliITSTrackV1.cxx
1 //   ITS Track Class
2 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
3 // It contain all the usefull information for the track and the method to calculate, modify or extract them
4 // The track is mainly caracterized by the state vector of elements (fX0, fX1, fX2, fX3, fX4) and the
5 // corresponding covariance matrix of elements (C00, C10, ..... C44) that is triangular
6 //
7 #include <iostream.h>
8 #include <TMath.h>
9 #include <TVector.h>
10 #include <TObjArray.h>
11 #include "AliRun.h"
12 #include "AliITSRad.h"
13 #include "AliITSTrackV1.h"
14 #include "AliGenerator.h"
15 #include "AliMagF.h"
16
17
18 ClassImp(AliITSTrackV1)
19
20 AliITSTrackV1::AliITSTrackV1() {
21 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
22 // default constructor   
23  
24   fC00=fC10=fC11=fC20=fC21=fC22=fC30=fC31=fC32=fC33=fC40=fC41=fC42=fC43=fC44=0.;
25   flistCluster = new TObjArray; 
26   fNumClustInTrack =0;
27   fChi2=-1;
28   flabel =0; 
29   fVertex.ResizeTo(3); 
30   fErrorVertex.ResizeTo(3);
31   fLayer = -1; 
32   fClusterInTrack = new TMatrix(6,9);
33   Int_t i;
34   for(i=0; i<6; i++) (*fClusterInTrack)(i,6)=(*fClusterInTrack)(i,7)=
35                            (*fClusterInTrack)(i,8)=-1.;
36   frtrack=0.; 
37   fd2.ResizeTo(6);
38   ftgl2.ResizeTo(6); 
39   fdtgl.ResizeTo(6);
40   
41 //////////////////////////////////////// gets magnetic field factor ////////////////////////////////
42
43   AliMagF * fieldPointer = gAlice->Field();
44   fFieldFactor = (Double_t)fieldPointer->Factor();
45   //cout<< " field factor = "<<fFieldFactor<<"\n"; getchar();
46
47 /////////////////////////////////////////////////////////////////////////////////////////////////////////
48    
49 }
50
51
52  
53 AliITSTrackV1::AliITSTrackV1(const AliITSTrackV1 &cobj) {
54 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
55 // copy constructor    
56
57   fClusterInTrack = new TMatrix(6,9);
58   Int_t i;
59   for(i=0; i<6; i++) (*fClusterInTrack)(i,6)=(*fClusterInTrack)(i,7)=
60                            (*fClusterInTrack)(i,8)=-1.;
61   flistCluster = new TObjArray; 
62   fVertex.ResizeTo(3); 
63   fErrorVertex.ResizeTo(3);
64   fVertex = cobj.fVertex;
65   fErrorVertex = cobj.fErrorVertex;
66   flabel = cobj.flabel;
67   fLayer=cobj.fLayer;
68   fTPCtrack = cobj.fTPCtrack;
69   fNumClustInTrack = cobj.fNumClustInTrack; 
70   fChi2= cobj.fChi2;
71   fX0=cobj.fX0; fX1=cobj.fX1; fX2=cobj.fX2; fX3=cobj.fX3; fX4=cobj.fX4; 
72   frtrack=cobj.frtrack;
73   fDv=cobj.fDv;
74   fZv=cobj.fZv;
75   fsigmaDv=cobj.fsigmaDv;
76   fsigmaZv=cobj.fsigmaZv;
77   fd2.ResizeTo(6);
78   ftgl2.ResizeTo(6); 
79   fdtgl.ResizeTo(6); 
80   fd2=cobj.fd2;
81   ftgl2=cobj.ftgl2;
82   fdtgl=cobj.fdtgl;
83     
84   fC00=cobj.fC00; fC10=cobj.fC10; fC11=cobj.fC11; fC20=cobj.fC20; fC21=cobj.fC21;
85   fC22=cobj.fC22; fC30=cobj.fC30; fC31=cobj.fC31; fC32=cobj.fC32; fC33=cobj.fC33; 
86   fC40=cobj.fC40; fC41=cobj.fC41; fC42=cobj.fC42; fC43=cobj.fC43; fC44=cobj.fC44; 
87  
88   *fClusterInTrack = *cobj.fClusterInTrack;
89   
90   fFieldFactor=cobj.fFieldFactor;
91  
92   for(i=0; i<cobj.flistCluster->GetSize(); i++) 
93     flistCluster->AddLast(cobj.flistCluster->At(i));
94  
95 }
96
97 AliITSTrackV1::AliITSTrackV1(AliTPCtrack &obj)
98
99 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
100 // special constructor to convert a TPC track into an ITS track
101
102 //////////////////////////////////////// gets magnetic field factor ////////////////////////////////
103
104   AliMagF * fieldPointer = gAlice->Field();
105   fFieldFactor = (Double_t)fieldPointer->Factor();
106   //cout<< " field factor = "<<fFieldFactor<<"\n"; getchar();
107
108 /////////////////////////////////////////////////////////////////////////////////////////////////////////
109
110
111   fTPCtrack = &obj;
112   fVertex.ResizeTo(3); 
113   fErrorVertex.ResizeTo(3);
114
115   fd2.ResizeTo(6);
116   ftgl2.ResizeTo(6); 
117   fdtgl.ResizeTo(6); 
118   AliGenerator *gener = gAlice->Generator();
119   Float_t vxg,vyg,vzg;
120   gener->GetOrigin(vxg,vyg,vzg);
121   
122  
123   fVertex(0)=(Double_t)vxg;
124   fVertex(1)=(Double_t)vyg; 
125   fVertex(2)=(Double_t)vzg;    
126   
127   fLayer = 7;
128   //fmCovariance = new TMatrix(5,5);
129   fClusterInTrack = new TMatrix(6,9);
130   
131   Int_t i;
132   for(i=0; i<6; i++) (*fClusterInTrack)(i,6)=(*fClusterInTrack)(i,7)=
133                            (*fClusterInTrack)(i,8)=-1.;
134   flistCluster = new TObjArray; 
135   fNumClustInTrack = 0; 
136   LmTPC(); 
137
138 }
139
140
141 AliITSTrackV1::~AliITSTrackV1() {
142 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
143 //destructor
144  
145   if(flistCluster) delete flistCluster; 
146   if(fClusterInTrack) delete fClusterInTrack;
147
148 }     
149
150 void AliITSTrackV1::PutCElements(Double_t C00, Double_t C10, Double_t C11, Double_t C20, Double_t C21, 
151 Double_t C22, Double_t C30, Double_t C31, Double_t C32, Double_t C33, Double_t C40, 
152 Double_t C41, Double_t C42, Double_t C43, Double_t C44){
153 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
154 // puts elements of covariance matrix
155
156   fC00=C00; fC10=C10; fC11=C11; fC20=C20; fC21=C21; fC22=C22; fC30=C30; fC31=C31; fC32=C32; fC33=C33;
157   fC40=C40; fC41=C41; fC42=C42; fC43=C43; fC44=C44; 
158 }
159   
160 void AliITSTrackV1::GetCElements(Double_t &C00, Double_t &C10, Double_t &C11, Double_t &C20, Double_t &C21, 
161 Double_t &C22, Double_t &C30, Double_t &C31, Double_t &C32, Double_t &C33, Double_t &C40, 
162 Double_t &C41, Double_t &C42, Double_t &C43, Double_t &C44) const {
163 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
164 // gets elements of covariance matrix
165
166
167   C00=fC00; C10=fC10; C11=fC11; C20=fC20; C21=fC21; C22=fC22; C30=fC30; C31=fC31; C32=fC32; C33=fC33;
168   C40=fC40; C41=fC41; C42=fC42; C43=fC43; C44=fC44; 
169
170 }
171
172 void AliITSTrackV1::GetXElements(Double_t &X0, Double_t &X1, Double_t &X2, Double_t &X3, Double_t &X4) const {
173 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
174 // gets the elements of the state vector
175   X0=fX0; X1=fX1; X2=fX2; X3=fX3; X4=fX4; 
176 }
177
178 void AliITSTrackV1::PutXElements(Double_t X0, Double_t X1, Double_t X2, Double_t X3, Double_t X4){
179 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
180 // puts the elements of the state vector
181   fX0=X0; fX1=X1; fX2=X2; fX3=X3; fX4=X4; 
182 }   
183
184 void AliITSTrackV1::LmTPC() {
185 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
186 // Transform the TPC state vector from TPC-local to master and build a new state vector ITS-type
187 // The covariance matrix is also modified accordingly      
188
189    
190   Double_t alpha = fTPCtrack->GetAlpha(); 
191  
192   //printf("LmTPC: alpha %f\n",alpha); 
193
194   Double_t yTPC   = fTPCtrack->GetY();
195   Double_t zTPC   = fTPCtrack->GetZ();
196   Double_t cTPC   = fTPCtrack->GetC();
197   Double_t etaTPC = fTPCtrack->GetEta();
198   Double_t tglTPC = fTPCtrack->GetTgl();
199   
200
201   Double_t xm, ym, zm;
202   Double_t sina = TMath::Sin(alpha);
203   Double_t cosa = TMath::Cos(alpha);
204   Double_t xl= fTPCtrack->GetX();
205   xm = xl * cosa - yTPC*sina;
206   ym = xl * sina + yTPC*cosa;
207   zm = zTPC;  
208   //cout<<" xl e alpha = "<<xl<<" "<<alpha<<"\n"; getchar(); 
209
210   Double_t x0m,y0m;
211
212   ///////////////////////////////////// determine yo //////////////////////////////////////////////////
213   
214   Double_t vxl=fVertex(0)*cosa+fVertex(1)*sina;
215   Double_t vyl= -fVertex(0)*sina+fVertex(1)*cosa;
216   Double_t xo,yo, signy;
217   Double_t r = 1./cTPC;
218   xo =  etaTPC / cTPC;
219   fxoTPC=xo;
220   Double_t yo1, yo2, diffsq1, diffsq2;  
221   yo1 = yTPC +  TMath::Sqrt(r*r - (xl-xo)*(xl-xo)); 
222   yo2 = yTPC -  TMath::Sqrt(r*r - (xl-xo)*(xl-xo));   
223   diffsq1=TMath::Abs((yo1-vyl)*(yo1-vyl)+(xo-vxl)*(xo-vxl)-r*r);
224   diffsq2=TMath::Abs((yo2-vyl)*(yo2- vyl)+(xo-vxl)*(xo-vxl)-r*r);
225   if(diffsq1<diffsq2) {yo=yo1; signy=1.;} else {yo=yo2; signy=-1.;};
226   
227   ////////////////////////////////////////////////////////////////////////////////////////////////////
228
229   x0m = xo * cosa - yo * sina;
230   y0m = xo * sina + yo * cosa;  
231
232   frtrack=TMath::Sqrt(xm*xm+ym*ym);
233   Double_t pigre=TMath::Pi();
234   Double_t phi=0.0;
235   if(ym == 0. || xm == 0.) {
236     if(ym == 0. && xm == 0.) {cout << "  Error in AliITSTrackV1::LmTPC x=0 and y=0 \n"; getchar();}
237     if(ym ==0. && xm>0.) phi=0.;
238     if(ym==0. && xm<0.) phi=pigre;
239     if(xm==0 && ym>0.) phi=pigre/2.;
240     if(xm==0 && ym<0.) phi=1.5*pigre;   
241   }
242   else {
243     if (xm>0. && ym>0.) phi=TMath::ATan(ym/xm);
244     if (xm<0. && ym>0.) phi=pigre+TMath::ATan(ym/xm);
245     if (xm<0. && ym<0.) phi=pigre+TMath::ATan(ym/xm); 
246     if (xm>0. && ym<0.) phi=(2*pigre)+TMath::ATan(ym/xm);     
247   };
248   if(phi<0. || phi>(2*pigre)) {cout<<"attention error on phi in  AliITSTrackV1:LmTPC \n"; getchar();}  
249
250   fX0=phi;
251   fX1=zm;
252   fX3=tglTPC;
253   fX4=cTPC;
254
255   
256   Double_t dd=TMath::Sqrt((x0m-fVertex(0))*(x0m-fVertex(0))+(y0m-fVertex(1))*(y0m-fVertex(1)));
257   Double_t signdd;
258   if (r>0) signdd=1.; else signdd=-1.;
259   fX2=signdd*dd-r;
260   //cout<<" fvertex = "<<fVertex(0)<<" "<<fVertex(1)<<" "<<fVertex(2)<<"\n";
261       
262   Double_t cov[15];
263   fTPCtrack->GetCovariance(cov);
264
265   Double_t  dfidy, dDdy, dDdC, dDdeta;
266
267   dfidy=(xm*cosa+ym*sina)/(frtrack*frtrack);
268   dDdy=signdd*((y0m-fVertex(1))*cosa-(x0m-fVertex(0))*sina)/dd;
269   Double_t dyodr=signy*(r+(xl-xo)*etaTPC)/TMath::Sqrt(r*r-(xl-xo)*(xl-xo));
270   Double_t dyomdr=sina*etaTPC+cosa*dyodr;
271   Double_t dxomdr=cosa*etaTPC-sina*dyodr;
272   Double_t ddddR=((x0m-fVertex(0))*dxomdr+(y0m-fVertex(1))*dyomdr)/dd;
273   dDdC=-r*r*(signdd*ddddR-1.);
274   Double_t dyoldxol=signy*(xl-xo)/TMath::Sqrt(r*r-(xl-xo)*(xl-xo));
275   Double_t dxomdeta=r*(cosa-sina*dyoldxol);
276   Double_t dyomdeta=r*(sina+cosa*dyoldxol);
277   dDdeta=signdd*((x0m-fVertex(0))*dxomdeta+(y0m-fVertex(1))*dyomdeta)/dd;
278   
279   Double_t f00=dfidy;
280   Double_t f20=dDdy;    
281   Double_t f22=dDdC;    
282   Double_t f23=dDdeta;
283   
284   Double_t t00=cov[0]*f00;
285   Double_t t02=cov[0]*f20+cov[6]*f22+cov[3]*f23;
286   Double_t t20=cov[6]*f00;
287   Double_t t22=cov[6]*f20+cov[9]*f22+cov[8]*f23;
288   
289   fC00=f00*t00;
290   fC10=cov[1]*f00;
291   fC11=cov[2];
292   fC20=f20*t00+f22*t20+f23*cov[3]*f00;
293   fC21=f20*cov[1]+f22*cov[7]+f23*cov[4];
294   fC22=f20*t02+f22*t22+f23*(cov[3]*f20+cov[8]*f22+cov[5]*f23);
295   fC30=cov[10]*f00;
296   fC31=cov[11];
297   fC32=cov[10]*f20+cov[13]*f22+cov[12]*f23;
298   fC33=cov[14];
299   fC40=t20;
300   fC41=cov[7];
301   fC42=t22;
302   fC43=cov[13];
303   fC44=cov[9];
304   
305   //cout<<" C32 e C44 = "<<fC32<<" "<<fC44<<"\n"; getchar();
306    
307 }
308
309
310 AliITSTrackV1 &AliITSTrackV1::operator=(AliITSTrackV1 obj) {
311 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
312 // assignement operator
313
314   delete flistCluster;
315   delete fClusterInTrack;
316   fClusterInTrack = new TMatrix(6,9);
317   flistCluster = new TObjArray; 
318   flabel = obj.flabel;
319   fTPCtrack = obj.fTPCtrack; 
320   fNumClustInTrack = obj.fNumClustInTrack; 
321   fChi2= obj.fChi2;
322   fVertex=obj.fVertex;
323   fErrorVertex=obj.fErrorVertex;
324   fX0=obj.fX0; fX1=obj.fX1; fX2=obj.fX2; fX3=obj.fX3; fX4=obj.fX4;  
325   fLayer=obj.fLayer;
326   frtrack=obj.frtrack;
327   fDv=obj.fDv;
328   fZv=obj.fZv;
329   fsigmaDv=obj.fsigmaDv;
330   fsigmaZv=obj.fsigmaZv;
331   fd2=obj.fd2;
332   ftgl2=obj.ftgl2; 
333   fdtgl=obj.fdtgl; 
334   
335   fC00=obj.fC00; fC10=obj.fC10; fC11=obj.fC11; fC20=obj.fC20; fC21=obj.fC21;
336   fC22=obj.fC22; fC30=obj.fC30; fC31=obj.fC31; fC32=obj.fC32; fC33=obj.fC33; 
337   fC40=obj.fC40; fC41=obj.fC41; fC42=obj.fC42; fC43=obj.fC43; fC44=obj.fC44; 
338   
339   
340   *fClusterInTrack = *obj.fClusterInTrack;
341   Int_t i;
342   for(i=0; i<obj.flistCluster->GetSize(); i++) flistCluster->AddLast(obj.flistCluster->At(i));
343
344   return *this;
345   
346 }
347
348 void AliITSTrackV1::PutCluster(Int_t layerc, TVector vecclust) {
349 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
350 // put information for clusters
351   
352   (*fClusterInTrack)(layerc,0) = vecclust(0);
353   (*fClusterInTrack)(layerc,1) = vecclust(1);
354   (*fClusterInTrack)(layerc,2) = vecclust(2);
355   (*fClusterInTrack)(layerc,3) = vecclust(3);
356   (*fClusterInTrack)(layerc,4) = vecclust(4);  
357   (*fClusterInTrack)(layerc,5) = vecclust(5);
358   (*fClusterInTrack)(layerc,6) = vecclust(6);
359   (*fClusterInTrack)(layerc,7) = vecclust(7);
360   (*fClusterInTrack)(layerc,8) = vecclust(8);
361
362 }
363
364
365 void AliITSTrackV1::GetClusters() {
366 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
367 // prints the clusters belonging to the current track
368
369   TMatrix a(*fClusterInTrack);
370   TMatrix b(6,3);
371   Int_t i;
372   for(i=0;i<6; i++){
373    b(i,0)=a(i,6); b(i,1)=a(i,7); b(i,2)=a(i,8); 
374   }
375   a.Print();
376  // b.Print(); 
377  
378 }
379
380
381 TVector AliITSTrackV1::GetLabTrack(Int_t lay) {
382 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
383 // gets the label of the track
384
385   TVector vecLabel(3);
386   vecLabel(0)=( (Float_t) (*fClusterInTrack)(lay,6) );
387   vecLabel(1)=( (Float_t) (*fClusterInTrack)(lay,7) );
388   vecLabel(2)=( (Float_t) (*fClusterInTrack)(lay,8) );
389   return vecLabel;
390 }
391
392 void AliITSTrackV1::Search(TVector VecTotLabref, Long_t &labref, Int_t &freq){
393 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
394 // define label
395
396   Int_t vecfreq[18];
397   
398   Int_t i,j;
399   for(i=0; i<18; i++) vecfreq[i]=0;
400   
401   for(i=0; i<18; i++) {
402     for(j=0; j<18; j++) {
403      // if(VecTotLabref(i) == 0.) VecTotLabref(i)=-3.;  //commentato il 5-3-2001
404       if( (VecTotLabref(i)>=0.) && (VecTotLabref(i)==VecTotLabref(j)) ) vecfreq[i]++;    
405     }  
406   }
407   Int_t imax=-1000;
408   Long_t  labdefault= (Long_t)1000000.;
409   freq=0;
410   for(i=0; i<18; i++) {
411     if(vecfreq[i]>freq) {freq=vecfreq[i]; imax=i;}  
412   }
413   if(imax<0) labref=labdefault; else labref=(Long_t) VecTotLabref(imax);
414
415
416
417 void AliITSTrackV1::Propagation(Double_t rk) {
418 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
419 //Propagation of track 
420   Double_t duepi=2.*TMath::Pi();
421   Double_t rkm1=frtrack;
422 //cout<<" rk e rkm1 dentro Propagation   "<<rk<<" "<<rkm1<<"\n";
423   
424   //
425   Double_t aAk=ArgA(rk), aAkm1=ArgA(rkm1);
426   Double_t ak=Arga(rk), akm1=Arga(rkm1);
427   fX0+=TMath::ASin(aAk)-TMath::ASin(aAkm1);
428
429   if(fX0>duepi) fX0-=duepi;
430   if(fX0<0.) fX0+=duepi;
431         
432   Double_t tgl=fX3;
433   Double_t c=fX4;
434   Double_t d=fX2;
435   Double_t cy=c/2;
436   fX1+=tgl/cy*(TMath::ASin(ak)-TMath::ASin(akm1));  
437   frtrack=rk;
438
439         
440   Double_t bk=ArgB(rk), bkm1=ArgB(rkm1);        
441   Double_t ck=ArgC(rk), ckm1=ArgC(rkm1);  
442
443   Double_t f02=ck/TMath::Sqrt(1.-aAk*aAk) - ckm1/TMath::Sqrt(1.-aAkm1*aAkm1);
444   Double_t f04=bk/TMath::Sqrt(1.-aAk*aAk) - bkm1/TMath::Sqrt(1.-aAkm1*aAkm1);   
445   Double_t f12=tgl*d*(1./rk - 1./rkm1);
446   Double_t f13=rk - rkm1;
447   
448
449   Double_t c00=fC00;
450   Double_t c10=fC10;
451   Double_t c11=fC11;
452   Double_t c20=fC20;
453   Double_t c21=fC21;
454   Double_t c22=fC22;
455   Double_t c30=fC30;
456   Double_t c31=fC31;
457   Double_t c32=fC32;
458   Double_t c33=fC33;
459   Double_t c40=fC40;
460   Double_t c41=fC41;
461   Double_t c42=fC42;
462   Double_t c43=fC43;
463   Double_t c44=fC44;
464   
465   Double_t r10=c10+c21*f02+c41*f04;
466   Double_t r20=c20+c22*f02+c42*f04;
467   Double_t r30=c30+c32*f02+c43*f04;
468   Double_t r40=c40+c42*f02+c44*f04;
469   Double_t r21=c21+c22*f12+c32*f13;
470   Double_t r31=c31+c32*f12+c33*f13;                                                                     
471   Double_t r41=c41+c42*f12+c43*f13;
472
473   fC00=c00+c20*f02+c40*f04+f02*r20+f04*r40;
474   fC10=r10+f12*r20+f13*r30;
475   fC11=c11+c21*f12+c31*f13+f12*r21+f13*r31;
476   fC20=r20;
477   fC21=r21;
478   fC30=r30;
479   fC31=r31;
480   fC40=r40;
481   fC41=r41;
482   
483 }
484
485 void AliITSTrackV1::AddEL(AliITSRad *rl, Double_t signdE, Bool_t flagtot, Double_t mass) {
486 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
487 //  add energy loss
488
489   TVector s(6);  
490   s(0)=0.0026+0.00283; s(1)=0.018; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
491   //s(0)=0.0026+0.00283*2.; s(1)=0.018*2.; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
492   //0.00277 is added in the first layer to take into account the energy loss in the beam pipe
493
494   //for(int k=0; k<6; k++) cout<<s(k)<<" "; cout<<"\n";
495   Int_t k;
496   for(k=0; k<6; k++) s(k)=s(k)*1.6;    //forint
497  
498   
499   Double_t phi=fX0;
500   
501   if(phi<0.174 ) s(5)=s(5)+0.012; 
502   if(phi>6.1 ) s(5)=s(5)+0.012; // to take into account rail 
503   if(phi>2.96 && phi<3.31 ) s(5)=s(5)+0.012;   
504   
505    
506   Double_t tgl=fX3;
507   Double_t theta=((TMath::Pi())/2.)-TMath::ATan(tgl);
508   //phi*=180./TMath::Pi();
509   //theta*=180./TMath::Pi();
510   //Double_t rad90=(TMath::Pi())/2.;
511   Double_t rad40=(TMath::Pi())*40./180.;
512   Double_t rad100=(TMath::Pi())*100/180;
513   Double_t rad360=(TMath::Pi())*2.;
514   Int_t imax=rl->Getimax();
515   Int_t jmax=rl->Getjmax();
516   Int_t i=(Int_t) ( (theta-rad40)/rad100*imax);
517   Int_t j=(Int_t) ( phi/rad360*jmax );
518   //Int_t i=(Int_t)( ((theta-((TMath::Pi())/4.))/((TMath::Pi())/2.))*imax );
519   //Int_t j=(Int_t)( (phi/((TMath::Pi())*2.))*jmax );
520   if(i<0) i=0;
521   if(i>=imax) i=imax-1;
522   if(j<0) j=0;
523   if(j>=jmax) j=jmax-1;
524   /*
525   s(0) = 0.0028/TMath::Sin(theta)+( rl->GetRadMatrix1() )(i,j);   // 0.0028 takes into account the beam pipe
526   s(1) = ( rl->GetRadMatrix2() )(i,j);
527   s(2) = ( rl->GetRadMatrix3() )(i,j);
528   s(3) = ( rl->GetRadMatrix4() )(i,j);
529   s(4) = ( rl->GetRadMatrix5() )(i,j);
530   s(5) = ( rl->GetRadMatrix6() )(i,j);
531   
532
533   */  
534   
535   //for(int k=0; k<6; k++) cout<<s(k)<<" "; getchar();
536     
537   //if(phi>60) {cout<<" phi = "<<phi<<"\n"; getchar();}
538   //if(theta<45 || theta>135) {cout<<" theta = "<<theta<<"\n"; getchar();}
539   //cout<<" dentro AddEl: phi, theta = "<<phi<<" "<<theta<<"\n"; getchar(); 
540     
541   Double_t cl=1.+fX3*fX3;  // cl=1/(cosl)**2 = 1 + (tgl)**2 
542   Double_t sqcl=TMath::Sqrt(cl);
543   Double_t pt=GetPt();
544      
545   Double_t p2=pt*pt*cl;
546   Double_t e=TMath::Sqrt(p2+mass*mass);
547   Double_t beta2=p2/(p2+mass*mass);
548   
549   Double_t dE;
550   if(flagtot) {
551     Double_t stot=s(0)+s(1)+s(2)+s(3)+s(4)+s(5);
552     dE=0.153/beta2*(log(5940*beta2/(1-beta2)) - beta2)*stot*21.82*sqcl;
553   } else {
554     dE=0.153/beta2*(log(5940*beta2/(1-beta2)) - beta2)*s(fLayer-1)*21.82*sqcl;
555   }   
556   dE=signdE*dE/1000.; 
557           
558   e+=dE;
559   Double_t p=TMath::Sqrt(e*e-mass*mass);   
560   Double_t sign=1.;
561   if(fX4 < 0.) sign=-1.; 
562   pt=sign*p/sqcl; 
563   Double_t cc=(0.299792458*0.2*fFieldFactor)/(pt*100.);
564   fX4=cc;
565   
566 }
567
568 void  AliITSTrackV1::Correct(Double_t rk) {
569 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
570 // correct track to take into account real geometry detector
571
572   Double_t duepi=2.*TMath::Pi();
573   Double_t rkm1=frtrack;
574   Double_t aAk=ArgA(rk), aAkm1=ArgA(rkm1);
575   Double_t ak=Arga(rk), akm1=Arga(rkm1);
576
577   fX0+=TMath::ASin(aAk)-TMath::ASin(aAkm1);
578   if(fX0>duepi) fX0-=duepi;
579   if(fX0<0.) fX0+=duepi;
580         
581   Double_t tgl=fX3;
582   Double_t c=fX4;
583   Double_t cy=c/2;
584   fX1+=tgl/cy*(TMath::ASin(ak)-TMath::ASin(akm1));
585   frtrack=rk;
586                 
587 }
588
589 void AliITSTrackV1::AddMS(AliITSRad *rl) {
590 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it       
591 //////////   Modification of the covariance matrix to take into account multiple scattering  ///////////     
592    
593   TVector s(6);
594
595   s(0)=0.0026+0.00283; s(1)=0.018; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
596  //  s(0)=0.0026+0.00283*2.; s(1)=0.018*2.; s(2)=0.0094; s(3)=0.0095; s(4)=0.0091; s(5)=0.0087;
597 //0.00277 is added in the first layer to take into account the energy loss in the beam pipe 
598
599   Int_t k;
600   for(k=0; k<6; k++) s(k)=s(k)*1.6;   // forint
601   
602   Double_t phi=fX0;
603  if(phi<0.174 ) s(5)=s(5)+0.012; //aggiunta provvisoria
604   if(phi>6.1 ) s(5)=s(5)+0.012; //aggiunta provvisoria  
605    if(phi>2.96 && phi< 3.31) s(5)=s(5)+0.012; //aggiunta provvisoria 
606       
607   Double_t tgl=fX3;
608   Double_t theta=((TMath::Pi())/2.)-TMath::ATan(tgl);
609   Double_t rad40=(TMath::Pi())*40./180.;      // rivedere
610   Double_t rad100=(TMath::Pi())*100/180;
611   Double_t rad360=(TMath::Pi())*2.;
612   Int_t imax=rl->Getimax();
613   Int_t jmax=rl->Getjmax();
614   Int_t i=(Int_t) ( (theta-rad40)/rad100*imax);
615   Int_t j=(Int_t) ( phi/rad360*jmax);
616
617   if(i<0) i=0;
618   if(i>=imax) i=imax-1;
619   if(j<0) j=0;
620   if(j>=jmax) j=jmax-1;
621   /*
622   s(0) = 0.0028/TMath::Sin(theta)+( rl->GetRadMatrix1() )(i,j);   // 0.0028 takes into account the beam pipe
623   s(1) = ( rl->GetRadMatrix2() )(i,j);
624   s(2) = ( rl->GetRadMatrix3() )(i,j);
625   s(3) = ( rl->GetRadMatrix4() )(i,j);
626   s(4) = ( rl->GetRadMatrix5() )(i,j);
627   s(5) = ( rl->GetRadMatrix6() )(i,j);
628    */   
629   Double_t mass=0.1396;
630   Int_t layer=(Int_t)GetLayer();
631   
632   Double_t cosl=TMath::Cos(TMath::ATan(tgl));   
633   Double_t d=fX2;
634   Double_t c=fX4;
635   Double_t cy=c/2.;
636   Double_t q20=1./(cosl*cosl);   
637   Double_t q30=c*tgl;
638    
639   Double_t q40=cy*(frtrack*frtrack-d*d)/(1.+ 2.*cy*d);
640   Double_t dd=d+cy*d*d-cy*frtrack*frtrack;
641   Double_t dprova=frtrack*frtrack - dd*dd;
642   Double_t q41=0.;
643   if(dprova>0.) q41=-1./cosl*TMath::Sqrt(dprova)/(1.+ 2.*cy*d);
644                 
645   Double_t p2=(GetPt()*GetPt())/(cosl*cosl);
646   Double_t beta2=p2/(p2+mass*mass);
647   Double_t theta2=14.1*14.1/(beta2*p2*1.e6)*(s(layer-1)/cosl);
648
649   fC22+=theta2*(q40*q40+q41*q41);
650   fC32+=theta2*q20*q40;
651   fC33+=theta2*q20*q20;
652   fC42+=theta2*q30*q40;
653   fC43+=theta2*q30*q20;
654   fC44+=theta2*q30*q30;
655     
656 }
657 void AliITSTrackV1::PrimaryTrack(AliITSRad *rl) {
658 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
659 // calculation of part of covariance matrix for vertex constraint
660
661   Double_t rlayer[6];
662         
663   rlayer[0]=4.; rlayer[1]=7.;  rlayer[2]=14.9;  rlayer[3]=23.8;  
664   rlayer[4]=39.1;  rlayer[5]=43.6;
665         
666   Double_t cy=fX4/2.;   
667   Double_t tgl=(fX1-fZv)*cy/TMath::ASin(cy*frtrack);
668   Double_t frtrack=1.;  
669   fX0=0.;
670   fX1=frtrack*tgl; 
671   fX2=fDv;
672   fX3=tgl; 
673
674   fC00=fC10=fC11=fC20=fC21=fC22=fC30=fC31=fC32=fC33=fC40=fC41=fC42=fC43=0.;
675
676   AddEL(rl,1.,1);
677   fLayer=0;
678   Int_t i;
679   for (i=0; i<6; i++) {
680     Propagation(rlayer[i]);
681     fLayer++;
682     fd2(i)=fC22;
683     ftgl2(i)=fC33;
684     fdtgl(i)=fC32; 
685     AddMS(rl);    
686     AddEL(rl,-1,0);        
687   }             
688 }       
689  
690   
691 Int_t AliITSTrackV1::DoNotCross(Double_t rk) const{
692 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it  
693 // determine if the track cross a layer
694
695   Double_t c=fX4;
696   Double_t d=fX2;
697   Double_t cy=c/2.;
698   return (TMath::Abs((cy*rk+(1.+cy*d)*d/rk)/(1.+2.*cy*d))>=1.)?1:0;
699 }
700   
701
702 Double_t AliITSTrackV1::ArgA(Double_t rk) const {
703 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
704 // quantity usefull in Propagation  
705   Double_t c=fX4;
706   Double_t d=fX2;
707   Double_t cy=c/2.;
708   Double_t arg=(cy*rk + (1 + cy*d)*d/rk)/(1.+ 2.*cy*d);         
709   if (TMath::Abs(arg) < 1.) return arg;
710   //cout<<"class AliITSTrack: ArgA out of range !\n";/* getchar();*/
711   return (arg>0) ? 0.99999999999 : -0.9999999999;
712 }
713    
714 Double_t AliITSTrackV1::Arga(Double_t rk) const {
715 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
716 // quantity usefull in Propagation 
717   Double_t c=fX4;
718   Double_t d=fX2;
719   Double_t cy=c/2.;
720   Double_t arg=(rk*rk - d*d)/(1.+ 2.*cy*d);             
721   if (arg<0.) {/*cout<<"class AliITSTrack: Arga out of range !\n";*/ arg=0.;} 
722   return cy*TMath::Sqrt(arg);
723 }
724         
725 Double_t AliITSTrackV1::ArgB(Double_t rk) const {
726 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
727 // quantity usefull in Propagation 
728   Double_t c=fX4;
729   Double_t d=fX2;
730   Double_t cy=c/2.;        
731   return (rk*rk - d*d)/(rk*(1.+ 2.*cy*d)*(1.+ 2.*cy*d));
732 }
733    
734 Double_t AliITSTrackV1::ArgC(Double_t rk) const {
735 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
736 // quantity usefull in Propagation 
737   Double_t c=fX4;
738   Double_t d=fX2;
739   Double_t cy=c/2.;             
740   return  (1./rk - 2.*cy*ArgA(rk)/(1.+ 2.*cy*d));
741 }
742
743
744