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