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 <TObjectTable.h>
   7
   8 #include "AliITStracking.h"
   9 #include "AliRun.h"
  10 #include "AliITS.h"
  11 #include "AliITSgeom.h"
  12 #include "AliITSRecPoint.h"
  13 #include "AliITSRad.h"
  14 #include "AliITStrack.h"
  15 #include "AliITSgeoinfo.h"
  16
  17 ClassImp(AliITStracking)
  18
  19
  20 AliITStracking::AliITStracking(TList *trackITSlist, AliITStrack *reference,
  21                 AliITS *aliITS, TObjArray *rpoints, Double_t Ptref, Int_t **vettid, Bool_t flagvert,
  22                                          AliITSRad *rl, AliITSgeoinfo *geoinfo) {
  23 ///////////////////////   This function perform the tracking in ITS detectors /////////////////////
  24 ///////////////////////     reference is a pointer to the final best track    /////////////////////
  25 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
  26 // The authors  thank   Mariana Bondila to have help them to resolve some problems.  July-2000
  27
  28   //Rlayer[0]=4.; Rlayer[1]=7.;  Rlayer[2]=14.9;  Rlayer[3]=23.8;  Rlayer[4]=39.1;  Rlayer[5]=43.6; //vecchio
  29
  30   Int_t index;
  31   for(index =0; index<trackITSlist->GetSize(); index++) {
  32     AliITStrack *trackITS = (AliITStrack *) trackITSlist->At(index);
  33
  34     if((*trackITS).GetLayer()==7) reference->SetChi2(10.223e140);
  35    // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n";
  36    //  cout<<"fvtrack =" <<"\n";
  37    //  cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n";
  38    //  cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n";
  39    //  cout<< " Pt = "<<(*trackITS).GetPt()<<"\n";
  40    //  getchar();
  41     Double_t Chi2Now, Chi2Ref;
  42     if((*trackITS).GetLayer()==1 ) {
  43       Chi2Now = trackITS->GetChi2();
  44       Float_t NumClustNow = trackITS->GetNumClust();
  45       if(trackITS->GetNumClust()) Chi2Now /= (Double_t )trackITS->GetNumClust();
  46       Chi2Ref = reference->GetChi2();
  47       Float_t NumClustRef = reference->GetNumClust();
  48       if(reference->GetNumClust()) Chi2Ref /= (Double_t )reference->GetNumClust();
  49       //cout<<" Chi2Now and Chi2Ref = "<<Chi2Now<<" "<<Chi2Ref<<"\n";
  50                 if( NumClustNow > NumClustRef ) {*reference = *trackITS;}
  51       if((NumClustNow == NumClustRef )&& (Chi2Now < Chi2Ref))  {*reference = *trackITS;}
  52       continue;
  53     }
  54     Float_t NumClustNow = trackITS->GetNumClust();
  55     if(NumClustNow) {
  56       Chi2Now = trackITS->GetChi2();
  57       Chi2Now/=NumClustNow;
  58       //cout<<" Chi2Now =  "<<Chi2Now<<"\n";
  59     /*
  60     // if(Ptref > 0.6 && Chi2Now > 20.) continue;
  61     if(Ptref > 0.6 && Chi2Now > 30.) continue;
  62     if((Ptref <= 0.6 && Ptref>0.2)&& Chi2Now > 15.) continue;
  63      // if(Chi2Now>5.) continue;
  64       //if(Chi2Now>15.) continue;
  65      // if(Ptref <= 0.2 && Chi2Now > 10.) continue;
  66      if(Ptref <= 0.2 && Chi2Now > 8.) continue;
  67      */
  68      if(Ptref > 1.0 && Chi2Now > 30.) continue;
  69      if((Ptref >= 0.6 && Ptref<=1.0) && Chi2Now > 40.) continue;
  70      if((Ptref <= 0.6 && Ptref>0.2)&& Chi2Now > 40.) continue;
  71      if(Ptref <= 0.2 && Chi2Now > 8.) continue;
  72     }
  73
  74     Int_t layerInit = (*trackITS).GetLayer();
  75     Int_t layernew = layerInit - 2;  // -1 for new layer, -1 for matrix index
  76
  77     //Int_t NLadder[]= {20, 40, 14, 22, 34, 38};   //vecchio
  78     //Int_t NDetector[]= {4,  4,   6,  8, 23, 26}; //vecchio
  79
  80     TList listoftrack;
  81     Int_t ladp, ladm, detp,detm,ladinters,detinters;
  82     Int_t layerfin=layerInit-1;
  83     //Double_t Rfin=Rlayer[layerfin-1];  // vecchio
  84          Double_t Rfin=geoinfo->Avrad[layerfin-1];  // nuovo
  85     // cout<<"Prima di intersection \n";
  86
  87     Int_t  outinters=NewIntersection(*trackITS, Rfin, layerfin, ladinters, detinters, geoinfo);
  88
  89    // cout<<" outinters = "<<outinters<<"\n";
  90    //  cout<<" Layer ladder detector intersection ="<<layerfin<<" "<<ladinters<<" "<<detinters<<"\n";
  91    //  cout << " phiinters zinters = "<<(*trackITS)(0) << " "<<(*trackITS)(1)<<"\n"; getchar();
  92
  93     if(outinters==-1) continue;
  94
  95     Int_t flaghit=0;
  96     if(outinters==0){
  97       TVector Touclad(9), Toucdet(9);
  98       Int_t lycur=layerfin;
  99       ladp=ladinters+1;
 100       ladm=ladinters-1;
 101       //if(ladm <= 0) ladm=NLadder[layerfin-1];  //vecchio
 102                 if(ladm <= 0) ladm=geoinfo->Nlad[layerfin-1];    //nuovo
 103       //if(ladp > NLadder[layerfin-1]) ladp=1;          //vecchio
 104                 if(ladp > geoinfo->Nlad[layerfin-1]) ladp=1;  //nuovo
 105       detp=detinters+1;
 106       detm=detinters-1;
 107       Int_t idetot=1;
 108       Touclad(0)=ladinters; Touclad(1)=ladm; Touclad(2)=ladp;
 109       Touclad(3)=ladinters; Touclad(4)=ladm; Touclad(5)=ladp;
 110       Touclad(6)=ladinters; Touclad(7)=ladm; Touclad(8)=ladp;
 111       Toucdet(0)=detinters; Toucdet(1)=detinters; Toucdet(2)=detinters;
 112       //if(detm > 0 && detp <= NDetector[layerfin-1]) {   //vecchio
 113                 if(detm > 0 && detp <= geoinfo->Ndet[layerfin-1]) {     //nuovo
 114         idetot=9;
 115         Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
 116         Toucdet(6)=detp; Toucdet(7)=detp; Toucdet(8)=detp;
 117       }
 118
 119       //if(detm > 0 && detp > NDetector[layerfin-1]) {  //vecchio
 120                 if(detm > 0 && detp > geoinfo->Ndet[layerfin-1]) {   //nuovo
 121         idetot=6;
 122         Toucdet(3)=detm; Toucdet(4)=detm; Toucdet(5)=detm;
 123       }
 124
 125       //if(detm <= 0 && detp <= NDetector[layerfin-1]) {  //vecchio
 126                 if(detm <= 0 && detp <= geoinfo->Ndet[layerfin-1]) {   //nuovo
 127         idetot=6;
 128         Toucdet(3)=detp; Toucdet(4)=detp; Toucdet(5)=detp;
 129       }
 130       Int_t iriv;
 131       for (iriv=0; iriv<idetot; iriv++) {  //for on detectors
 132         AliITSgeom *g1 = aliITS->GetITSgeom();
 133         TVector CTF(9);
 134         g1->GetCenterThetaPhi(layerInit-1,(Int_t)Touclad(iriv),(Int_t)Toucdet(iriv),CTF);
 135
 136         // cout<<" layer, ladder, det, xo, yo, zo = "<<layerInit-1<<" "<<(Int_t)Touclad(iriv)<<
 137         // " "<<(Int_t)Toucdet(iriv)<<" "<<CTF(0)<<" "<<CTF(1)<<" "<<CTF(2)<< " "<<CTF(6)<<"\n"; getchar();
 138
 139         ////////////////////////////////////////////////////////////////////////////////////////////////
 140
 141         /*** Rec points sorted by module *****/
 142         /**************************************/
 143
 144         Int_t index;
 145         AliITSRecPoint *recp;
 146         AliITSgeom *geom = aliITS->GetITSgeom();
 147         index = geom->GetModuleIndex(lycur,Touclad(iriv),Toucdet(iriv));
 148         Int_t lay,lad,det;
 149         geom->GetModuleId(index,lay,lad,det);
 150         aliITS->ResetRecPoints();
 151         //gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty
 152         gAlice->TreeR()->GetEvent(index); //first entry in TreeR is empty
 153
 154         Int_t npoints=rpoints->GetEntries();
 155         Int_t *indlist=new Int_t[npoints+1];
 156         Int_t counter=0;
 157         Int_t ind;
 158         for (ind=0; ind<=npoints; ind++) {
 159           indlist[ind]=-1;
 160                if (*(vettid[index]+ind)==0) {
 161               indlist[counter]=ind;
 162                    counter++;
 163               }
 164         }
 165
 166         ind=-1;
 167
 168         for(;;) {
 169           ind++;
 170           if(indlist[ind] < 0) recp=0;
 171                else recp = (AliITSRecPoint*)rpoints->UncheckedAt(indlist[ind]);
 172
 173           if((!recp)  )  break;
 174           TVector cluster(3),vecclust(9);
 175           vecclust(6)=vecclust(7)=vecclust(8)=-1.;
 176           Double_t sigma[2];
 177             //modificata 8-3-2001
 178           // set veclust in global
 179           Float_t global[3], local[3];
 180           local[0]=recp->GetX();
 181           local[1]=0.;
 182           local[2]= recp->GetZ();
 183           AliITSgeom *g1 = aliITS->GetITSgeom();
 184           Int_t play = lycur;
 185           Int_t plad = TMath::Nint(Touclad(iriv));
 186           Int_t pdet = TMath::Nint(Toucdet(iriv));
 187           g1->LtoG(play,plad,pdet,local,global);
 188
 189           vecclust(0)=global[0];
 190           vecclust(1)=global[1];
 191           vecclust(2)=global[2];
 192
 193           /*
 194           ////  modificato 8-3-2001
 195           vecclust(0)=recp->GetRhit();;
 196           vecclust(1)=recp->Getphi();
 197           vecclust(2)=recp->GetZglobal();
 198           */
 199
 200           vecclust(3) = (float)recp->fTracks[0];
 201           vecclust(4) = (float)indlist[ind];
 202           vecclust(5) = (float)index;
 203           vecclust(6) = (float)recp->fTracks[0];
 204           vecclust(7) = (float)recp->fTracks[1];
 205           vecclust(8) = (float)recp->fTracks[2];
 206
 207           sigma[0] = (Double_t)  recp->GetSigmaX2();
 208           sigma[1] = (Double_t) recp->GetSigmaZ2();
 209            //commentato 8-3-2001
 210          //now we are in r,phi,z in global
 211           cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit
 212           cluster(1) = PhiDef(vecclust(0),vecclust(1));    // phi hit
 213           cluster(2) = vecclust(2);                   // z hit
 214
 215           /*
 216          //modificato 8-3-2001
 217           cluster(0) = vecclust(0);//r hit
 218           cluster(1) = vecclust(1);    // phi hit
 219           cluster(2) = vecclust(2);                   // z hit
 220           */
 221          // cout<<" layer = "<<play<<"\n";
 222          // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" "
 223          // <<vecclust(2)<<"\n"; getchar();
 224           //cluster(1)= cluster(1)-trackITS->Getalphaprov();  //provvisorio;
 225                          //if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi(); //provvisorio
 226                          //cout<<" cluster(1) dopo = "<<cluster(1)<< " alphaprov = "<<trackITS->Getalphaprov()<<"\n";
 227           Float_t sigmatotphi, sigmatotz;
 228
 229           //Float_t epsphi=3.2, epsz=3.;
 230             Float_t epsphi=5.0, epsz=5.0;
 231           if(Ptref<0.2) {epsphi=3.; epsz=3.;}
 232
 233           Double_t Rtrack=(*trackITS).Getrtrack();
 234           Double_t sigmaphi=sigma[0]/(Rtrack*Rtrack);
 235           sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
 236
 237           sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
 238   //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
 239   //if(vecclust(3)==481) getchar();
 240                if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi());
 241                if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi());
 242           if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue;
 243                         // cout<<" supero sigmaphi \n";
 244           AliITStrack *newTrack = new AliITStrack((*trackITS));
 245           (*newTrack).SetLayer((*trackITS).GetLayer()-1);
 246
 247                if (TMath::Abs(Rtrack-cluster(0))/Rtrack>1e-6)
 248                                                (*newTrack).Correct(Double_t(cluster(0)));
 249                         //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<<        (*newTrack).GetZ()<<"\n";
 250           if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){
 251              delete newTrack;
 252              continue;}
 253
 254           if(iriv == 0) flaghit=1;
 255
 256           (*newTrack).AddMS(rl);  // add the multiple scattering matrix to the covariance matrix
 257           (*newTrack).AddEL(rl,1.,0);
 258
 259           Double_t sigmanew[2];
 260           sigmanew[0]= sigmaphi;
 261           sigmanew[1]=sigma[1];
 262
 263           if(flagvert)
 264             KalmanFilterVert(newTrack,cluster,sigmanew);
 265           else
 266             KalmanFilter(newTrack,cluster,sigmanew);
 267
 268
 269           (*newTrack).PutCluster(layernew, vecclust);
 270            newTrack->AddClustInTrack();
 271
 272            listoftrack.AddLast(newTrack);
 273
 274         }   // end of for(;;) on rec points
 275
 276         delete [] indlist;
 277
 278       }  // end of for on detectors
 279
 280     }//end if(outinters==0)
 281
 282     if(flaghit==0 || outinters==-2) {
 283       AliITStrack *newTrack = new AliITStrack(*trackITS);
 284       (*newTrack).SetLayer((*trackITS).GetLayer()-1);
 285       (*newTrack).AddMS(rl);  // add the multiple scattering matrix to the covariance matrix
 286       (*newTrack).AddEL(rl,1.,0);
 287
 288       listoftrack.AddLast(newTrack);
 289     }
 290
 291
 292     //gObjectTable->Print();   // stampa memoria
 293
 294     AliITStracking(&listoftrack, reference, aliITS, rpoints,Ptref,vettid,flagvert,rl, geoinfo);
 295     listoftrack.Delete();
 296   } // end of for on tracks
 297
 298   //gObjectTable->Print();   // stampa memoria
 299
 300 }
 301
 302 Int_t AliITStracking::NewIntersection(AliITStrack &track, Double_t rk,Int_t layer, Int_t &ladder,
 303 Int_t &detector, AliITSgeoinfo *geoinfo) {
 304 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 305 // Found the intersection and the detector
 306
 307   if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;}
 308   track.Propagation(rk);
 309   Double_t zinters=track.GetZ();
 310   Double_t phinters=track.Getphi();
 311   //phinters = phinters+track.Getalphaprov(); //provvisorio
 312   //if(phinters>2.*3.14) phinters=phinters-2.*3.14; //provvisorio
 313   //cout<<"zinters = "<<zinters<<"  phinters = "<<phinters<<"\n";
 314
 315   //////////////////////////////////      limits for Geometry 5      /////////////////////////////
 316
 317   //Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
 318   //Int_t NDetector[]= {4,  4,   6,  8, 23, 26};
 319
 320   //Float_t Detx[]= {0.64, 0.64, 3.509, 3.509, 3.65, 3.65 };
 321   //Float_t Detz[]= {4.19, 4.19, 3.75 , 3.75 , 2   , 2    };
 322
 323   ////////////////////////////////////////////////////////////////////////////////////////////////
 324
 325   TVector det(9);
 326   TVector ListDet(2);
 327   TVector DistZCenter(2);
 328   AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
 329
 330   Int_t iz=0;
 331   Double_t epsz=1.2;
 332   Double_t epszpixel=0.05;
 333
 334   Int_t iD;
 335   //for(iD = 1; iD<= NDetector[layer-1]; iD++) {   //vecchio
 336   for(iD = 1; iD<= geoinfo->Ndet[layer-1]; iD++) {
 337     g1->GetCenterThetaPhi(layer,1,iD,det);
 338     //Double_t zmin=det(2)-Detz[layer-1];   //vecchio
 339          Double_t zmin=det(2)-geoinfo->Detz[layer-1];   //nuovo
 340     //if(iD==1) zmin=det(2)-(Detz[layer-1])*epsz;               //vecchio
 341          if(iD==1) zmin=det(2)-(geoinfo->Detz[layer-1])*epsz;           //nuovo
 342     //Double_t zmax=det(2)+Detz[layer-1];    //vecchio
 343          Double_t zmax=det(2)+geoinfo->Detz[layer-1];    //nuovo
 344     //if(iD==NDetector[layer-1]) zmax=det(2)+(Detz[layer-1])*epsz;  //vecchio
 345          if(iD==geoinfo->Ndet[layer-1]) zmax=det(2)+(geoinfo->Detz[layer-1])*epsz;   //nuovo
 346     //added to take into account problem on drift
 347     if(layer == 4 || layer==3) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
 348     //cout<<"zmin zinters zmax det(2)= "<<zmin<<" "<<zinters<<" "<<zmax<<" "<<det(2)<<"\n";
 349     if(zinters > zmin && zinters <= zmax) {
 350       if(iz>1) {cout<< " Errore su iz in NewIntersection \n"; getchar();}
 351       else {
 352         ListDet(iz)= iD; DistZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
 353       }
 354     }
 355   }
 356
 357   if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;}
 358   detector=Int_t (ListDet(0));
 359   if(iz>1 && (DistZCenter(0)>DistZCenter(1)))   detector=Int_t (ListDet(1));
 360
 361   AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
 362   Float_t global[3];
 363   Float_t local[3];
 364   TVector ListLad(2);
 365   TVector DistphiCenter(2);
 366   Int_t ip=0;
 367   Double_t pigre=TMath::Pi();
 368
 369   Int_t iLd;
 370   //for(iLd = 1; iLd<= NLadder[layer-1]; iLd++) {   //vecchio
 371   for(iLd = 1; iLd<= geoinfo->Nlad[layer-1]; iLd++) {  //nuovo
 372           g1->GetCenterThetaPhi(layer,iLd,detector,det);
 373   Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));
 374   // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
 375   Double_t xmin,ymin,xmax,ymax;
 376  // Double_t phiconfr=0.0;
 377   //cout<<" phiconfr inizio =  "<<phiconfr <<"\n"; getchar();
 378   local[1]=local[2]=0.;
 379   //local[0]= -(Detx[layer-1]);  //vecchio
 380   local[0]= -(geoinfo->Detx[layer-1]);    //nuovo
 381   //if(layer==1)    local[0]= (Detx[layer-1]);  //take into account different reference system  //vecchio
 382   if(layer==1)    local[0]= (geoinfo->Detx[layer-1]);  //take into account different reference system   //nuovo
 383   g2->LtoG(layer,iLd,detector,local,global);
 384   xmax=global[0]; ymax=global[1];
 385   //local[0]= (Detx[layer-1]);   //vecchio
 386   local[0]= (geoinfo->Detx[layer-1]);   //nuovo
 387   //if(layer==1)    local[0]= -(Detx[layer-1]);  //take into account different reference system //vecchio
 388   if(layer==1)    local[0]= -(geoinfo->Detx[layer-1]);  //take into account different reference system //nuovo
 389   g2->LtoG(layer,iLd,detector,local,global);
 390   xmin=global[0]; ymin=global[1];
 391   Double_t phimin=PhiDef(xmin,ymin);
 392   Double_t phimax=PhiDef(xmax,ymax);
 393   //cout<<" xmin ymin = "<<xmin<<" "<<ymin<<"\n";
 394   // cout<<" xmax ymax = "<<xmax<<" "<<ymax<<"\n";
 395   // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n";
 396
 397   Double_t phiconfr=phinters;
 398  if(phimin>phimax ){
 399      if(phimin <5.5) {cout<<" Error in NewIntersection for phi \n"; getchar();}
 400     phimin=phimin-(2.*pigre);
 401     if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
 402     if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
 403   }
 404   //  cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar();
 405   if(phiconfr>phimin && phiconfr<= phimax) {
 406     if(ip>1) {
 407       cout<< " Errore su ip in NewIntersection \n"; getchar();
 408     }
 409       else  {
 410         ListLad(ip)= iLd; DistphiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
 411       }
 412     }
 413   }
 414   if(ip==0) { cout<< " No detector along phi \n"; getchar();}
 415   ladder=Int_t (ListLad(0));
 416   if(ip>1 && (DistphiCenter(0)>DistphiCenter(1)))   ladder=Int_t (ListLad(1));
 417
 418   return 0;
 419 }
 420
 421
 422 Double_t AliITStracking::PhiDef(Double_t x, Double_t y){
 423   Double_t pigre= TMath::Pi();
 424   Double_t phi=0.0;
 425   if(y == 0. || x == 0.) {
 426     if(y == 0. && x == 0.) {
 427       cout << "  Error in AliITStracking::PhiDef x=0 and y=0 \n"; getchar();
 428     }
 429     if(y==0. && x>0.) phi=0.;
 430     if(y==0. && x<0.) phi=pigre;
 431     if(x==0 && y>0.) phi=pigre/2.;
 432     if(x==0 && y<0.) phi=1.5*pigre;
 433   }
 434     else {
 435       if (x>0. && y>0.) phi=TMath::ATan(y/x);
 436       if (x<0. && y>0.) phi=pigre+TMath::ATan(y/x);
 437       if (x<0. && y<0.) phi=pigre+TMath::ATan(y/x);
 438       if (x>0. && y<0.) phi=(2.*pigre)+TMath::ATan(y/x);
 439     }
 440   if(phi<0. || phi>(2*pigre)) {
 441     cout<<" Error on phi in  AliITStracking::PhiDef \n"; getchar();
 442   }
 443   return phi;
 444 }
 445
 446
 447 void AliITStracking::KalmanFilter(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
 448 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 449 // Kalman filter without vertex constraint
 450
 451
 452   ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
 453
 454   Double_t m[2];
 455   Double_t rk,phik,zk;
 456   rk=cluster(0);   phik=cluster(1);  zk=cluster(2);
 457   m[0]=phik;    m[1]=zk;
 458
 459   ///////////////////////////////////// Evaluation of the error matrix V  ///////////////////////////////
 460
 461   Double_t V00=sigma[0];
 462   Double_t V11=sigma[1];
 463
 464   ///////////////////////////////////////////////////////////////////////////////////////////
 465
 466
 467   Double_t Cin00,Cin10,Cin20,Cin30,Cin40,Cin11,Cin21,Cin31,Cin41,Cin22,Cin32,Cin42,Cin33,Cin43,Cin44;
 468
 469   newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,
 470                          Cin41,Cin42,Cin43,Cin44); //get C matrix
 471
 472   Double_t Rold00=Cin00+V00;
 473   Double_t Rold10=Cin10;
 474   Double_t Rold11=Cin11+V11;
 475
 476 //////////////////////////////////// R matrix inversion  ///////////////////////////////////////////////
 477
 478   Double_t det=Rold00*Rold11-Rold10*Rold10;
 479   Double_t R00=Rold11/det;
 480   Double_t R10=-Rold10/det;
 481   Double_t R11=Rold00/det;
 482
 483 ////////////////////////////////////////////////////////////////////////////////////////////////////////
 484
 485   Double_t K00=Cin00*R00+Cin10*R10;
 486   Double_t K01=Cin00*R10+Cin10*R11;
 487   Double_t K10=Cin10*R00+Cin11*R10;
 488   Double_t K11=Cin10*R10+Cin11*R11;
 489   Double_t K20=Cin20*R00+Cin21*R10;
 490   Double_t K21=Cin20*R10+Cin21*R11;
 491   Double_t K30=Cin30*R00+Cin31*R10;
 492   Double_t K31=Cin30*R10+Cin31*R11;
 493   Double_t K40=Cin40*R00+Cin41*R10;
 494   Double_t K41=Cin40*R10+Cin41*R11;
 495
 496   Double_t X0,X1,X2,X3,X4;
 497   newTrack->GetXElements(X0,X1,X2,X3,X4);     // get the state vector
 498
 499   Double_t savex0=X0, savex1=X1;
 500
 501   X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1);
 502   X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1);
 503   X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1);
 504   X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1);
 505   X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1);
 506
 507   Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44;
 508
 509   C00=Cin00-K00*Cin00-K01*Cin10;
 510   C10=Cin10-K00*Cin10-K01*Cin11;
 511   C20=Cin20-K00*Cin20-K01*Cin21;
 512   C30=Cin30-K00*Cin30-K01*Cin31;
 513   C40=Cin40-K00*Cin40-K01*Cin41;
 514
 515   C11=Cin11-K10*Cin10-K11*Cin11;
 516   C21=Cin21-K10*Cin20-K11*Cin21;
 517   C31=Cin31-K10*Cin30-K11*Cin31;
 518   C41=Cin41-K10*Cin40-K11*Cin41;
 519
 520   C22=Cin22-K20*Cin20-K21*Cin21;
 521   C32=Cin32-K20*Cin30-K21*Cin31;
 522   C42=Cin42-K20*Cin40-K21*Cin41;
 523
 524   C33=Cin33-K30*Cin30-K31*Cin31;
 525   C43=Cin43-K30*Cin40-K31*Cin41;
 526
 527   C44=Cin44-K40*Cin40-K41*Cin41;
 528
 529   newTrack->PutXElements(X0,X1,X2,X3,X4);               // put the new state vector
 530
 531   newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the
 532                                                                                        // new cov matrix
 533   Double_t VMCold00=V00-C00;
 534   Double_t VMCold10=-C10;
 535   Double_t VMCold11=V11-C11;
 536
 537 ///////////////////////////////////// Matrix VMC inversion  ////////////////////////////////////////////////
 538
 539   det=VMCold00*VMCold11-VMCold10*VMCold10;
 540   Double_t VMC00=VMCold11/det;
 541   Double_t VMC10=-VMCold10/det;
 542   Double_t VMC11=VMCold00/det;
 543
 544 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
 545
 546   Double_t chi2=(m[0]-X0)*( VMC00*(m[0]-X0) + 2.*VMC10*(m[1]-X1) ) +
 547                 (m[1]-X1)*VMC11*(m[1]-X1);
 548
 549   newTrack->SetChi2(newTrack->GetChi2()+chi2);
 550
 551 }
 552
 553
 554 void AliITStracking::KalmanFilterVert(AliITStrack *newTrack,TVector &cluster,Double_t sigma[2]){
 555 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 556 // Kalman filter with vertex constraint
 557
 558   ////////////////////////////// Evaluation of the measurement vector m ///////////////
 559
 560   Double_t m[4];
 561   Double_t rk,phik,zk;
 562   rk=cluster(0);   phik=cluster(1);  zk=cluster(2);
 563   m[0]=phik;    m[1]=zk;
 564
 565   Double_t CC=(*newTrack).GetC();
 566   Double_t Zv=(*newTrack).GetZv();
 567   Double_t Dv=(*newTrack).GetDv();
 568   Double_t Cy=CC/2.;
 569   Double_t tgl= (zk-Zv)*Cy/TMath::ASin(Cy*rk);
 570   m[2]=Dv;    m[3]=tgl;
 571
 572   ///////////////////////////////////// Evaluation of the error matrix V  //////////////
 573   Int_t Layer=newTrack->GetLayer();
 574   Double_t V00=sigma[0];
 575   Double_t V11=sigma[1];
 576   Double_t V31=sigma[1]/rk;
 577   Double_t sigmaDv=newTrack->GetsigmaDv();
 578   Double_t V22=sigmaDv*sigmaDv  + newTrack->Getd2(Layer-1);
 579   Double_t V32=newTrack->Getdtgl(Layer-1);
 580   Double_t sigmaZv=newTrack->GetsigmaZv();
 581   Double_t V33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
 582   ///////////////////////////////////////////////////////////////////////////////////////
 583
 584   Double_t Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,Cin41,Cin42,Cin43,Cin44;
 585
 586   newTrack->GetCElements(Cin00,Cin10,Cin11,Cin20,Cin21,Cin22,Cin30,Cin31,Cin32,Cin33,Cin40,
 587                          Cin41,Cin42,Cin43,Cin44); //get C matrix
 588
 589   Double_t R[4][4];
 590   R[0][0]=Cin00+V00;
 591   R[1][0]=Cin10;
 592   R[2][0]=Cin20;
 593   R[3][0]=Cin30;
 594   R[1][1]=Cin11+V11;
 595   R[2][1]=Cin21;
 596   R[3][1]=Cin31+sigma[1]/rk;
 597   R[2][2]=Cin22+sigmaDv*sigmaDv+newTrack->Getd2(Layer-1);
 598   R[3][2]=Cin32+newTrack->Getdtgl(Layer-1);
 599   R[3][3]=Cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(Layer-1);
 600
 601   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];
 602   R[2][3]=R[3][2];
 603
 604 /////////////////////  Matrix R inversion ////////////////////////////////////////////
 605
 606   const Int_t n=4;
 607   Double_t big, hold;
 608   Double_t d=1.;
 609   Int_t ll[n],mm[n];
 610
 611   Int_t i,j,k;
 612
 613   for(k=0; k<n; k++) {
 614     ll[k]=k;
 615     mm[k]=k;
 616     big=R[k][k];
 617     for(j=k; j<n ; j++) {
 618       for (i=j; i<n; i++) {
 619         if(TMath::Abs(big) < TMath::Abs(R[i][j]) ) { big=R[i][j]; ll[k]=i; mm[k]=j; }
 620       }
 621     }
 622 //
 623     j= ll[k];
 624     if(j > k) {
 625       for(i=0; i<n; i++) { hold=-R[k][i]; R[k][i]=R[j][i]; R[j][i]=hold; }
 626
 627     }
 628 //
 629     i=mm[k];
 630     if(i > k ) {
 631       for(j=0; j<n; j++) { hold=-R[j][k]; R[j][k]=R[j][i]; R[j][i]=hold; }
 632     }
 633 //
 634     if(!big) {
 635       d=0.;
 636       cout << "Singular matrix\n";
 637     }
 638     for(i=0; i<n; i++) {
 639       if(i == k) { continue; }
 640       R[i][k]=R[i][k]/(-big);
 641     }
 642 //
 643     for(i=0; i<n; i++) {
 644       hold=R[i][k];
 645       for(j=0; j<n; j++) {
 646         if(i == k || j == k) { continue; }
 647         R[i][j]=hold*R[k][j]+R[i][j];
 648       }
 649     }
 650 //
 651     for(j=0; j<n; j++) {
 652       if(j == k) { continue; }
 653       R[k][j]=R[k][j]/big;
 654     }
 655 //
 656     d=d*big;
 657 //
 658     R[k][k]=1./big;
 659   }
 660 //
 661   for(k=n-1; k>=0; k--) {
 662     i=ll[k];
 663     if(i > k) {
 664       for (j=0; j<n; j++) {hold=R[j][k]; R[j][k]=-R[j][i]; R[j][i]=hold;}
 665     }
 666     j=mm[k];
 667     if(j > k) {
 668       for (i=0; i<n; i++) {hold=R[k][i]; R[k][i]=-R[j][i]; R[j][i]=hold;}
 669       }
 670   }
 671 //////////////////////////////////////////////////////////////////////////////////
 672
 673
 674   Double_t K00=Cin00*R[0][0]+Cin10*R[1][0]+Cin20*R[2][0]+Cin30*R[3][0];
 675   Double_t K01=Cin00*R[1][0]+Cin10*R[1][1]+Cin20*R[2][1]+Cin30*R[3][1];
 676   Double_t K02=Cin00*R[2][0]+Cin10*R[2][1]+Cin20*R[2][2]+Cin30*R[3][2];
 677   Double_t K03=Cin00*R[3][0]+Cin10*R[3][1]+Cin20*R[3][2]+Cin30*R[3][3];
 678   Double_t K10=Cin10*R[0][0]+Cin11*R[1][0]+Cin21*R[2][0]+Cin31*R[3][0];
 679   Double_t K11=Cin10*R[1][0]+Cin11*R[1][1]+Cin21*R[2][1]+Cin31*R[3][1];
 680   Double_t K12=Cin10*R[2][0]+Cin11*R[2][1]+Cin21*R[2][2]+Cin31*R[3][2];
 681   Double_t K13=Cin10*R[3][0]+Cin11*R[3][1]+Cin21*R[3][2]+Cin31*R[3][3];
 682   Double_t K20=Cin20*R[0][0]+Cin21*R[1][0]+Cin22*R[2][0]+Cin32*R[3][0];
 683   Double_t K21=Cin20*R[1][0]+Cin21*R[1][1]+Cin22*R[2][1]+Cin32*R[3][1];
 684   Double_t K22=Cin20*R[2][0]+Cin21*R[2][1]+Cin22*R[2][2]+Cin32*R[3][2];
 685   Double_t K23=Cin20*R[3][0]+Cin21*R[3][1]+Cin22*R[3][2]+Cin32*R[3][3];
 686   Double_t K30=Cin30*R[0][0]+Cin31*R[1][0]+Cin32*R[2][0]+Cin33*R[3][0];
 687   Double_t K31=Cin30*R[1][0]+Cin31*R[1][1]+Cin32*R[2][1]+Cin33*R[3][1];
 688   Double_t K32=Cin30*R[2][0]+Cin31*R[2][1]+Cin32*R[2][2]+Cin33*R[3][2];
 689   Double_t K33=Cin30*R[3][0]+Cin31*R[3][1]+Cin32*R[3][2]+Cin33*R[3][3];
 690   Double_t K40=Cin40*R[0][0]+Cin41*R[1][0]+Cin42*R[2][0]+Cin43*R[3][0];
 691   Double_t K41=Cin40*R[1][0]+Cin41*R[1][1]+Cin42*R[2][1]+Cin43*R[3][1];
 692   Double_t K42=Cin40*R[2][0]+Cin41*R[2][1]+Cin42*R[2][2]+Cin43*R[3][2];
 693   Double_t K43=Cin40*R[3][0]+Cin41*R[3][1]+Cin42*R[3][2]+Cin43*R[3][3];
 694
 695   Double_t X0,X1,X2,X3,X4;
 696   newTrack->GetXElements(X0,X1,X2,X3,X4);     // get the state vector
 697
 698   Double_t savex0=X0, savex1=X1, savex2=X2, savex3=X3;
 699
 700   X0+=K00*(m[0]-savex0)+K01*(m[1]-savex1)+K02*(m[2]-savex2)+
 701       K03*(m[3]-savex3);
 702   X1+=K10*(m[0]-savex0)+K11*(m[1]-savex1)+K12*(m[2]-savex2)+
 703       K13*(m[3]-savex3);
 704   X2+=K20*(m[0]-savex0)+K21*(m[1]-savex1)+K22*(m[2]-savex2)+
 705       K23*(m[3]-savex3);
 706   X3+=K30*(m[0]-savex0)+K31*(m[1]-savex1)+K32*(m[2]-savex2)+
 707       K33*(m[3]-savex3);
 708   X4+=K40*(m[0]-savex0)+K41*(m[1]-savex1)+K42*(m[2]-savex2)+
 709       K43*(m[3]-savex3);
 710
 711   Double_t C00,C10,C20,C30,C40,C11,C21,C31,C41,C22,C32,C42,C33,C43,C44;
 712
 713   C00=Cin00-K00*Cin00-K01*Cin10-K02*Cin20-K03*Cin30;
 714   C10=Cin10-K00*Cin10-K01*Cin11-K02*Cin21-K03*Cin31;
 715   C20=Cin20-K00*Cin20-K01*Cin21-K02*Cin22-K03*Cin32;
 716   C30=Cin30-K00*Cin30-K01*Cin31-K02*Cin32-K03*Cin33;
 717   C40=Cin40-K00*Cin40-K01*Cin41-K02*Cin42-K03*Cin43;
 718
 719   C11=Cin11-K10*Cin10-K11*Cin11-K12*Cin21-K13*Cin31;
 720   C21=Cin21-K10*Cin20-K11*Cin21-K12*Cin22-K13*Cin32;
 721   C31=Cin31-K10*Cin30-K11*Cin31-K12*Cin32-K13*Cin33;
 722   C41=Cin41-K10*Cin40-K11*Cin41-K12*Cin42-K13*Cin43;
 723
 724   C22=Cin22-K20*Cin20-K21*Cin21-K22*Cin22-K23*Cin32;
 725   C32=Cin32-K20*Cin30-K21*Cin31-K22*Cin32-K23*Cin33;
 726   C42=Cin42-K20*Cin40-K21*Cin41-K22*Cin42-K23*Cin43;
 727
 728   C33=Cin33-K30*Cin30-K31*Cin31-K32*Cin32-K33*Cin33;
 729   C43=Cin43-K30*Cin40-K31*Cin41-K32*Cin42-K33*Cin43;
 730
 731   C44=Cin44-K40*Cin40-K41*Cin41-K42*Cin42-K43*Cin43;
 732
 733   newTrack->PutXElements(X0,X1,X2,X3,X4);               // put the new state vector
 734
 735   newTrack->PutCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44); // put in track the
 736                                                                                        // new cov matrix
 737
 738   Double_t VMC[4][4];
 739
 740   VMC[0][0]=V00-C00; VMC[1][0]=-C10; VMC[2][0]=-C20; VMC[3][0]=-C30;
 741   VMC[1][1]=V11-C11; VMC[2][1]=-C21; VMC[3][1]=V31-C31;
 742   VMC[2][2]=V22-C22; VMC[3][2]=V32-C32;
 743   VMC[3][3]=V33-C33;
 744
 745   VMC[0][1]=VMC[1][0]; VMC[0][2]=VMC[2][0]; VMC[0][3]=VMC[3][0];
 746   VMC[1][2]=VMC[2][1]; VMC[1][3]=VMC[3][1];
 747   VMC[2][3]=VMC[3][2];
 748
 749
 750 /////////////////////// VMC matrix inversion ///////////////////////////////////
 751
 752   d=1.;
 753
 754   for(k=0; k<n; k++) {
 755     ll[k]=k;
 756     mm[k]=k;
 757     big=VMC[k][k];
 758     for(j=k; j<n ; j++) {
 759       for (i=j; i<n; i++) {
 760         if(TMath::Abs(big) < TMath::Abs(VMC[i][j]) ) { big=VMC[i][j]; ll[k]=i; mm[k]=j; }
 761       }
 762     }
 763 //
 764     j= ll[k];
 765     if(j > k) {
 766       for(i=0; i<n; i++) { hold=-VMC[k][i]; VMC[k][i]=VMC[j][i]; VMC[j][i]=hold; }
 767
 768     }
 769 //
 770     i=mm[k];
 771     if(i > k ) {
 772       for(j=0; j<n; j++) { hold=-VMC[j][k]; VMC[j][k]=VMC[j][i]; VMC[j][i]=hold; }
 773     }
 774 //
 775     if(!big) {
 776       d=0.;
 777       cout << "Singular matrix\n";
 778     }
 779     for(i=0; i<n; i++) {
 780       if(i == k) { continue; }
 781       VMC[i][k]=VMC[i][k]/(-big);
 782     }
 783 //
 784     for(i=0; i<n; i++) {
 785       hold=VMC[i][k];
 786       for(j=0; j<n; j++) {
 787         if(i == k || j == k) { continue; }
 788         VMC[i][j]=hold*VMC[k][j]+VMC[i][j];
 789       }
 790     }
 791 //
 792     for(j=0; j<n; j++) {
 793       if(j == k) { continue; }
 794       VMC[k][j]=VMC[k][j]/big;
 795     }
 796 //
 797     d=d*big;
 798 //
 799     VMC[k][k]=1./big;
 800   }
 801 //
 802   for(k=n-1; k>=0; k--) {
 803     i=ll[k];
 804     if(i > k) {
 805       for (j=0; j<n; j++) {hold=VMC[j][k]; VMC[j][k]=-VMC[j][i]; VMC[j][i]=hold;}
 806     }
 807     j=mm[k];
 808     if(j > k) {
 809       for (i=0; i<n; i++) {hold=VMC[k][i]; VMC[k][i]=-VMC[j][i]; VMC[j][i]=hold;}
 810       }
 811   }
 812
 813
 814 ////////////////////////////////////////////////////////////////////////////////
 815
 816   Double_t chi2=(m[0]-X0)*( VMC[0][0]*(m[0]-X0) + 2.*VMC[1][0]*(m[1]-X1) +
 817                    2.*VMC[2][0]*(m[2]-X2)+ 2.*VMC[3][0]*(m[3]-X3) ) +
 818                 (m[1]-X1)* ( VMC[1][1]*(m[1]-X1) + 2.*VMC[2][1]*(m[2]-X2)+
 819                    2.*VMC[3][1]*(m[3]-X3) ) +
 820                 (m[2]-X2)* ( VMC[2][2]*(m[2]-X2)+ 2.*VMC[3][2]*(m[3]-X3) ) +
 821                 (m[3]-X3)*VMC[3][3]*(m[3]-X3);
 822
 823   newTrack->SetChi2(newTrack->GetChi2()+chi2);
 824
 825 }
 826