ITS/AliITStracking.cxx

   1 #include <iostream.h>
   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
  16 ClassImp(AliITStracking)
  17
  18
  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
  25
  26   Rlayer[0]=4.; Rlayer[1]=7.;  Rlayer[2]=14.9;  Rlayer[3]=23.8;  Rlayer[4]=39.1;  Rlayer[5]=43.6;
  27
  28   Int_t index;
  29   for(index =0; index<trackITSlist->GetSize(); index++) {
  30     AliITStrack *trackITS = (AliITStrack *) trackITSlist->At(index);
  31
  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";
  38    //  getchar();
  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     }
  52     Float_t NumClustNow = trackITS->GetNumClust();
  53     if(NumClustNow) {
  54       Chi2Now = trackITS->GetChi2();
  55       Chi2Now/=NumClustNow;
  56       //cout<<" Chi2Now =  "<<Chi2Now<<"\n";
  57     /*
  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;
  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;
  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};
  76     Int_t NDetector[]= {4,  4,   6,  8, 23, 26};
  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";
  83
  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
  92     Int_t flaghit=0;
  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;
 121       }
 122       Int_t iriv;
 123       for (iriv=0; iriv<idetot; iriv++) {  //for on detectors
 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();
 143         //gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty
 144         gAlice->TreeR()->GetEvent(index); //first entry in TreeR is empty
 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;
 159
 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];
 189
 190           sigma[0] = (Double_t)  recp->GetSigmaX2();
 191                          sigma[1] = (Double_t) recp->GetSigmaZ2();
 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
 196          // cout<<" layer = "<<play<<"\n";
 197          // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" "
 198          // <<vecclust(2)<<"\n"; getchar();
 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.;
 205             Float_t epsphi=5.0, epsz=5.0;
 206           if(Ptref<0.2) {epsphi=3.; epsz=3.;}
 207
 208           Double_t Rtrack=(*trackITS).Getrtrack();
 209           Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack);
 210           sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
 211
 212           sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
 213   //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
 214   //if(vecclust(3)==481) getchar();
 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;
 230
 231           (*newTrack).AddMS(rl);  // add the multiple scattering matrix to the covariance matrix
 232           (*newTrack).AddEL(rl,1.,0);
 233
 234           Double_t sigmanew[2];
 235           sigmanew[0]= sigmaphi;
 236           sigmanew[1]=sigma[1];
 237           //cout<<" Chiamo Kalman \n"; getchar();
 238
 239           if(flagvert)
 240             KalmanFilterVert(newTrack,cluster,sigmanew);
 241           else
 242             KalmanFilter(newTrack,cluster,sigmanew);
 243
 244
 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
 255
 256     }//end if(outinters==0)
 257
 258     if(flaghit==0 || outinters==-2) {
 259       AliITStrack *newTrack = new AliITStrack(*trackITS);
 260       (*newTrack).SetLayer((*trackITS).GetLayer()-1);
 261       (*newTrack).AddMS(rl);  // add the multiple scattering matrix to the covariance matrix
 262       (*newTrack).AddEL(rl,1.,0);
 263
 264       listoftrack.AddLast(newTrack);
 265     }
 266
 267
 268     //gObjectTable->Print();   // stampa memoria
 269
 270     AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert,rl);
 271     listoftrack.Delete();
 272   } // end of for on tracks
 273
 274   //gObjectTable->Print();   // stampa memoria
 275
 276 }
 277
 278
 279 Int_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};
 294   Int_t NDetector[]= {4,  4,   6,  8, 23, 26};
 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
 308   Int_t iD;
 309   for(iD = 1; iD<= NDetector[layer-1]; iD++) {
 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;
 315     //added to take into account problem on drift
 316     if(layer == 4 || layer==3) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
 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();
 337
 338   Int_t iLd;
 339   for(iLd = 1; iLd<= NLadder[layer-1]; iLd++) {
 340           g1->GetCenterThetaPhi(layer,iLd,detector,det);
 341   Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));
 342   // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
 343   Double_t xmin,ymin,xmax,ymax;
 344  // Double_t phiconfr=0.0;
 345   //cout<<" phiconfr inizio =  "<<phiconfr <<"\n"; getchar();
 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";
 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();}
 364     phimin=phimin-(2.*pigre);
 365     if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
 366     if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
 367   }
 368   //  cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar();
 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
 386 Double_t AliITStracking::PhiDef(Double_t x, Double_t y){
 387   Double_t pigre= TMath::Pi();
 388   Double_t phi=0.0;
 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
 411 void 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
 473 void 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();
 515
 516
 517
 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;
 523   TVector  x=newTrack->GetVector();
 524   TVector savex=x;
 525   x*=H; x-=m;
 526   x*=-1; x*=K; x+=savex;
 527   TMatrix saveC=C;
 528   C.Mult(K,tmp); C-=saveC; C*=-1;
 529
 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