]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITStracking.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Set high water mark after last particle.
[u/mrichter/AliRoot.git] / ITS / AliITStracking.cxx
CommitLineData
dce1757b 1#include <TMath.h>
2#include <TList.h>
3#include <TTree.h>
4#include <TVector.h>
5#include <TMatrix.h>
6#include <TObjectTable.h>
7
8#include "AliITStracking.h"
9#include "AliRun.h"
10#include "AliITS.h"
11#include "AliITSgeom.h"
12#include "AliITSRecPoint.h"
13#include "AliITStrack.h"
14
15ClassImp(AliITStracking)
16
17
18AliITStracking::AliITStracking(TList *trackITSlist, AliITStrack *reference,
19 AliITS *aliITS, TObjArray *rpoints, Double_t Ptref, Int_t **vettid, Bool_t flagvert) {
20/////////////////////// This function perform the tracking in ITS detectors /////////////////////
21/////////////////////// reference is a pointer to the final best track /////////////////////
22//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
23// The authors thank Mariana Bondila to have help them to resolve some problems. July-2000
24
dce1757b 25 Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6;
cc6e75dd 26
27 Int_t index;
28 for(index =0; index<trackITSlist->GetSize(); index++) {
dce1757b 29 AliITStrack *trackITS = (AliITStrack *) trackITSlist->At(index);
cc6e75dd 30
dce1757b 31 if((*trackITS).GetLayer()==7) reference->SetChi2(10.223e140);
cc6e75dd 32 // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n";
33 // cout<<"fvtrack =" <<"\n";
34 // cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n";
35 // cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n";
36 // cout<< " Pt = "<<(*trackITS).GetPt()<<"\n";
37 // getchar();
dce1757b 38 Double_t Chi2Now, Chi2Ref;
39 if((*trackITS).GetLayer()==1 ) {
40 Chi2Now = trackITS->GetChi2();
41 Float_t NumClustNow = trackITS->GetNumClust();
42 if(trackITS->GetNumClust()) Chi2Now /= (Double_t )trackITS->GetNumClust();
43 Chi2Ref = reference->GetChi2();
44 Float_t NumClustRef = reference->GetNumClust();
45 if(reference->GetNumClust()) Chi2Ref /= (Double_t )reference->GetNumClust();
46 //cout<<" Chi2Now and Chi2Ref = "<<Chi2Now<<" "<<Chi2Ref<<"\n";
47 if( NumClustNow > NumClustRef ) {*reference = *trackITS;}
48 if((NumClustNow == NumClustRef )&& (Chi2Now < Chi2Ref)) {*reference = *trackITS;}
49 continue;
50 }
dce1757b 51 Float_t NumClustNow = trackITS->GetNumClust();
52 if(NumClustNow) {
53 Chi2Now = trackITS->GetChi2();
54 Chi2Now/=NumClustNow;
55 //cout<<" Chi2Now = "<<Chi2Now<<"\n";
56
57 // if(Ptref > 0.6 && Chi2Now > 20.) continue;
58 if(Ptref > 0.6 && Chi2Now > 30.) continue;
59 if((Ptref <= 0.6 && Ptref>0.2)&& Chi2Now > 15.) continue;
60 // if(Chi2Now>5.) continue;
61 //if(Chi2Now>15.) continue;
62 // if(Ptref <= 0.2 && Chi2Now > 10.) continue;
63 if(Ptref <= 0.2 && Chi2Now > 8.) continue;
64 }
65
66 Int_t layerInit = (*trackITS).GetLayer();
67 Int_t layernew = layerInit - 2; // -1 for new layer, -1 for matrix index
68
69 Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
70 Int_t NDetector[]= {4, 4, 5, 8, 23, 26};
71
72 TList listoftrack;
73 Int_t ladp, ladm, detp,detm,ladinters,detinters;
74 Int_t layerfin=layerInit-1;
75 Double_t Rfin=Rlayer[layerfin-1];
76 // cout<<"Prima di intersection \n";
cc6e75dd 77
dce1757b 78 Int_t outinters=NewIntersection(*trackITS, Rfin, layerfin, ladinters, detinters);
79
80 // cout<<" outinters = "<<outinters<<"\n";
81 // cout<<" Layer ladder detector intersection ="<<layerfin<<" "<<ladinters<<" "<<detinters<<"\n";
82 // cout << " phiinters zinters = "<<(*trackITS)(0) << " "<<(*trackITS)(1)<<"\n"; getchar();
83
84 if(outinters==-1) continue;
85
cc6e75dd 86 Int_t flaghit=0;
dce1757b 87 if(outinters==0){
88 TVector Touclad(9), Toucdet(9);
89 Int_t lycur=layerfin;
90 ladp=ladinters+1;
91 ladm=ladinters-1;
92 if(ladm <= 0) ladm=NLadder[layerfin-1];
93 if(ladp > NLadder[layerfin-1]) ladp=1;
94 detp=detinters+1;
95 detm=detinters-1;
96 Int_t idetot=1;
97 Touclad(0)=ladinters; Touclad(1)=ladm; Touclad(2)=ladp;
98 Touclad(3)=ladinters; Touclad(4)=ladm; Touclad(5)=ladp;
99 Touclad(6)=ladinters; Touclad(7)=ladm; Touclad(8)=ladp;
100 Toucdet(0)=detinters; Toucdet(1)=detinters; Toucdet(2)=detinters;
101 if(detm > 0 && detp <= NDetector[layerfin-1]) {
102 idetot=9;
103 Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
104 Toucdet(6)=detp; Toucdet(7)=detp; Toucdet(8)=detp;
105 }
106
107 if(detm > 0 && detp > NDetector[layerfin-1]) {
108 idetot=6;
109 Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
110 }
111
112 if(detm <= 0 && detp <= NDetector[layerfin-1]) {
113 idetot=6;
114 Toucdet(3)=detp; Toucdet(4)=detp; Toucdet(5)=detp;
cc6e75dd 115 }
116 Int_t iriv;
117 for (iriv=0; iriv<idetot; iriv++) { //for on detectors
dce1757b 118 AliITSgeom *g1 = aliITS->GetITSgeom();
119 TVector CTF(9);
120 g1->GetCenterThetaPhi(layerInit-1,(Int_t)Touclad(iriv),(Int_t)Toucdet(iriv),CTF);
121
122 // cout<<" layer, ladder, det, xo, yo, zo = "<<layerInit-1<<" "<<(Int_t)Touclad(iriv)<<
123 // " "<<(Int_t)Toucdet(iriv)<<" "<<CTF(0)<<" "<<CTF(1)<<" "<<CTF(2)<< " "<<CTF(6)<<"\n"; getchar();
124
125 ////////////////////////////////////////////////////////////////////////////////////////////////
126
127 /*** Rec points sorted by module *****/
128 /**************************************/
129
130 Int_t index;
131 AliITSRecPoint *recp;
132 AliITSgeom *geom = aliITS->GetITSgeom();
133 index = geom->GetModuleIndex(lycur,Touclad(iriv),Toucdet(iriv));
134 Int_t lay,lad,det;
135 geom->GetModuleId(index,lay,lad,det);
136 aliITS->ResetRecPoints();
137 gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty
138
139 Int_t npoints=rpoints->GetEntries();
140 Int_t *indlist=new Int_t[npoints+1];
141 Int_t counter=0;
142 Int_t ind;
143 for (ind=0; ind<=npoints; ind++) {
144 indlist[ind]=-1;
145 if (*(vettid[index]+ind)==0) {
146 indlist[counter]=ind;
147 counter++;
148 }
149 }
150
151 ind=-1;
cc6e75dd 152
dce1757b 153 for(;;) {
154 ind++;
155 if(indlist[ind] < 0) recp=0;
156 else recp = (AliITSRecPoint*)rpoints->UncheckedAt(indlist[ind]);
157
158 if((!recp) ) break;
159 TVector cluster(3),vecclust(9);
160 vecclust(6)=vecclust(7)=vecclust(8)=-1.;
161 Double_t sigma[2];
162 // set veclust in global
163 Float_t global[3], local[3];
164 local[0]=recp->GetX();
165 local[1]=0.;
166 local[2]= recp->GetZ();
167 AliITSgeom *g1 = aliITS->GetITSgeom();
168 Int_t play = lycur;
169 Int_t plad = TMath::Nint(Touclad(iriv));
170 Int_t pdet = TMath::Nint(Toucdet(iriv));
171 g1->LtoG(play,plad,pdet,local,global);
172
173 vecclust(0)=global[0];
174 vecclust(1)=global[1];
175 vecclust(2)=global[2];
176 vecclust(3) = (float)recp->fTracks[0];
177 vecclust(4) = (float)indlist[ind];
178 vecclust(5) = (float)index;
179 vecclust(6) = (float)recp->fTracks[0];
180 vecclust(7) = (float)recp->fTracks[1];
181 vecclust(8) = (float)recp->fTracks[2];
cc6e75dd 182
183 sigma[0] = (Double_t) recp->GetSigmaX2();
184 sigma[1] = (Double_t) recp->GetSigmaZ2();
dce1757b 185 //now we are in r,phi,z in global
186 cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit
187 cluster(1) = PhiDef(vecclust(0),vecclust(1)); // phi hit
188 cluster(2) = vecclust(2); // z hit
cc6e75dd 189 // cout<<" layer = "<<play<<"\n";
190 // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" "
191 // <<vecclust(2)<<"\n"; getchar();
dce1757b 192 //cluster(1)= cluster(1)-trackITS->Getalphaprov(); //provvisorio;
193 //if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi(); //provvisorio
194 //cout<<" cluster(1) dopo = "<<cluster(1)<< " alphaprov = "<<trackITS->Getalphaprov()<<"\n";
195 Float_t sigmatotphi, sigmatotz;
196
197 //Float_t epsphi=3.2, epsz=3.;
198 Float_t epsphi=3.2, epsz=3.0;
199 //if(Ptref<0.2) {epsphi=3.; epsz=3.;}
200
201 Double_t Rtrack=(*trackITS).Getrtrack();
202 Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack);
cc6e75dd 203 sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
204
dce1757b 205 sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
206 //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
cc6e75dd 207 //if(vecclust(3)==481) getchar();
dce1757b 208 if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi());
209 if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi());
210 if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue;
211 // cout<<" supero sigmaphi \n";
212 AliITStrack *newTrack = new AliITStrack((*trackITS));
213 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
214
215 if (TMath::Abs(Rtrack-cluster(0))/Rtrack>1e-6)
216 (*newTrack).Correct(Double_t(cluster(0)));
217 //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<< (*newTrack).GetZ()<<"\n";
218 if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){
219 delete newTrack;
220 continue;}
221
222 if(iriv == 0) flaghit=1;
cc6e75dd 223
224 (*newTrack).AddMS(); // add the multiple scattering matrix to the covariance matrix
dce1757b 225 (*newTrack).AddEL(1.,0);
cc6e75dd 226
dce1757b 227 Double_t sigmanew[2];
228 sigmanew[0]= sigmaphi;
229 sigmanew[1]=sigma[1];
230 //cout<<" Chiamo Kalman \n"; getchar();
cc6e75dd 231
dce1757b 232 if(flagvert)
233 KalmanFilterVert(newTrack,cluster,sigmanew);
234 else
235 KalmanFilter(newTrack,cluster,sigmanew);
cc6e75dd 236
dce1757b 237
dce1757b 238 (*newTrack).PutCluster(layernew, vecclust);
239 newTrack->AddClustInTrack();
240
241 listoftrack.AddLast(newTrack);
242
243 } // end of for(;;) on rec points
244
245 delete [] indlist;
246
247 } // end of for on detectors
cc6e75dd 248
dce1757b 249 }//end if(outinters==0)
cc6e75dd 250
dce1757b 251 if(flaghit==0 || outinters==-2) {
252 AliITStrack *newTrack = new AliITStrack(*trackITS);
253 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
254 (*newTrack).AddMS(); // add the multiple scattering matrix to the covariance matrix
255 (*newTrack).AddEL(1.,0);
256
257 listoftrack.AddLast(newTrack);
258 }
259
260
261 //gObjectTable->Print(); // stampa memoria
cc6e75dd 262
dce1757b 263 AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert);
264 listoftrack.Delete();
265 } // end of for on tracks
cc6e75dd 266
dce1757b 267 //gObjectTable->Print(); // stampa memoria
268
269}
270
271
272Int_t AliITStracking::NewIntersection(AliITStrack &track, Double_t rk,Int_t layer, Int_t &ladder, Int_t &detector) {
273//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
274// Found the intersection and the detector
275
276 if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;}
277 track.Propagation(rk);
278 Double_t zinters=track.GetZ();
279 Double_t phinters=track.Getphi();
280 //phinters = phinters+track.Getalphaprov(); //provvisorio
281 //if(phinters>2.*3.14) phinters=phinters-2.*3.14; //provvisorio
282 //cout<<"zinters = "<<zinters<<" phinters = "<<phinters<<"\n";
283
284 ////////////////////////////////// limits for Geometry 5 /////////////////////////////
285
286 Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
287 Int_t NDetector[]= {4, 4, 5, 8, 23, 26};
288
289 Float_t Detx[]= {0.64, 0.64, 3.509, 3.509, 3.65, 3.65 };
290 Float_t Detz[]= {4.19, 4.19, 3.75 , 3.75 , 2 , 2 };
291 ////////////////////////////////////////////////////////////////////////////////////////////////
292
293 TVector det(9);
294 TVector ListDet(2);
295 TVector DistZCenter(2);
296 AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
297 Int_t iz=0;
298 Double_t epsz=1.2;
299 Double_t epszpixel=0.05;
300
cc6e75dd 301 Int_t iD;
302 for(iD = 1; iD<= NDetector[layer-1]; iD++) {
dce1757b 303 g1->GetCenterThetaPhi(layer,1,iD,det);
304 Double_t zmin=det(2)-Detz[layer-1];
305 if(iD==1) zmin=det(2)-(Detz[layer-1])*epsz;
306 Double_t zmax=det(2)+Detz[layer-1];
307 if(iD==NDetector[layer-1]) zmax=det(2)+(Detz[layer-1])*epsz;
3d7c7ec8 308 //added to take into account problem on drift
dce1757b 309 if(layer == 3 || layer==2) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
310 //cout<<"zmin zinters zmax det(2)= "<<zmin<<" "<<zinters<<" "<<zmax<<" "<<det(2)<<"\n";
311 if(zinters > zmin && zinters <= zmax) {
312 if(iz>1) {cout<< " Errore su iz in NewIntersection \n"; getchar();}
313 else {
314 ListDet(iz)= iD; DistZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
315 }
316 }
317 }
318
319 if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;}
320 detector=Int_t (ListDet(0));
321 if(iz>1 && (DistZCenter(0)>DistZCenter(1))) detector=Int_t (ListDet(1));
322
323 AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
324 Float_t global[3];
325 Float_t local[3];
326 TVector ListLad(2);
327 TVector DistphiCenter(2);
328 Int_t ip=0;
329 Double_t pigre=TMath::Pi();
cc6e75dd 330
331 Int_t iLd;
332 for(iLd = 1; iLd<= NLadder[layer-1]; iLd++) {
dce1757b 333 g1->GetCenterThetaPhi(layer,iLd,detector,det);
334 Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));
dce1757b 335 // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
336 Double_t xmin,ymin,xmax,ymax;
3d7c7ec8 337 // Double_t phiconfr=0.0;
338 //cout<<" phiconfr inizio = "<<phiconfr <<"\n"; getchar();
dce1757b 339 local[1]=local[2]=0.;
340 local[0]= -(Detx[layer-1]);
341 if(layer==1) local[0]= (Detx[layer-1]); //take into account different reference system
342 g2->LtoG(layer,iLd,detector,local,global);
343 xmax=global[0]; ymax=global[1];
344 local[0]= (Detx[layer-1]);
345 if(layer==1) local[0]= -(Detx[layer-1]); //take into account different reference system
346 g2->LtoG(layer,iLd,detector,local,global);
347 xmin=global[0]; ymin=global[1];
348 Double_t phimin=PhiDef(xmin,ymin);
349 Double_t phimax=PhiDef(xmax,ymax);
350 //cout<<" xmin ymin = "<<xmin<<" "<<ymin<<"\n";
351 // cout<<" xmax ymax = "<<xmax<<" "<<ymax<<"\n";
3d7c7ec8 352 // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n";
353
354 Double_t phiconfr=phinters;
355 if(phimin>phimax ){
356 if(phimin <5.5) {cout<<" Error in NewIntersection for phi \n"; getchar();}
dce1757b 357 phimin=phimin-(2.*pigre);
3d7c7ec8 358 if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
dce1757b 359 if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
3d7c7ec8 360 }
361 // cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar();
dce1757b 362 if(phiconfr>phimin && phiconfr<= phimax) {
363 if(ip>1) {
364 cout<< " Errore su ip in NewIntersection \n"; getchar();
365 }
366 else {
367 ListLad(ip)= iLd; DistphiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
368 }
369 }
370 }
371 if(ip==0) { cout<< " No detector along phi \n"; getchar();}
372 ladder=Int_t (ListLad(0));
373 if(ip>1 && (DistphiCenter(0)>DistphiCenter(1))) ladder=Int_t (ListLad(1));
374
375 return 0;
376}
377
378
379Double_t AliITStracking::PhiDef(Double_t x, Double_t y){
380 Double_t pigre= TMath::Pi();
6204061f 381 Double_t phi=0.0;
dce1757b 382 if(y == 0. || x == 0.) {
383 if(y == 0. && x == 0.) {
384 cout << " Error in AliITStracking::PhiDef x=0 and y=0 \n"; getchar();
385 }
386 if(y==0. && x>0.) phi=0.;
387 if(y==0. && x<0.) phi=pigre;
388 if(x==0 && y>0.) phi=pigre/2.;
389 if(x==0 && y<0.) phi=1.5*pigre;
390 }
391 else {
392 if (x>0. && y>0.) phi=TMath::ATan(y/x);
393 if (x<0. && y>0.) phi=pigre+TMath::ATan(y/x);
394 if (x<0. && y<0.) phi=pigre+TMath::ATan(y/x);
395 if (x>0. && y<0.) phi=(2.*pigre)+TMath::ATan(y/x);
396 }
397 if(phi<0. || phi>(2*pigre)) {
398 cout<<" Error on phi in AliITStracking::PhiDef \n"; getchar();
399 }
400 return phi;
401}
402
403
404void AliITStracking::KalmanFilter(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
405//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
406// Kalman filter without vertex constraint
407
408 TMatrix H(2,5); H.UnitMatrix();
409 TMatrix Ht(TMatrix::kTransposed, H);
410
411
412 ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
413
414 TVector m(2);
415 Double_t rk,phik,zk;
416 rk=cluster(0); phik=cluster(1); zk=cluster(2);
417 m(0)=phik; m(1)=zk;
418 // cout<<" r and m = "<<rk<<" "<<m(0)<<" "<<m(1)<<"\n";
419 ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
420
421 TMatrix V(2,2);
422 V(0,1)=0.; V(1,0)=0.;
423 V(0,0)=sigma[0];
424 V(1,1)=sigma[1];
425 ///////////////////////////////////////////////////////////////////////////////////////////
426
427 TMatrix C=newTrack->GetCMatrix();
428 TMatrix tmp(H,TMatrix::kMult,C);
429 TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
430
431 R.Invert();
432 // cout<<" R prima = \n";
433 // R.Print(); getchar();
434
435 TMatrix K(C,TMatrix::kMult,Ht); K*=R;
436
437 TVector x=newTrack->GetVector();
438
439 TVector savex=x;
440 x*=H; x-=m;
441 x*=-1.; x*=K; x+=savex;
442 TMatrix saveC=C;
443 C.Mult(K,tmp); C-=saveC; C*=-1;
444 newTrack->GetVector()=x;
445 newTrack->GetCMatrix()=C;
446
447 TVector res= newTrack->GetVector();
448 //cout<<" res = "<<res(0)<<" "<<res(1)<<" "<<res(2)<<" "<<res(3)<<" "<<res(4)<<"\n";
449 res*=H; res-=m; res*=-1.;
450 TMatrix Cn=newTrack->GetCMatrix();
451 TMatrix tmpn(H,TMatrix::kMult,Cn);
452 TMatrix Rn(tmpn,TMatrix::kMult,Ht); Rn-=V; Rn*=-1.;
453
454 Rn.Invert();
455 TVector r=res; res*=Rn;
456 //cout<<" R dopo = \n";
457 //Rn.Print(); getchar();
458 Double_t chi2= r*res;
459 //cout<<"chi2 ="<<chi2<<"\n"; getchar();
460
461 newTrack->SetChi2(newTrack->GetChi2()+chi2);
462
463}
464
465
466void AliITStracking::KalmanFilterVert(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
467//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
468// Kalman filter with vertex constraint
469 TMatrix H(4,5); H.UnitMatrix();
470 TMatrix Ht(TMatrix::kTransposed, H);
471
472 ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
473
474 TVector m(4);
475 Double_t rk,phik,zk;
476 rk=cluster(0); phik=cluster(1); zk=cluster(2);
477 m(0)=phik; m(1)=zk;
478
479 Double_t CC=(*newTrack).GetC();
480 Double_t Zv=(*newTrack).GetZv();
481 Double_t Dv=(*newTrack).GetDv();
482 // cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
483 Double_t Cy=CC/2.;
484 Double_t tgl= (zk-Zv)*Cy/TMath::ASin(Cy*rk);
485 m(2)=Dv; m(3)=tgl;
486 // cout<<" m = \n";
487 // cout<<m(0)<<" "<<m(1)<<" "<<m(2)<<" "<<m(3)<<"\n";
488 ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
489
490 TMatrix V(4,4);
491 V(0,1)=0.; V(1,0)=0.;
492 Int_t Layer=newTrack->GetLayer();
493 V(0,0)=sigma[0];
494 V(1,1)=sigma[1];
495 V(1,3)=sigma[1]/rk; V(3,1)=V(1,3);
496 Double_t sigmaDv=newTrack->GetsigmaDv();
497 V(2,2)=sigmaDv*sigmaDv + newTrack->Getd2(Layer-1);
498 V(2,3)=newTrack->Getdtgl(Layer-1); V(3,2)=V(2,3);
499 Double_t sigmaZv=newTrack->GetsigmaZv();
500 V(3,3)=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
501 ///////////////////////////////////////////////////////////////////////////////////////////
502
503 //cout<<" d2 tgl2 dtgl = "<<newTrack->Getd2(Layer-1)<<" "<<newTrack->Gettgl2(Layer-1)<<" "<<newTrack->Getdtgl(Layer-1)<<"\n";
504 // cout<<" V = \n";
505 // V.Print(); getchar();
506
507 TMatrix C=newTrack->GetCMatrix();
cc6e75dd 508
509
510
dce1757b 511 TMatrix tmp(H,TMatrix::kMult,C);
512 TMatrix R(tmp,TMatrix::kMult,Ht); R+=V;
513
514 R.Invert();
515 TMatrix K(C,TMatrix::kMult,Ht); K*=R;
cc6e75dd 516 TVector x=newTrack->GetVector();
dce1757b 517 TVector savex=x;
518 x*=H; x-=m;
519 x*=-1; x*=K; x+=savex;
520 TMatrix saveC=C;
cc6e75dd 521 C.Mult(K,tmp); C-=saveC; C*=-1;
522
dce1757b 523 newTrack->GetVector()=x;
524 newTrack->GetCMatrix()=C;
525 TVector res= newTrack->GetVector();
526 //cout<<" res = "<<res(0)<<" "<<res(1)<<" "<<res(2)<<" "<<res(3)<<" "<<res(4)<<"\n";
527 res*=H; res-=m; res*=-1.;
528 TMatrix Cn=newTrack->GetCMatrix();
529 TMatrix tmpn(H,TMatrix::kMult,Cn);
530 TMatrix Rn(tmpn,TMatrix::kMult,Ht); Rn-=V; Rn*=-1.;
531
532 Rn.Invert();
533 TVector r=res; res*=Rn;
534 Double_t chi2= r*res;
535
536 newTrack->SetChi2(newTrack->GetChi2()+chi2);
537
538}
539