7 #include <TObjectTable.h>
9 #include "AliITStracking.h"
12 #include "AliITSgeom.h"
13 #include "AliITSRecPoint.h"
14 #include "AliITStrack.h"
16 ClassImp(AliITStracking)
19 AliITStracking::AliITStracking(TList *trackITSlist, AliITStrack *reference,
20 AliITS *aliITS, TObjArray *rpoints, Double_t Ptref, Int_t **vettid, Bool_t flagvert, AliITSRad *rl) {
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
26 Rlayer[0]=4.; Rlayer[1]=7.; Rlayer[2]=14.9; Rlayer[3]=23.8; Rlayer[4]=39.1; Rlayer[5]=43.6;
29 for(index =0; index<trackITSlist->GetSize(); index++) {
30 AliITStrack *trackITS = (AliITStrack *) trackITSlist->At(index);
32 if((*trackITS).GetLayer()==7) reference->SetChi2(10.223e140);
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";
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;}
52 Float_t NumClustNow = trackITS->GetNumClust();
54 Chi2Now = trackITS->GetChi2();
56 //cout<<" Chi2Now = "<<Chi2Now<<"\n";
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;
64 if(Ptref <= 0.2 && Chi2Now > 8.) continue;
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;
72 Int_t layerInit = (*trackITS).GetLayer();
73 Int_t layernew = layerInit - 2; // -1 for new layer, -1 for matrix index
75 Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
76 Int_t NDetector[]= {4, 4, 6, 8, 23, 26};
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";
84 Int_t outinters=NewIntersection(*trackITS, Rfin, layerfin, ladinters, detinters);
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();
90 if(outinters==-1) continue;
94 TVector Touclad(9), Toucdet(9);
98 if(ladm <= 0) ladm=NLadder[layerfin-1];
99 if(ladp > NLadder[layerfin-1]) ladp=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]) {
109 Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
110 Toucdet(6)=detp; Toucdet(7)=detp; Toucdet(8)=detp;
113 if(detm > 0 && detp > NDetector[layerfin-1]) {
115 Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
118 if(detm <= 0 && detp <= NDetector[layerfin-1]) {
120 Toucdet(3)=detp; Toucdet(4)=detp; Toucdet(5)=detp;
123 for (iriv=0; iriv<idetot; iriv++) { //for on detectors
124 AliITSgeom *g1 = aliITS->GetITSgeom();
126 g1->GetCenterThetaPhi(layerInit-1,(Int_t)Touclad(iriv),(Int_t)Toucdet(iriv),CTF);
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();
131 ////////////////////////////////////////////////////////////////////////////////////////////////
133 /*** Rec points sorted by module *****/
134 /**************************************/
137 AliITSRecPoint *recp;
138 AliITSgeom *geom = aliITS->GetITSgeom();
139 index = geom->GetModuleIndex(lycur,Touclad(iriv),Toucdet(iriv));
141 geom->GetModuleId(index,lay,lad,det);
142 aliITS->ResetRecPoints();
143 //gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty
144 gAlice->TreeR()->GetEvent(index); //first entry in TreeR is empty
146 Int_t npoints=rpoints->GetEntries();
147 Int_t *indlist=new Int_t[npoints+1];
150 for (ind=0; ind<=npoints; ind++) {
152 if (*(vettid[index]+ind)==0) {
153 indlist[counter]=ind;
162 if(indlist[ind] < 0) recp=0;
163 else recp = (AliITSRecPoint*)rpoints->UncheckedAt(indlist[ind]);
166 TVector cluster(3),vecclust(9);
167 vecclust(6)=vecclust(7)=vecclust(8)=-1.;
169 //modificata 8-3-2001
170 // set veclust in global
171 Float_t global[3], local[3];
172 local[0]=recp->GetX();
174 local[2]= recp->GetZ();
175 AliITSgeom *g1 = aliITS->GetITSgeom();
177 Int_t plad = TMath::Nint(Touclad(iriv));
178 Int_t pdet = TMath::Nint(Toucdet(iriv));
179 g1->LtoG(play,plad,pdet,local,global);
181 vecclust(0)=global[0];
182 vecclust(1)=global[1];
183 vecclust(2)=global[2];
186 //// modificato 8-3-2001
187 vecclust(0)=recp->GetRhit();;
188 vecclust(1)=recp->Getphi();
189 vecclust(2)=recp->GetZglobal();
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];
199 sigma[0] = (Double_t) recp->GetSigmaX2();
200 sigma[1] = (Double_t) recp->GetSigmaZ2();
201 //commentato 8-3-2001
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
205 cluster(2) = vecclust(2); // z hit
208 //modificato 8-3-2001
209 cluster(0) = vecclust(0);//r hit
210 cluster(1) = vecclust(1); // phi hit
211 cluster(2) = vecclust(2); // z hit
213 // cout<<" layer = "<<play<<"\n";
214 // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" "
215 // <<vecclust(2)<<"\n"; getchar();
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;
221 //Float_t epsphi=3.2, epsz=3.;
222 Float_t epsphi=5.0, epsz=5.0;
223 if(Ptref<0.2) {epsphi=3.; epsz=3.;}
225 Double_t Rtrack=(*trackITS).Getrtrack();
226 Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack);
227 sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
229 sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
230 //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
231 //if(vecclust(3)==481) getchar();
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);
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){
246 if(iriv == 0) flaghit=1;
248 (*newTrack).AddMS(rl); // add the multiple scattering matrix to the covariance matrix
249 (*newTrack).AddEL(rl,1.,0);
251 Double_t sigmanew[2];
252 sigmanew[0]= sigmaphi;
253 sigmanew[1]=sigma[1];
256 KalmanFilterVert(newTrack,cluster,sigmanew);
258 KalmanFilter(newTrack,cluster,sigmanew);
261 (*newTrack).PutCluster(layernew, vecclust);
262 newTrack->AddClustInTrack();
264 listoftrack.AddLast(newTrack);
266 } // end of for(;;) on rec points
270 } // end of for on detectors
272 }//end if(outinters==0)
274 if(flaghit==0 || outinters==-2) {
275 AliITStrack *newTrack = new AliITStrack(*trackITS);
276 (*newTrack).SetLayer((*trackITS).GetLayer()-1);
277 (*newTrack).AddMS(rl); // add the multiple scattering matrix to the covariance matrix
278 (*newTrack).AddEL(rl,1.,0);
280 listoftrack.AddLast(newTrack);
284 //gObjectTable->Print(); // stampa memoria
286 AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert,rl);
287 listoftrack.Delete();
288 } // end of for on tracks
290 //gObjectTable->Print(); // stampa memoria
295 Int_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
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";
307 ////////////////////////////////// limits for Geometry 5 /////////////////////////////
309 Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
310 Int_t NDetector[]= {4, 4, 6, 8, 23, 26};
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 ////////////////////////////////////////////////////////////////////////////////////////////////
318 TVector DistZCenter(2);
319 AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
322 Double_t epszpixel=0.05;
325 for(iD = 1; iD<= NDetector[layer-1]; iD++) {
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;
331 //added to take into account problem on drift
332 if(layer == 4 || layer==3) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
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();}
337 ListDet(iz)= iD; DistZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
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));
346 AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
350 TVector DistphiCenter(2);
352 Double_t pigre=TMath::Pi();
355 for(iLd = 1; iLd<= NLadder[layer-1]; iLd++) {
356 g1->GetCenterThetaPhi(layer,iLd,detector,det);
357 Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));
358 // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
359 Double_t xmin,ymin,xmax,ymax;
360 // Double_t phiconfr=0.0;
361 //cout<<" phiconfr inizio = "<<phiconfr <<"\n"; getchar();
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";
375 // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n";
377 Double_t phiconfr=phinters;
379 if(phimin <5.5) {cout<<" Error in NewIntersection for phi \n"; getchar();}
380 phimin=phimin-(2.*pigre);
381 if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
382 if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
384 // cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar();
385 if(phiconfr>phimin && phiconfr<= phimax) {
387 cout<< " Errore su ip in NewIntersection \n"; getchar();
390 ListLad(ip)= iLd; DistphiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
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));
402 Double_t AliITStracking::PhiDef(Double_t x, Double_t y){
403 Double_t pigre= TMath::Pi();
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();
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;
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);
420 if(phi<0. || phi>(2*pigre)) {
421 cout<<" Error on phi in AliITStracking::PhiDef \n"; getchar();
427 void 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
432 ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
436 rk=cluster(0); phik=cluster(1); zk=cluster(2);
439 ///////////////////////////////////// Evaluation of the error matrix V ///////////////////////////////
441 Double_t V00=sigma[0];
442 Double_t V11=sigma[1];
444 ///////////////////////////////////////////////////////////////////////////////////////////
447 Double_t Cin00,Cin10,Cin20,Cin30,Cin40,Cin11,Cin21,Cin31,Cin41,Cin22,Cin32,Cin42,Cin33,Cin43,Cin44;
449 newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,
450 Cin41,Cin42,Cin43,Cin44); //get C matrix
452 Double_t Rold00=Cin00+V00;
453 Double_t Rold10=Cin10;
454 Double_t Rold11=Cin11+V11;
456 //////////////////////////////////// R matrix inversion ///////////////////////////////////////////////
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;
463 ////////////////////////////////////////////////////////////////////////////////////////////////////////
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;
476 Double_t X0,X1,X2,X3,X4;
477 newTrack->GetXElements(X0,X1,X2,X3,X4); // get the state vector
479 Double_t savex0=X0, savex1=X1;
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);
487 Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44;
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;
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;
500 C22=Cin22-K20*Cin20-K21*Cin21;
501 C32=Cin32-K20*Cin30-K21*Cin31;
502 C42=Cin42-K20*Cin40-K21*Cin41;
504 C33=Cin33-K30*Cin30-K31*Cin31;
505 C43=Cin43-K30*Cin40-K31*Cin41;
507 C44=Cin44-K40*Cin40-K41*Cin41;
509 newTrack->PutXElements(X0,X1,X2,X3,X4); // put the new state vector
511 newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the
513 Double_t VMCold00=V00-C00;
514 Double_t VMCold10=-C10;
515 Double_t VMCold11=V11-C11;
517 ///////////////////////////////////// Matrix VMC inversion ////////////////////////////////////////////////
519 det=VMCold00*VMCold11-VMCold10*VMCold10;
520 Double_t VMC00=VMCold11/det;
521 Double_t VMC10=-VMCold10/det;
522 Double_t VMC11=VMCold00/det;
524 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
526 Double_t chi2=(m[0]-X0)*( VMC00*(m[0]-X0) + 2.*VMC10*(m[1]-X1) ) +
527 (m[1]-X1)*VMC11*(m[1]-X1);
529 newTrack->SetChi2(newTrack->GetChi2()+chi2);
534 void 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
538 ////////////////////////////// Evaluation of the measurement vector m ///////////////
542 rk=cluster(0); phik=cluster(1); zk=cluster(2);
545 Double_t CC=(*newTrack).GetC();
546 Double_t Zv=(*newTrack).GetZv();
547 Double_t Dv=(*newTrack).GetDv();
549 Double_t tgl= (zk-Zv)*Cy/TMath::ASin(Cy*rk);
552 ///////////////////////////////////// Evaluation of the error matrix V //////////////
553 Int_t Layer=newTrack->GetLayer();
554 Double_t V00=sigma[0];
555 Double_t V11=sigma[1];
556 Double_t V31=sigma[1]/rk;
557 Double_t sigmaDv=newTrack->GetsigmaDv();
558 Double_t V22=sigmaDv*sigmaDv + newTrack->Getd2(Layer-1);
559 Double_t V32=newTrack->Getdtgl(Layer-1);
560 Double_t sigmaZv=newTrack->GetsigmaZv();
561 Double_t V33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
562 ///////////////////////////////////////////////////////////////////////////////////////
564 Double_t Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,Cin41,Cin42,Cin43,Cin44;
566 newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,
567 Cin41,Cin42,Cin43,Cin44); //get C matrix
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);
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];
584 ///////////////////// Matrix R inversion ////////////////////////////////////////////
597 for(j=k; j<n ; j++) {
598 for (i=j; i<n; i++) {
599 if(TMath::Abs(big) < TMath::Abs(R[i][j]) ) { big=R[i][j]; ll[k]=i; mm[k]=j; }
605 for(i=0; i<n; i++) { hold=-R[k][i]; R[k][i]=R[j][i]; R[j][i]=hold; }
611 for(j=0; j<n; j++) { hold=-R[j][k]; R[j][k]=R[j][i]; R[j][i]=hold; }
616 cout << "Singular matrix\n";
619 if(i == k) { continue; }
620 R[i][k]=R[i][k]/(-big);
626 if(i == k || j == k) { continue; }
627 R[i][j]=hold*R[k][j]+R[i][j];
632 if(j == k) { continue; }
641 for(k=n-1; k>=0; k--) {
644 for (j=0; j<n; j++) {hold=R[j][k]; R[j][k]=-R[j][i]; R[j][i]=hold;}
648 for (i=0; i<n; i++) {hold=R[k][i]; R[k][i]=-R[j][i]; R[j][i]=hold;}
651 //////////////////////////////////////////////////////////////////////////////////
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];
675 Double_t X0,X1,X2,X3,X4;
676 newTrack->GetXElements(X0,X1,X2,X3,X4); // get the state vector
678 Double_t savex0=X0, savex1=X1, savex2=X2, savex3=X3;
680 X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1)+K02*(m[2]-savex2)+
682 X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1)+K12*(m[2]-savex2)+
684 X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1)+K22*(m[2]-savex2)+
686 X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1)+K32*(m[2]-savex2)+
688 X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1)+K42*(m[2]-savex2)+
691 Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44;
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;
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;
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;
708 C33=Cin33-K30*Cin30-K31*Cin31-K32*Cin32-K33*Cin33;
709 C43=Cin43-K30*Cin40-K31*Cin41-K32*Cin42-K33*Cin43;
711 C44=Cin44-K40*Cin40-K41*Cin41-K42*Cin42-K43*Cin43;
713 newTrack->PutXElements(X0,X1,X2,X3,X4); // put the new state vector
715 newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the
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;
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];
730 /////////////////////// VMC matrix inversion ///////////////////////////////////
738 for(j=k; j<n ; j++) {
739 for (i=j; i<n; i++) {
740 if(TMath::Abs(big) < TMath::Abs(VMC[i][j]) ) { big=VMC[i][j]; ll[k]=i; mm[k]=j; }
746 for(i=0; i<n; i++) { hold=-VMC[k][i]; VMC[k][i]=VMC[j][i]; VMC[j][i]=hold; }
752 for(j=0; j<n; j++) { hold=-VMC[j][k]; VMC[j][k]=VMC[j][i]; VMC[j][i]=hold; }
757 cout << "Singular matrix\n";
760 if(i == k) { continue; }
761 VMC[i][k]=VMC[i][k]/(-big);
767 if(i == k || j == k) { continue; }
768 VMC[i][j]=hold*VMC[k][j]+VMC[i][j];
773 if(j == k) { continue; }
774 VMC[k][j]=VMC[k][j]/big;
782 for(k=n-1; k>=0; k--) {
785 for (j=0; j<n; j++) {hold=VMC[j][k]; VMC[j][k]=-VMC[j][i]; VMC[j][i]=hold;}
789 for (i=0; i<n; i++) {hold=VMC[k][i]; VMC[k][i]=-VMC[j][i]; VMC[j][i]=hold;}
794 ////////////////////////////////////////////////////////////////////////////////
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);
803 newTrack->SetChi2(newTrack->GetChi2()+chi2);