]> 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
The macros was working only in SSD mode. Now it works il All mode
[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.;
53e7090a 168 Double_t sigma[2];
169 //modificata 8-3-2001
dce1757b 170 // set veclust in global
171 Float_t global[3], local[3];
172 local[0]=recp->GetX();
173 local[1]=0.;
174 local[2]= recp->GetZ();
175 AliITSgeom *g1 = aliITS->GetITSgeom();
176 Int_t play = lycur;
177 Int_t plad = TMath::Nint(Touclad(iriv));
178 Int_t pdet = TMath::Nint(Toucdet(iriv));
179 g1->LtoG(play,plad,pdet,local,global);
53e7090a 180
dce1757b 181 vecclust(0)=global[0];
182 vecclust(1)=global[1];
53e7090a 183 vecclust(2)=global[2];
184
185 /*
186 //// modificato 8-3-2001
187 vecclust(0)=recp->GetRhit();;
188 vecclust(1)=recp->Getphi();
189 vecclust(2)=recp->GetZglobal();
190 */
191
dce1757b 192 vecclust(3) = (float)recp->fTracks[0];
193 vecclust(4) = (float)indlist[ind];
194 vecclust(5) = (float)index;
195 vecclust(6) = (float)recp->fTracks[0];
196 vecclust(7) = (float)recp->fTracks[1];
197 vecclust(8) = (float)recp->fTracks[2];
cc6e75dd 198
199 sigma[0] = (Double_t) recp->GetSigmaX2();
53e7090a 200 sigma[1] = (Double_t) recp->GetSigmaZ2();
201 //commentato 8-3-2001
dce1757b 202 //now we are in r,phi,z in global
203 cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit
204 cluster(1) = PhiDef(vecclust(0),vecclust(1)); // phi hit
53e7090a 205 cluster(2) = vecclust(2); // z hit
206
207 /*
208 //modificato 8-3-2001
209 cluster(0) = vecclust(0);//r hit
210 cluster(1) = vecclust(1); // phi hit
dce1757b 211 cluster(2) = vecclust(2); // z hit
53e7090a 212 */
cc6e75dd 213 // cout<<" layer = "<<play<<"\n";
214 // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" "
215 // <<vecclust(2)<<"\n"; getchar();
dce1757b 216 //cluster(1)= cluster(1)-trackITS->Getalphaprov(); //provvisorio;
217 //if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi(); //provvisorio
218 //cout<<" cluster(1) dopo = "<<cluster(1)<< " alphaprov = "<<trackITS->Getalphaprov()<<"\n";
219 Float_t sigmatotphi, sigmatotz;
220
221 //Float_t epsphi=3.2, epsz=3.;
8af13b4b 222 Float_t epsphi=5.0, epsz=5.0;
223 if(Ptref<0.2) {epsphi=3.; epsz=3.;}
dce1757b 224
225 Double_t Rtrack=(*trackITS).Getrtrack();
226 Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack);
cc6e75dd 227 sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
228
dce1757b 229 sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
230 //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
cc6e75dd 231 //if(vecclust(3)==481) getchar();
dce1757b 232 if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi());
233 if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi());
234 if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue;
235 // cout<<" supero sigmaphi \n";
236 AliITStrack *newTrack = new AliITStrack((*trackITS));
237 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
238
239 if (TMath::Abs(Rtrack-cluster(0))/Rtrack>1e-6)
240 (*newTrack).Correct(Double_t(cluster(0)));
241 //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<< (*newTrack).GetZ()<<"\n";
242 if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){
243 delete newTrack;
244 continue;}
245
246 if(iriv == 0) flaghit=1;
cc6e75dd 247
8af13b4b 248 (*newTrack).AddMS(rl); // add the multiple scattering matrix to the covariance matrix
249 (*newTrack).AddEL(rl,1.,0);
cc6e75dd 250
dce1757b 251 Double_t sigmanew[2];
252 sigmanew[0]= sigmaphi;
253 sigmanew[1]=sigma[1];
53e7090a 254
dce1757b 255 if(flagvert)
256 KalmanFilterVert(newTrack,cluster,sigmanew);
257 else
258 KalmanFilter(newTrack,cluster,sigmanew);
cc6e75dd 259
dce1757b 260
dce1757b 261 (*newTrack).PutCluster(layernew, vecclust);
262 newTrack->AddClustInTrack();
263
264 listoftrack.AddLast(newTrack);
265
266 } // end of for(;;) on rec points
267
268 delete [] indlist;
269
270 } // end of for on detectors
cc6e75dd 271
dce1757b 272 }//end if(outinters==0)
cc6e75dd 273
dce1757b 274 if(flaghit==0 || outinters==-2) {
275 AliITStrack *newTrack = new AliITStrack(*trackITS);
276 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
8af13b4b 277 (*newTrack).AddMS(rl); // add the multiple scattering matrix to the covariance matrix
278 (*newTrack).AddEL(rl,1.,0);
dce1757b 279
280 listoftrack.AddLast(newTrack);
281 }
282
283
284 //gObjectTable->Print(); // stampa memoria
53e7090a 285
8af13b4b 286 AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert,rl);
dce1757b 287 listoftrack.Delete();
288 } // end of for on tracks
cc6e75dd 289
dce1757b 290 //gObjectTable->Print(); // stampa memoria
291
292}
293
294
295Int_t AliITStracking::NewIntersection(AliITStrack &track, Double_t rk,Int_t layer, Int_t &ladder, Int_t &detector) {
296//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
297// Found the intersection and the detector
298
299 if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;}
300 track.Propagation(rk);
301 Double_t zinters=track.GetZ();
302 Double_t phinters=track.Getphi();
303 //phinters = phinters+track.Getalphaprov(); //provvisorio
304 //if(phinters>2.*3.14) phinters=phinters-2.*3.14; //provvisorio
305 //cout<<"zinters = "<<zinters<<" phinters = "<<phinters<<"\n";
306
307 ////////////////////////////////// limits for Geometry 5 /////////////////////////////
308
309 Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
8af13b4b 310 Int_t NDetector[]= {4, 4, 6, 8, 23, 26};
dce1757b 311
312 Float_t Detx[]= {0.64, 0.64, 3.509, 3.509, 3.65, 3.65 };
313 Float_t Detz[]= {4.19, 4.19, 3.75 , 3.75 , 2 , 2 };
314 ////////////////////////////////////////////////////////////////////////////////////////////////
315
316 TVector det(9);
317 TVector ListDet(2);
318 TVector DistZCenter(2);
319 AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
320 Int_t iz=0;
321 Double_t epsz=1.2;
322 Double_t epszpixel=0.05;
323
cc6e75dd 324 Int_t iD;
325 for(iD = 1; iD<= NDetector[layer-1]; iD++) {
dce1757b 326 g1->GetCenterThetaPhi(layer,1,iD,det);
327 Double_t zmin=det(2)-Detz[layer-1];
328 if(iD==1) zmin=det(2)-(Detz[layer-1])*epsz;
329 Double_t zmax=det(2)+Detz[layer-1];
330 if(iD==NDetector[layer-1]) zmax=det(2)+(Detz[layer-1])*epsz;
3d7c7ec8 331 //added to take into account problem on drift
99b69e4c 332 if(layer == 4 || layer==3) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
dce1757b 333 //cout<<"zmin zinters zmax det(2)= "<<zmin<<" "<<zinters<<" "<<zmax<<" "<<det(2)<<"\n";
334 if(zinters > zmin && zinters <= zmax) {
335 if(iz>1) {cout<< " Errore su iz in NewIntersection \n"; getchar();}
336 else {
337 ListDet(iz)= iD; DistZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
338 }
339 }
340 }
341
342 if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;}
343 detector=Int_t (ListDet(0));
344 if(iz>1 && (DistZCenter(0)>DistZCenter(1))) detector=Int_t (ListDet(1));
345
346 AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
347 Float_t global[3];
348 Float_t local[3];
349 TVector ListLad(2);
350 TVector DistphiCenter(2);
351 Int_t ip=0;
352 Double_t pigre=TMath::Pi();
cc6e75dd 353
354 Int_t iLd;
355 for(iLd = 1; iLd<= NLadder[layer-1]; iLd++) {
dce1757b 356 g1->GetCenterThetaPhi(layer,iLd,detector,det);
357 Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));
dce1757b 358 // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
359 Double_t xmin,ymin,xmax,ymax;
3d7c7ec8 360 // Double_t phiconfr=0.0;
361 //cout<<" phiconfr inizio = "<<phiconfr <<"\n"; getchar();
dce1757b 362 local[1]=local[2]=0.;
363 local[0]= -(Detx[layer-1]);
364 if(layer==1) local[0]= (Detx[layer-1]); //take into account different reference system
365 g2->LtoG(layer,iLd,detector,local,global);
366 xmax=global[0]; ymax=global[1];
367 local[0]= (Detx[layer-1]);
368 if(layer==1) local[0]= -(Detx[layer-1]); //take into account different reference system
369 g2->LtoG(layer,iLd,detector,local,global);
370 xmin=global[0]; ymin=global[1];
371 Double_t phimin=PhiDef(xmin,ymin);
372 Double_t phimax=PhiDef(xmax,ymax);
373 //cout<<" xmin ymin = "<<xmin<<" "<<ymin<<"\n";
374 // cout<<" xmax ymax = "<<xmax<<" "<<ymax<<"\n";
3d7c7ec8 375 // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n";
376
377 Double_t phiconfr=phinters;
378 if(phimin>phimax ){
379 if(phimin <5.5) {cout<<" Error in NewIntersection for phi \n"; getchar();}
dce1757b 380 phimin=phimin-(2.*pigre);
3d7c7ec8 381 if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
dce1757b 382 if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
3d7c7ec8 383 }
384 // cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar();
dce1757b 385 if(phiconfr>phimin && phiconfr<= phimax) {
386 if(ip>1) {
387 cout<< " Errore su ip in NewIntersection \n"; getchar();
388 }
389 else {
390 ListLad(ip)= iLd; DistphiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
391 }
392 }
393 }
394 if(ip==0) { cout<< " No detector along phi \n"; getchar();}
395 ladder=Int_t (ListLad(0));
396 if(ip>1 && (DistphiCenter(0)>DistphiCenter(1))) ladder=Int_t (ListLad(1));
397
398 return 0;
399}
400
401
402Double_t AliITStracking::PhiDef(Double_t x, Double_t y){
403 Double_t pigre= TMath::Pi();
6204061f 404 Double_t phi=0.0;
dce1757b 405 if(y == 0. || x == 0.) {
406 if(y == 0. && x == 0.) {
407 cout << " Error in AliITStracking::PhiDef x=0 and y=0 \n"; getchar();
408 }
409 if(y==0. && x>0.) phi=0.;
410 if(y==0. && x<0.) phi=pigre;
411 if(x==0 && y>0.) phi=pigre/2.;
412 if(x==0 && y<0.) phi=1.5*pigre;
413 }
414 else {
415 if (x>0. && y>0.) phi=TMath::ATan(y/x);
416 if (x<0. && y>0.) phi=pigre+TMath::ATan(y/x);
417 if (x<0. && y<0.) phi=pigre+TMath::ATan(y/x);
418 if (x>0. && y<0.) phi=(2.*pigre)+TMath::ATan(y/x);
419 }
420 if(phi<0. || phi>(2*pigre)) {
421 cout<<" Error on phi in AliITStracking::PhiDef \n"; getchar();
422 }
423 return phi;
424}
425
426
427void AliITStracking::KalmanFilter(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
428//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
429// Kalman filter without vertex constraint
430
dce1757b 431
432 ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
433
53e7090a 434 Double_t m[2];
dce1757b 435 Double_t rk,phik,zk;
436 rk=cluster(0); phik=cluster(1); zk=cluster(2);
53e7090a 437 m[0]=phik; m[1]=zk;
438
dce1757b 439 ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
440
53e7090a 441 Double_t V00=sigma[0];
442 Double_t V11=sigma[1];
443
dce1757b 444 ///////////////////////////////////////////////////////////////////////////////////////////
445
53e7090a 446
447 Double_t Cin00,Cin10,Cin20,Cin30,Cin40,Cin11,Cin21,Cin31,Cin41,Cin22,Cin32,Cin42,Cin33,Cin43,Cin44;
448
449 newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,
450 Cin41,Cin42,Cin43,Cin44); //get C matrix
451
452 Double_t Rold00=Cin00+V00;
453 Double_t Rold10=Cin10;
454 Double_t Rold11=Cin11+V11;
455
456//////////////////////////////////// R matrix inversion ///////////////////////////////////////////////
457
458 Double_t det=Rold00*Rold11-Rold10*Rold10;
459 Double_t R00=Rold11/det;
460 Double_t R10=-Rold10/det;
461 Double_t R11=Rold00/det;
462
463////////////////////////////////////////////////////////////////////////////////////////////////////////
464
465 Double_t K00=Cin00*R00+Cin10*R10;
466 Double_t K01=Cin00*R10+Cin10*R11;
467 Double_t K10=Cin10*R00+Cin11*R10;
468 Double_t K11=Cin10*R10+Cin11*R11;
469 Double_t K20=Cin20*R00+Cin21*R10;
470 Double_t K21=Cin20*R10+Cin21*R11;
471 Double_t K30=Cin30*R00+Cin31*R10;
472 Double_t K31=Cin30*R10+Cin31*R11;
473 Double_t K40=Cin40*R00+Cin41*R10;
474 Double_t K41=Cin40*R10+Cin41*R11;
475
476 Double_t X0,X1,X2,X3,X4;
477 newTrack->GetXElements(X0,X1,X2,X3,X4); // get the state vector
478
479 Double_t savex0=X0, savex1=X1;
480
481 X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1);
482 X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1);
483 X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1);
484 X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1);
485 X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1);
486
487 Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44;
488
489 C00=Cin00-K00*Cin00-K01*Cin10;
490 C10=Cin10-K00*Cin10-K01*Cin11;
491 C20=Cin20-K00*Cin20-K01*Cin21;
492 C30=Cin30-K00*Cin30-K01*Cin31;
493 C40=Cin40-K00*Cin40-K01*Cin41;
494
495 C11=Cin11-K10*Cin10-K11*Cin11;
496 C21=Cin21-K10*Cin20-K11*Cin21;
497 C31=Cin31-K10*Cin30-K11*Cin31;
498 C41=Cin41-K10*Cin40-K11*Cin41;
499
500 C22=Cin22-K20*Cin20-K21*Cin21;
501 C32=Cin32-K20*Cin30-K21*Cin31;
502 C42=Cin42-K20*Cin40-K21*Cin41;
503
504 C33=Cin33-K30*Cin30-K31*Cin31;
505 C43=Cin43-K30*Cin40-K31*Cin41;
506
507 C44=Cin44-K40*Cin40-K41*Cin41;
508
509 newTrack->PutXElements(X0,X1,X2,X3,X4); // put the new state vector
510
511 newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the
512 // new cov matrix
513 Double_t VMCold00=V00-C00;
514 Double_t VMCold10=-C10;
515 Double_t VMCold11=V11-C11;
516
517///////////////////////////////////// Matrix VMC inversion ////////////////////////////////////////////////
518
519 det=VMCold00*VMCold11-VMCold10*VMCold10;
520 Double_t VMC00=VMCold11/det;
521 Double_t VMC10=-VMCold10/det;
522 Double_t VMC11=VMCold00/det;
523
524////////////////////////////////////////////////////////////////////////////////////////////////////////////
525
526 Double_t chi2=(m[0]-X0)*( VMC00*(m[0]-X0) + 2.*VMC10*(m[1]-X1) ) +
527 (m[1]-X1)*VMC11*(m[1]-X1);
dce1757b 528
529 newTrack->SetChi2(newTrack->GetChi2()+chi2);
530
531}
532
533
534void AliITStracking::KalmanFilterVert(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
535//Origin A. Badala' and G.S. Pappalardo: e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
536// Kalman filter with vertex constraint
dce1757b 537
53e7090a 538 ////////////////////////////// Evaluation of the measurement vector m ///////////////
dce1757b 539
53e7090a 540 Double_t m[4];
dce1757b 541 Double_t rk,phik,zk;
542 rk=cluster(0); phik=cluster(1); zk=cluster(2);
53e7090a 543 m[0]=phik; m[1]=zk;
dce1757b 544
545 Double_t CC=(*newTrack).GetC();
546 Double_t Zv=(*newTrack).GetZv();
547 Double_t Dv=(*newTrack).GetDv();
dce1757b 548 Double_t Cy=CC/2.;
549 Double_t tgl= (zk-Zv)*Cy/TMath::ASin(Cy*rk);
53e7090a 550 m[2]=Dv; m[3]=tgl;
dce1757b 551
53e7090a 552 ///////////////////////////////////// Evaluation of the error matrix V //////////////
dce1757b 553 Int_t Layer=newTrack->GetLayer();
53e7090a 554 Double_t V00=sigma[0];
555 Double_t V11=sigma[1];
556 Double_t V31=sigma[1]/rk;
dce1757b 557 Double_t sigmaDv=newTrack->GetsigmaDv();
53e7090a 558 Double_t V22=sigmaDv*sigmaDv + newTrack->Getd2(Layer-1);
559 Double_t V32=newTrack->Getdtgl(Layer-1);
dce1757b 560 Double_t sigmaZv=newTrack->GetsigmaZv();
53e7090a 561 Double_t V33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
562 ///////////////////////////////////////////////////////////////////////////////////////
563
564 Double_t Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,Cin41,Cin42,Cin43,Cin44;
565
566 newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,
567 Cin41,Cin42,Cin43,Cin44); //get C matrix
568
569 Double_t R[4][4];
570 R[0][0]=Cin00+V00;
571 R[1][0]=Cin10;
572 R[2][0]=Cin20;
573 R[3][0]=Cin30;
574 R[1][1]=Cin11+V11;
575 R[2][1]=Cin21;
576 R[3][1]=Cin31+sigma[1]/rk;
577 R[2][2]=Cin22+sigmaDv*sigmaDv+newTrack->Getd2(Layer-1);
578 R[3][2]=Cin32+newTrack->Getdtgl(Layer-1);
579 R[3][3]=Cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
580
581 R[0][1]=R[1][0]; R[0][2]=R[2][0]; R[0][3]=R[3][0]; R[1][2]=R[2][1]; R[1][3]=R[3][1];
582 R[2][3]=R[3][2];
583
584///////////////////// Matrix R inversion ////////////////////////////////////////////
585
03454f1c 586 const Int_t n=4;
53e7090a 587 Double_t big, hold;
588 Double_t d=1.;
589 Int_t ll[n],mm[n];
03454f1c 590
591 Int_t i,j,k;
53e7090a 592
03454f1c 593 for(k=0; k<n; k++) {
53e7090a 594 ll[k]=k;
595 mm[k]=k;
596 big=R[k][k];
03454f1c 597 for(j=k; j<n ; j++) {
598 for (i=j; i<n; i++) {
53e7090a 599 if(TMath::Abs(big) < TMath::Abs(R[i][j]) ) { big=R[i][j]; ll[k]=i; mm[k]=j; }
600 }
601 }
602//
03454f1c 603 j= ll[k];
53e7090a 604 if(j > k) {
03454f1c 605 for(i=0; i<n; i++) { hold=-R[k][i]; R[k][i]=R[j][i]; R[j][i]=hold; }
53e7090a 606
607 }
608//
03454f1c 609 i=mm[k];
53e7090a 610 if(i > k ) {
03454f1c 611 for(j=0; j<n; j++) { hold=-R[j][k]; R[j][k]=R[j][i]; R[j][i]=hold; }
53e7090a 612 }
613//
614 if(!big) {
615 d=0.;
616 cout << "Singular matrix\n";
617 }
03454f1c 618 for(i=0; i<n; i++) {
53e7090a 619 if(i == k) { continue; }
620 R[i][k]=R[i][k]/(-big);
621 }
622//
03454f1c 623 for(i=0; i<n; i++) {
53e7090a 624 hold=R[i][k];
03454f1c 625 for(j=0; j<n; j++) {
53e7090a 626 if(i == k || j == k) { continue; }
627 R[i][j]=hold*R[k][j]+R[i][j];
628 }
629 }
630//
03454f1c 631 for(j=0; j<n; j++) {
53e7090a 632 if(j == k) { continue; }
633 R[k][j]=R[k][j]/big;
634 }
635//
636 d=d*big;
637//
638 R[k][k]=1./big;
639 }
640//
03454f1c 641 for(k=n-1; k>=0; k--) {
642 i=ll[k];
53e7090a 643 if(i > k) {
03454f1c 644 for (j=0; j<n; j++) {hold=R[j][k]; R[j][k]=-R[j][i]; R[j][i]=hold;}
53e7090a 645 }
03454f1c 646 j=mm[k];
53e7090a 647 if(j > k) {
03454f1c 648 for (i=0; i<n; i++) {hold=R[k][i]; R[k][i]=-R[j][i]; R[j][i]=hold;}
53e7090a 649 }
650 }
651//////////////////////////////////////////////////////////////////////////////////
652
653
654 Double_t K00=Cin00*R[0][0]+Cin10*R[1][0]+Cin20*R[2][0]+Cin30*R[3][0];
655 Double_t K01=Cin00*R[1][0]+Cin10*R[1][1]+Cin20*R[2][1]+Cin30*R[3][1];
656 Double_t K02=Cin00*R[2][0]+Cin10*R[2][1]+Cin20*R[2][2]+Cin30*R[3][2];
657 Double_t K03=Cin00*R[3][0]+Cin10*R[3][1]+Cin20*R[3][2]+Cin30*R[3][3];
658 Double_t K10=Cin10*R[0][0]+Cin11*R[1][0]+Cin21*R[2][0]+Cin31*R[3][0];
659 Double_t K11=Cin10*R[1][0]+Cin11*R[1][1]+Cin21*R[2][1]+Cin31*R[3][1];
660 Double_t K12=Cin10*R[2][0]+Cin11*R[2][1]+Cin21*R[2][2]+Cin31*R[3][2];
661 Double_t K13=Cin10*R[3][0]+Cin11*R[3][1]+Cin21*R[3][2]+Cin31*R[3][3];
662 Double_t K20=Cin20*R[0][0]+Cin21*R[1][0]+Cin22*R[2][0]+Cin32*R[3][0];
663 Double_t K21=Cin20*R[1][0]+Cin21*R[1][1]+Cin22*R[2][1]+Cin32*R[3][1];
664 Double_t K22=Cin20*R[2][0]+Cin21*R[2][1]+Cin22*R[2][2]+Cin32*R[3][2];
665 Double_t K23=Cin20*R[3][0]+Cin21*R[3][1]+Cin22*R[3][2]+Cin32*R[3][3];
666 Double_t K30=Cin30*R[0][0]+Cin31*R[1][0]+Cin32*R[2][0]+Cin33*R[3][0];
667 Double_t K31=Cin30*R[1][0]+Cin31*R[1][1]+Cin32*R[2][1]+Cin33*R[3][1];
668 Double_t K32=Cin30*R[2][0]+Cin31*R[2][1]+Cin32*R[2][2]+Cin33*R[3][2];
669 Double_t K33=Cin30*R[3][0]+Cin31*R[3][1]+Cin32*R[3][2]+Cin33*R[3][3];
670 Double_t K40=Cin40*R[0][0]+Cin41*R[1][0]+Cin42*R[2][0]+Cin43*R[3][0];
671 Double_t K41=Cin40*R[1][0]+Cin41*R[1][1]+Cin42*R[2][1]+Cin43*R[3][1];
672 Double_t K42=Cin40*R[2][0]+Cin41*R[2][1]+Cin42*R[2][2]+Cin43*R[3][2];
673 Double_t K43=Cin40*R[3][0]+Cin41*R[3][1]+Cin42*R[3][2]+Cin43*R[3][3];
674
675 Double_t X0,X1,X2,X3,X4;
676 newTrack->GetXElements(X0,X1,X2,X3,X4); // get the state vector
677
678 Double_t savex0=X0, savex1=X1, savex2=X2, savex3=X3;
679
680 X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1)+K02*(m[2]-savex2)+
681 K03*(m[3]-savex3);
682 X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1)+K12*(m[2]-savex2)+
683 K13*(m[3]-savex3);
684 X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1)+K22*(m[2]-savex2)+
685 K23*(m[3]-savex3);
686 X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1)+K32*(m[2]-savex2)+
687 K33*(m[3]-savex3);
688 X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1)+K42*(m[2]-savex2)+
689 K43*(m[3]-savex3);
690
691 Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44;
692
693 C00=Cin00-K00*Cin00-K01*Cin10-K02*Cin20-K03*Cin30;
694 C10=Cin10-K00*Cin10-K01*Cin11-K02*Cin21-K03*Cin31;
695 C20=Cin20-K00*Cin20-K01*Cin21-K02*Cin22-K03*Cin32;
696 C30=Cin30-K00*Cin30-K01*Cin31-K02*Cin32-K03*Cin33;
697 C40=Cin40-K00*Cin40-K01*Cin41-K02*Cin42-K03*Cin43;
698
699 C11=Cin11-K10*Cin10-K11*Cin11-K12*Cin21-K13*Cin31;
700 C21=Cin21-K10*Cin20-K11*Cin21-K12*Cin22-K13*Cin32;
701 C31=Cin31-K10*Cin30-K11*Cin31-K12*Cin32-K13*Cin33;
702 C41=Cin41-K10*Cin40-K11*Cin41-K12*Cin42-K13*Cin43;
703
704 C22=Cin22-K20*Cin20-K21*Cin21-K22*Cin22-K23*Cin32;
705 C32=Cin32-K20*Cin30-K21*Cin31-K22*Cin32-K23*Cin33;
706 C42=Cin42-K20*Cin40-K21*Cin41-K22*Cin42-K23*Cin43;
707
708 C33=Cin33-K30*Cin30-K31*Cin31-K32*Cin32-K33*Cin33;
709 C43=Cin43-K30*Cin40-K31*Cin41-K32*Cin42-K33*Cin43;
710
711 C44=Cin44-K40*Cin40-K41*Cin41-K42*Cin42-K43*Cin43;
712
713 newTrack->PutXElements(X0,X1,X2,X3,X4); // put the new state vector
714
715 newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the
716 // new cov matrix
717
718 Double_t VMC[4][4];
719
720 VMC[0][0]=V00-C00; VMC[1][0]=-C10; VMC[2][0]=-C20; VMC[3][0]=-C30;
721 VMC[1][1]=V11-C11; VMC[2][1]=-C21; VMC[3][1]=V31-C31;
722 VMC[2][2]=V22-C22; VMC[3][2]=V32-C32;
723 VMC[3][3]=V33-C33;
724
725 VMC[0][1]=VMC[1][0]; VMC[0][2]=VMC[2][0]; VMC[0][3]=VMC[3][0];
726 VMC[1][2]=VMC[2][1]; VMC[1][3]=VMC[3][1];
727 VMC[2][3]=VMC[3][2];
728
729
730/////////////////////// VMC matrix inversion ///////////////////////////////////
dce1757b 731
53e7090a 732 d=1.;
733
03454f1c 734 for(k=0; k<n; k++) {
53e7090a 735 ll[k]=k;
736 mm[k]=k;
737 big=VMC[k][k];
03454f1c 738 for(j=k; j<n ; j++) {
739 for (i=j; i<n; i++) {
53e7090a 740 if(TMath::Abs(big) < TMath::Abs(VMC[i][j]) ) { big=VMC[i][j]; ll[k]=i; mm[k]=j; }
741 }
742 }
743//
03454f1c 744 j= ll[k];
53e7090a 745 if(j > k) {
03454f1c 746 for(i=0; i<n; i++) { hold=-VMC[k][i]; VMC[k][i]=VMC[j][i]; VMC[j][i]=hold; }
53e7090a 747
748 }
749//
03454f1c 750 i=mm[k];
53e7090a 751 if(i > k ) {
03454f1c 752 for(j=0; j<n; j++) { hold=-VMC[j][k]; VMC[j][k]=VMC[j][i]; VMC[j][i]=hold; }
53e7090a 753 }
754//
755 if(!big) {
756 d=0.;
757 cout << "Singular matrix\n";
758 }
03454f1c 759 for(i=0; i<n; i++) {
53e7090a 760 if(i == k) { continue; }
761 VMC[i][k]=VMC[i][k]/(-big);
762 }
763//
03454f1c 764 for(i=0; i<n; i++) {
53e7090a 765 hold=VMC[i][k];
03454f1c 766 for(j=0; j<n; j++) {
53e7090a 767 if(i == k || j == k) { continue; }
768 VMC[i][j]=hold*VMC[k][j]+VMC[i][j];
769 }
770 }
771//
03454f1c 772 for(j=0; j<n; j++) {
53e7090a 773 if(j == k) { continue; }
774 VMC[k][j]=VMC[k][j]/big;
775 }
776//
777 d=d*big;
778//
779 VMC[k][k]=1./big;
780 }
781//
03454f1c 782 for(k=n-1; k>=0; k--) {
783 i=ll[k];
53e7090a 784 if(i > k) {
03454f1c 785 for (j=0; j<n; j++) {hold=VMC[j][k]; VMC[j][k]=-VMC[j][i]; VMC[j][i]=hold;}
53e7090a 786 }
03454f1c 787 j=mm[k];
53e7090a 788 if(j > k) {
03454f1c 789 for (i=0; i<n; i++) {hold=VMC[k][i]; VMC[k][i]=-VMC[j][i]; VMC[j][i]=hold;}
53e7090a 790 }
791 }
792
793
794////////////////////////////////////////////////////////////////////////////////
795
796 Double_t chi2=(m[0]-X0)*( VMC[0][0]*(m[0]-X0) + 2.*VMC[1][0]*(m[1]-X1) +
797 2.*VMC[2][0]*(m[2]-X2)+ 2.*VMC[3][0]*(m[3]-X3) ) +
798 (m[1]-X1)* ( VMC[1][1]*(m[1]-X1) + 2.*VMC[2][1]*(m[2]-X2)+
799 2.*VMC[3][1]*(m[3]-X3) ) +
800 (m[2]-X2)* ( VMC[2][2]*(m[2]-X2)+ 2.*VMC[3][2]*(m[3]-X3) ) +
801 (m[3]-X3)*VMC[3][3]*(m[3]-X3);
dce1757b 802
803 newTrack->SetChi2(newTrack->GetChi2()+chi2);
804
805}
806