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