ITS/AliITSTrackerV1.cxx

   1 //The purpose of this class is to permorm the ITS tracking.
   2 //The constructor has the task to inizialize some private members.
   3 //The method DoTracking is written to be called by a macro. It gets the event number, the minimum and maximum
   4 //order number of TPC tracks that are to be tracked trough the ITS, and the file where the recpoints are
   5 //registered.
   6 //The method Recursivetracking is a recursive function that performs the tracking trough the ITS
   7 //The method Intersection found the layer, ladder and detector whre the intersection take place and caluclate
   8 //the cohordinates of this intersection. It returns an integer that is 0 if the intersection has been found
   9 //successfully.
  10 //The two mwthods Kalmanfilter and kalmanfiltervert operate the kalmanfilter without and with the vertex
  11 //imposition respectively.
  12 //The authors  thank Mariana Bondila to have help them to resolve some problems.  July-2000
  13
  14
  15 #include <fstream.h>
  16 #include <TMath.h>
  17 #include <TBranch.h>
  18 #include <TVector.h>
  19 #include <TFile.h>
  20 #include <TTree.h>
  21 #include "TParticle.h"
  22 #include "AliRun.h"
  23 #include "AliITS.h"
  24 #include "AliITSsegmentationSSD.h"
  25 #include "AliITSgeomSPD.h"
  26 #include "AliITSgeomSDD.h"
  27 #include "AliITSgeomSSD.h"
  28 #include "AliITSgeom.h"
  29 #include "AliITSRecPoint.h"
  30 #include "stdlib.h"
  31 #include "AliKalmanTrack.h"
  32 #include "AliMagF.h"
  33 #include "AliITSTrackV1.h"
  34 #include "AliITSIOTrack.h"
  35 #include "AliITSRad.h"
  36 #include "../TPC/AliTPCtracker.h"
  37 #include "AliITSTrackerV1.h"
  38
  39
  40 ClassImp(AliITSTrackerV1)
  41
  42
  43 //________________________________________________________________
  44
  45 AliITSTrackerV1::AliITSTrackerV1(AliITS* IITTSS, Bool_t flag) {
  46 //Origin   A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
  47 // Class constructor. It does some initializations.
  48
  49   fITS = IITTSS;
  50   fPtref=0.;
  51   fflagvert=flag;
  52   Int_t imax=200,jmax=450;
  53   frl = new AliITSRad(imax,jmax);
  54
  55 }
  56
  57 AliITSTrackerV1::AliITSTrackerV1(const AliITSTrackerV1 &cobj) {
  58 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
  59 // copy constructor
  60
  61     *fITS = *cobj.fITS;
  62          *fresult = *cobj.fresult;
  63     fPtref = cobj.fPtref;
  64          //frecPoints = fITS->RecPoints();
  65          //*frecPoints = *cobj.frecPoints;
  66          **fvettid = **cobj.fvettid;
  67          fflagvert = cobj.fflagvert;
  68          Int_t imax=200,jmax=450;
  69          frl = new AliITSRad(imax,jmax);
  70          *frl = *cobj.frl;
  71          Int_t i;
  72          for(i=0; i<6; i++) {
  73            fNlad[i] = cobj.fNlad[i];
  74       fNdet[i] = cobj.fNdet[i];
  75            fAvrad[i] = cobj.fAvrad[i];
  76            fDetx[i] = cobj.fDetx[i];
  77            fDetz[i] = cobj.fDetz[i];
  78         }
  79
  80 }
  81
  82 AliITSTrackerV1 &AliITSTrackerV1::operator=(AliITSTrackerV1 obj) {
  83 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
  84 // assignement operator
  85
  86     *fITS = *obj.fITS;
  87          *fresult = *obj.fresult;
  88     fPtref = obj.fPtref;
  89          **fvettid = **obj.fvettid;
  90          fflagvert = obj.fflagvert;
  91          Int_t imax=200,jmax=450;
  92          frl = new AliITSRad(imax,jmax);
  93          *frl = *obj.frl;
  94          Int_t i;
  95          for(i=0; i<6; i++) {
  96            fNlad[i] = obj.fNlad[i];
  97       fNdet[i] = obj.fNdet[i];
  98            fAvrad[i] = obj.fAvrad[i];
  99            fDetx[i] = obj.fDetx[i];
 100            fDetz[i] = obj.fDetz[i];
 101         }
 102
 103   return *this;
 104
 105 }
 106
 107
 108 //________________________________________________________________
 109
 110
 111 void AliITSTrackerV1::DoTracking(Int_t evNumber, Int_t minTr, Int_t maxTr, TFile *file) {
 112 //Origin   A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 113 //The method needs the event number, the minimum and maximum order number of TPC tracks that
 114 //are to be tracked trough the ITS, and the file where the recpoints are registered.
 115 //The method can be called by a macro. It preforms the tracking for all good TPC tracks
 116
 117
 118   printf("begin DoTracking - file %p\n",file);
 119
 120 ///////////////////////////////////////  gets information on geometry ///////////////////////////////////
 121
 122   AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
 123
 124   Int_t ll=1, dd=1;
 125   TVector det(9);
 126
 127   //cout<<" nlad ed ndet \n";
 128   Int_t ia;                              // fuori
 129   for(ia=0; ia<6; ia++) {
 130     fNlad[ia]=g1->GetNladders(ia+1);
 131     fNdet[ia]=g1->GetNdetectors(ia+1);
 132          //cout<<fNlad[i]<<" "<<fNdet[i]<<"\n";
 133   }
 134   //getchar();
 135
 136   //cout<<" raggio medio = ";
 137   Int_t ib;                                  // fuori
 138   for(ib=0; ib<6; ib++) {
 139     g1->GetCenterThetaPhi(ib+1,ll,dd,det);
 140     fAvrad[ib]=TMath::Sqrt(det(0)*det(0)+det(1)*det(1));
 141          //cout<<fAvrad[ib]<<" ";
 142   }
 143   //cout<<"\n"; getchar();
 144
 145   fDetx[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDx();
 146   fDetz[0] = ((AliITSgeomSPD*)(g1->GetShape(1, ll, dd)))->GetDz();
 147
 148   fDetx[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDx();
 149   fDetz[1] = ((AliITSgeomSPD*)(g1->GetShape(2, ll, dd)))->GetDz();
 150
 151   fDetx[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDx();
 152   fDetz[2] = ((AliITSgeomSDD*)(g1->GetShape(3, ll, dd)))->GetDz();
 153
 154   fDetx[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDx();
 155   fDetz[3] = ((AliITSgeomSDD*)(g1->GetShape(4, ll, dd)))->GetDz();
 156
 157   fDetx[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDx();
 158   fDetz[4] = ((AliITSgeomSSD*)(g1->GetShape(5, ll, dd)))->GetDz();
 159
 160   fDetx[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDx();
 161   fDetz[5] = ((AliITSgeomSSD*)(g1->GetShape(6, ll, dd)))->GetDz();
 162
 163   //cout<<"    Detx     Detz\n";
 164   //for(Int_t la=0; la<6; la++) cout<<"    "<<fDetx[la]<<"     "<<fDetz[la]<<"\n";
 165   //getchar();
 166 //////////////////////////////////////////////////////////////////////////////////////////////////////////
 167
 168
 169  //Int_t imax=200,jmax=450;
 170  //frl = new AliITSRad(imax,jmax);
 171  //cout<<" dopo costruttore AliITSRad\n"; getchar();
 172
 173   struct GoodTrack {
 174     Int_t lab,code;
 175     Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
 176     Bool_t flag;
 177   };
 178
 179
 180   gAlice->GetEvent(0);
 181
 182     AliKalmanTrack *kkprov;
 183     kkprov->SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
 184
 185
 186   TFile *cf=TFile::Open("AliTPCclusters.root");
 187   AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
 188   if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
 189
 190    AliTPCtracker *tracker = new AliTPCtracker(digp);
 191
 192  // Load clusters
 193    tracker->LoadInnerSectors();
 194    tracker->LoadOuterSectors();
 195
 196
 197   GoodTrack gt[15000];
 198   Int_t ngood=0;
 199   ifstream in("itsgood_tracks");
 200
 201   cerr<<"Reading itsgood tracks...\n";
 202   while (in>>gt[ngood].lab>>gt[ngood].code
 203           >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
 204           >>gt[ngood].x  >>gt[ngood].y >>gt[ngood].z
 205           >>gt[ngood].pxg  >>gt[ngood].pyg >>gt[ngood].pzg
 206           >>gt[ngood].ptg >>gt[ngood].flag) {
 207     ngood++;
 208     cerr<<ngood<<'\r';
 209     if (ngood==15000) {
 210       cerr<<"Too many good tracks !\n";
 211       break;
 212     }
 213   }
 214   if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
 215
 216
 217 // Load tracks
 218   TFile *tf=TFile::Open("AliTPCtracks.root");
 219   if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
 220   TObjArray tracks(200000);
 221    TTree *tracktree=(TTree*)tf->Get("TPCf");
 222    if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n";}
 223   TBranch *tbranch=tracktree->GetBranch("tracks");
 224   Int_t nentr=(Int_t)tracktree->GetEntries();
 225   Int_t kk;
 226
 227    AliTPCtrack *ioTrackTPC=0;
 228   for (kk=0; kk<nentr; kk++) {
 229     ioTrackTPC=new AliTPCtrack;
 230     tbranch->SetAddress(&ioTrackTPC);
 231     tracktree->GetEvent(kk);
 232     tracker->CookLabel(ioTrackTPC,0.1);
 233     tracks.AddLast(ioTrackTPC);
 234   }
 235    delete tracker;
 236   tf->Close();
 237
 238
 239   Int_t nt = tracks.GetEntriesFast();
 240   cerr<<"Number of found tracks "<<nt<<endl;
 241
 242   TVector dataOut(9);
 243   Int_t kkk=0;
 244
 245   Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
 246
 247   //////////////////////////////  good tracks definition in TPC  ////////////////////////////////
 248
 249   ofstream out1 ("AliITSTrag.out");
 250   Int_t i;
 251   for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
 252   out1.close();
 253
 254
 255   TVector vec(5);
 256   TTree *tr=gAlice->TreeR();
 257   Int_t nent=(Int_t)tr->GetEntries();
 258   //TClonesArray  *recPoints = RecPoints();      // nuova eliminata
 259   //TClonesArray  *recPoints = ITS->RecPoints();  // nuova
 260   //TObjArray  *
 261   frecPoints = fITS->RecPoints();  // nuovissima   tolta
 262
 263   Int_t numbpoints;
 264   Int_t totalpoints=0;
 265   Int_t *np = new Int_t[nent];
 266   //Int_t **
 267   fvettid = new Int_t* [nent];
 268   Int_t mod;
 269
 270   for (mod=0; mod<nent; mod++) {
 271     fvettid[mod]=0;
 272     fITS->ResetRecPoints();  // nuova
 273     //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
 274     gAlice->TreeR()->GetEvent(mod); //first entry in TreeR is empty
 275     numbpoints = frecPoints->GetEntries();
 276     totalpoints+=numbpoints;
 277     np[mod] = numbpoints;
 278   //cout<<" mod = "<<mod<<"   numbpoints = "<<numbpoints<<"\n"; getchar();
 279     fvettid[mod] = new Int_t[numbpoints];
 280     Int_t ii;
 281     for (ii=0;ii<numbpoints; ii++) *(fvettid[mod]+ii)=0;
 282   }
 283
 284   AliTPCtrack *track=0;  // sono qui
 285
 286
 287   if(minTr < 0) {minTr = 0; maxTr = nt-1;}
 288
 289 /*
 290   ///////////////////////////////// Definition of vertex end its error ////////////////////////////
 291   ////////////////////////// In the future it will be given by a method ///////////////////////////
 292   Double_t Vx=0.;
 293   Double_t Vy=0.;
 294   Double_t Vz=0.;
 295
 296   Float_t sigmavx=0.0050;      // 50  microns
 297   Float_t sigmavy=0.0050;      // 50  microns
 298   Float_t sigmavz=0.010;       // 100 microns
 299
 300   //Vx+=gRandom->Gaus(0,sigmavx);  Vy+=gRandom->Gaus(0,sigmavy);  Vz+=gRandom->Gaus(0,sigmavz);
 301   TVector vertex(3), ervertex(3)
 302   vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
 303   ervertex(0)=sigmavx;  ervertex(1)=sigmavy;  ervertex(2)=sigmavz;
 304   /////////////////////////////////////////////////////////////////////////////////////////////////
 305 */
 306
 307
 308   TTree tracktree1("TreeT","Tree with ITS tracks");
 309   AliITSIOTrack *ioTrack=0;
 310   tracktree1.Branch("ITStracks","AliITSIOTrack",&ioTrack,32000,0);
 311
 312   ofstream out ("AliITSTra.out");
 313
 314
 315   Int_t j;
 316   for (j=minTr; j<=maxTr; j++) {
 317     track=(AliTPCtrack*)tracks.UncheckedAt(j);
 318     Int_t flaglab=0;
 319     if (!track) continue;
 320     ////// elimination of not good tracks ////////////
 321     Int_t ilab=TMath::Abs(track->GetLabel());
 322     Int_t iii;
 323     for (iii=0;iii<ngood;iii++) {
 324          //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n"; getchar();
 325       if (ilab==gt[iii].lab) {
 326         flaglab=1;
 327         ptg=gt[iii].ptg;
 328         pxg=gt[iii].pxg;
 329         pyg=gt[iii].pyg;
 330         pzg=gt[iii].pzg;
 331         break;
 332       }
 333     }
 334          //cout<<" j flaglab =  " <<j<<" "<<flaglab<<"\n";  getchar();
 335     if (!flaglab) continue;
 336          //cout<<" j =  " <<j<<"\n";  getchar();
 337
 338          //////   propagation to the end of TPC //////////////
 339     Double_t xk=77.415;
 340     track->PropagateTo(xk, 28.94, 1.204e-3);     //Ne
 341          xk -=0.01;
 342     track->PropagateTo(xk, 44.77, 1.71);         //Tedlar
 343          xk -=0.04;
 344     track->PropagateTo(xk, 44.86, 1.45);         //kevlar
 345          xk -=2.0;
 346     track->PropagateTo(xk, 41.28, 0.029);        //Nomex
 347     xk-=16;
 348     track->PropagateTo(xk,36.2,1.98e-3); //C02
 349          xk -=0.01;
 350     track->PropagateTo(xk, 24.01, 2.7);  //Al
 351          xk -=0.01;
 352     track->PropagateTo(xk, 44.77, 1.71);         //Tedlar
 353          xk -=0.04;
 354     track->PropagateTo(xk, 44.86, 1.45);         //kevlar
 355          xk -=0.5;
 356     track->PropagateTo(xk, 41.28, 0.029);        //Nomex
 357
 358     ///////////////////////////////////////////////////////////////
 359
 360
 361          AliITSTrackV1 trackITS(*track);
 362
 363          if(fresult) delete fresult;
 364          fresult = new AliITSTrackV1(trackITS);
 365
 366          AliITSTrackV1 primaryTrack(trackITS);
 367
 368
 369          TVector vgeant(3);
 370          vgeant=(*fresult).GetVertex();
 371
 372   // Definition of dv and zv for vertex constraint
 373      Double_t sigmaDv=0.0050;  Double_t sigmaZv=0.010;
 374     //Double_t sigmaDv=0.0015;  Double_t sigmaZv=0.0015;
 375         Double_t uniform= gRandom->Uniform();
 376         Double_t signdv;
 377         if(uniform<=0.5) signdv=-1.;
 378            else
 379                  signdv=1.;
 380
 381         Double_t vr=TMath::Sqrt(vgeant(0)*vgeant(0)+ vgeant(1)*vgeant(1));
 382           Double_t dv=gRandom->Gaus(signdv*vr,(Float_t)sigmaDv);
 383     Double_t zv=gRandom->Gaus(vgeant(2),(Float_t)sigmaZv);
 384
 385   //cout<<" Dv e Zv = "<<dv<<" "<<zv<<"\n";
 386     trackITS.SetDv(dv);  trackITS.SetZv(zv);
 387     trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv);
 388     (*fresult).SetDv(dv);  (*fresult).SetZv(zv);
 389     (*fresult).SetsigmaDv(sigmaDv); (*fresult).SetsigmaZv(sigmaZv);
 390     primaryTrack.SetDv(dv);  primaryTrack.SetZv(zv);
 391     primaryTrack.SetsigmaDv(sigmaDv); primaryTrack.SetsigmaZv(sigmaZv);
 392
 393     primaryTrack.PrimaryTrack(frl);
 394     TVector  d2=primaryTrack.Getd2();
 395     TVector  tgl2=primaryTrack.Gettgl2();
 396     TVector  dtgl=primaryTrack.Getdtgl();
 397     trackITS.Setd2(d2); trackITS.Settgl2(tgl2);  trackITS.Setdtgl(dtgl);
 398     (*fresult).Setd2(d2); (*fresult).Settgl2(tgl2);  (*fresult).Setdtgl(dtgl);
 399          /*
 400     trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
 401     (*result).SetVertex(vertex);   (*result).SetErrorVertex(ervertex);
 402     */
 403
 404
 405   TList *list= new TList();
 406
 407   list->AddLast(&trackITS);
 408
 409   fPtref=TMath::Abs( (trackITS).GetPt() );
 410   //cout << "\n Pt = " << fPtref <<"\n";  //stampa
 411
 412   RecursiveTracking(list);  // nuova ITS
 413   list->Delete();
 414   delete list;
 415
 416   Int_t itot=-1;
 417   TVector vecTotLabRef(18);
 418   Int_t lay, k;
 419   for(lay=5; lay>=0; lay--) {
 420     TVector vecLabRef(3);
 421     vecLabRef=(*fresult).GetLabTrack(lay);
 422     Float_t clustZ=(*fresult).GetZclusterTrack( lay);
 423     for(k=0; k<3; k++){
 424                 Int_t lpp=(Int_t)vecLabRef(k);
 425                 if(lpp>=0) {
 426                   TParticle *p=(TParticle*) gAlice->Particle(lpp);
 427                   Int_t pcode=p->GetPdgCode();
 428                   if(pcode==11) vecLabRef(k)=p->GetFirstMother();
 429                 }
 430     itot++; vecTotLabRef(itot)=vecLabRef(k);
 431     if(vecLabRef(k)==0. && clustZ == 0.) vecTotLabRef(itot) =-3.; }
 432   }
 433   Long_t labref;
 434   Int_t freq;
 435   (*fresult).Search(vecTotLabRef, labref, freq);
 436
 437
 438   //if(freq < 6) labref=-labref;        // cinque - sei
 439   if(freq < 5) labref=-labref;        // cinque - sei
 440   (*fresult).SetLabel(labref);
 441
 442     // cout<<" progressive track number = "<<j<<"\r";
 443    // cout<<j<<"\r";
 444     Int_t numOfCluster=(*fresult).GetNumClust();
 445     //cout<<" progressive track number = "<<j<<"\n";    // stampa
 446     Long_t labITS=(*fresult).GetLabel();
 447     //cout << " ITS track label = " << labITS << "\n";  // stampa
 448     int lab=track->GetLabel();
 449     //cout << " TPC track label = " << lab <<"\n";      // stampa
 450
 451
 452 //propagation to vertex
 453
 454     Double_t rbeam=3.;
 455
 456     (*fresult).Propagation(rbeam);
 457
 458          Double_t c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44;
 459          (*fresult).GetCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44);
 460
 461     Double_t pt=TMath::Abs((*fresult).GetPt());
 462     Double_t dr=(*fresult).GetD();
 463     Double_t z=(*fresult).GetZ();
 464     Double_t tgl=(*fresult).GetTgl();
 465     Double_t c=(*fresult).GetC();
 466     Double_t cy=c/2.;
 467     Double_t dz=z-(tgl/cy)*TMath::ASin((*fresult).Arga(rbeam));
 468          dz-=vgeant(2);
 469
 470          // cout<<" dr e dz alla fine = "<<dr<<" "<<dz<<"\n"; getchar();
 471     Double_t phi=(*fresult).Getphi();
 472     Double_t phivertex = phi - TMath::ASin((*fresult).ArgA(rbeam));
 473     Double_t duepi=2.*TMath::Pi();
 474     if(phivertex>duepi) phivertex-=duepi;
 475     if(phivertex<0.) phivertex+=duepi;
 476     Double_t dtot=TMath::Sqrt(dr*dr+dz*dz);
 477
 478 //////////////////////////////////////////////////////////////////////////////////////////
 479
 480     Int_t idmodule,idpoint;
 481          if(numOfCluster >=5)  {            // cinque - sei
 482          //if(numOfCluster ==6)  {            // cinque - sei
 483
 484
 485       AliITSIOTrack outTrack;
 486
 487       ioTrack=&outTrack;
 488
 489       ioTrack->SetStatePhi(phi);
 490       ioTrack->SetStateZ(z);
 491       ioTrack->SetStateD(dr);
 492       ioTrack->SetStateTgl(tgl);
 493       ioTrack->SetStateC(c);
 494                 Double_t radius=(*fresult).Getrtrack();
 495                 ioTrack->SetRadius(radius);
 496                 Int_t charge;
 497                 if(c>0.) charge=-1;  else charge=1;
 498                 ioTrack->SetCharge(charge);
 499
 500
 501
 502       ioTrack->SetCovMatrix(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44);
 503
 504       Double_t px=pt*TMath::Cos(phivertex);
 505       Double_t py=pt*TMath::Sin(phivertex);
 506       Double_t pz=pt*tgl;
 507
 508       Double_t xtrack=dr*TMath::Sin(phivertex);
 509       Double_t ytrack=dr*TMath::Cos(phivertex);
 510       Double_t ztrack=dz+vgeant(2);
 511
 512
 513       ioTrack->SetPx(px);
 514       ioTrack->SetPy(py);
 515       ioTrack->SetPz(pz);
 516       ioTrack->SetX(xtrack);
 517       ioTrack->SetY(ytrack);
 518       ioTrack->SetZ(ztrack);
 519       ioTrack->SetLabel(labITS);
 520
 521       Int_t il;
 522                 for(il=0;il<6; il++){
 523                   ioTrack->SetIdPoint(il,(*fresult).GetIdPoint(il));
 524                   ioTrack->SetIdModule(il,(*fresult).GetIdModule(il));
 525                 }
 526       tracktree1.Fill();
 527
 528    //cout<<" labITS = "<<labITS<<"\n";
 529         //cout<<" phi z dr tgl c = "<<phi<<" "<<z<<" "<<dr<<" "<<tgl<<" "<<c<<"\n";  getchar();
 530
 531      dataOut(kkk) = ptg; kkk++; dataOut(kkk)=labITS; kkk++; dataOut(kkk)=lab; kkk++;
 532
 533       for (il=0;il<6;il++) {
 534         idpoint=(*fresult).GetIdPoint(il);
 535         idmodule=(*fresult).GetIdModule(il);
 536         *(fvettid[idmodule]+idpoint)=1;
 537         ioTrack->SetIdPoint(il,idpoint);
 538         ioTrack->SetIdModule(il,idmodule);
 539       }
 540
 541     //  cout<<"  +++++++++++++  pt e ptg = "<<pt<<" "<<ptg<<"  ++++++++++\n";
 542
 543        ///////////////////////////////
 544       Double_t difpt= (pt-ptg)/ptg*100.;
 545       dataOut(kkk)=difpt; kkk++;
 546       Double_t lambdag=TMath::ATan(pzg/ptg);
 547       Double_t   lam=TMath::ATan(tgl);
 548       Double_t diflam = (lam - lambdag)*1000.;
 549       dataOut(kkk) = diflam; kkk++;
 550       Double_t phig=TMath::ATan2(pyg,pxg);  if(phig<0) phig=2.*TMath::Pi()+phig;
 551       Double_t phi=phivertex;
 552
 553       Double_t difphi = (phi - phig)*1000.;
 554       dataOut(kkk)=difphi; kkk++;
 555       dataOut(kkk)=dtot*1.e4; kkk++;
 556       dataOut(kkk)=dr*1.e4; kkk++;
 557       dataOut(kkk)=dz*1.e4; kkk++;
 558       Int_t r;
 559       for (r=0; r<9; r++) { out<<dataOut(r)<<" ";}
 560       out<<"\n";
 561       kkk=0;
 562
 563
 564     } // end if on numOfCluster
 565   //gObjectTable->Print();    // stampa memoria
 566   }  //  end for (int j=minTr; j<=maxTr; j++)
 567
 568   out.close();
 569
 570
 571   static Bool_t first=kTRUE;
 572   static TFile *tfile;
 573
 574         if(first) {
 575             tfile=new TFile("itstracks.root","RECREATE");
 576             //cout<<"I have opened itstracks.root file "<<endl;
 577         }
 578         first=kFALSE;
 579         tfile->cd();
 580         tfile->ls();
 581
 582    char hname[30];
 583    sprintf(hname,"TreeT%d",evNumber);
 584
 585   tracktree1.Write(hname);
 586
 587
 588
 589             TTree *fAli=gAlice->TreeK();
 590             TFile *fileAli=0;
 591
 592             if (fAli) fileAli =fAli->GetCurrentFile();
 593             fileAli->cd();
 594
 595   ////////////////////////////////////////////////////////////////////////////////////////////////
 596
 597   printf("delete vectors\n");
 598   if(np) delete [] np;
 599   if(fvettid) delete [] fvettid;
 600   if(fresult) delete fresult;
 601
 602 }
 603
 604
 605
 606 void AliITSTrackerV1::RecursiveTracking(TList *trackITSlist) {
 607
 608 ///////////////////////   This function perform the recursive tracking in ITS detectors /////////////////////
 609 ///////////////////////     reference is a pointer to the final best track    /////////////////////
 610 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 611 // The authors  thank   Mariana Bondila to have help them to resolve some problems.  July-2000
 612
 613   //Rlayer[0]=4.; Rlayer[1]=7.;  Rlayer[2]=14.9;  Rlayer[3]=23.8;  Rlayer[4]=39.1;  Rlayer[5]=43.6; //vecchio
 614
 615   Int_t index;
 616   for(index =0; index<trackITSlist->GetSize(); index++) {
 617     AliITSTrackV1 *trackITS = (AliITSTrackV1 *) trackITSlist->At(index);
 618
 619     if((*trackITS).GetLayer()==7) fresult->SetChi2(10.223e140);
 620    // cout <<" Layer inizio = "<<(*trackITS).GetLayer()<<"\n";
 621    //  cout<<"fvtrack =" <<"\n";
 622    //  cout << (*trackITS)(0) << " "<<(*trackITS)(1)<<" "<<(*trackITS)(2)<<" "<<(*trackITS)(3)<<" "<<(*trackITS)(4)<<"\n";
 623    //  cout<< " rtrack = "<<(*trackITS).Getrtrack()<<"\n";
 624    //  cout<< " Pt = "<<(*trackITS).GetPt()<<"\n";
 625    //  getchar();
 626     Double_t chi2Now, chi2Ref;
 627     if((*trackITS).GetLayer()==1 ) {
 628       chi2Now = trackITS->GetChi2();
 629       Float_t numClustNow = trackITS->GetNumClust();
 630       if(trackITS->GetNumClust()) chi2Now /= (Double_t )trackITS->GetNumClust();
 631       chi2Ref = fresult->GetChi2();
 632       Float_t numClustRef = fresult->GetNumClust();
 633       if(fresult->GetNumClust()) chi2Ref /= (Double_t )fresult->GetNumClust();
 634       //cout<<" chi2Now and chi2Ref = "<<chi2Now<<" "<<chi2Ref<<"\n";
 635                 if( numClustNow > numClustRef ) {*fresult = *trackITS;}
 636       if((numClustNow == numClustRef )&& (chi2Now < chi2Ref))  {*fresult = *trackITS;}
 637       continue;
 638     }
 639     Float_t numClustNow = trackITS->GetNumClust();
 640     if(numClustNow) {
 641       chi2Now = trackITS->GetChi2();
 642       chi2Now/=numClustNow;
 643       //cout<<" chi2Now =  "<<chi2Now<<"\n";
 644     /*
 645     // if(Ptref > 0.6 && chi2Now > 20.) continue;
 646     if(Ptref > 0.6 && chi2Now > 30.) continue;
 647     if((Ptref <= 0.6 && Ptref>0.2)&& chi2Now > 15.) continue;
 648      // if(chi2Now>5.) continue;
 649       //if(chi2Now>15.) continue;
 650      // if(Ptref <= 0.2 && chi2Now > 10.) continue;
 651      if(Ptref <= 0.2 && chi2Now > 8.) continue;
 652      */
 653      if(fPtref > 1.0 && chi2Now > 30.) continue;
 654      if((fPtref >= 0.6 && fPtref<=1.0) && chi2Now > 40.) continue;
 655      if((fPtref <= 0.6 && fPtref>0.2)&& chi2Now > 40.) continue;
 656      if(fPtref <= 0.2 && chi2Now > 8.) continue;
 657     }
 658
 659     Int_t layerInit = (*trackITS).GetLayer();
 660     Int_t layernew = layerInit - 2;  // -1 for new layer, -1 for matrix index
 661
 662     //Int_t NLadder[]= {20, 40, 14, 22, 34, 38};   //vecchio
 663     //Int_t NDetector[]= {4,  4,   6,  8, 23, 26}; //vecchio
 664
 665     TList listoftrack;
 666     Int_t ladp, ladm, detp,detm,ladinters,detinters;
 667     Int_t layerfin=layerInit-1;
 668          Double_t rFin=fAvrad[layerfin-1];
 669     // cout<<"Prima di intersection \n";
 670
 671     Int_t  outinters=Intersection(*trackITS, rFin, layerfin, ladinters, detinters);
 672
 673    // cout<<" outinters = "<<outinters<<"\n";
 674    //  cout<<" Layer ladder detector intersection ="<<layerfin<<" "<<ladinters<<" "<<detinters<<"\n";
 675    //  cout << " phiinters zinters = "<<(*trackITS)(0) << " "<<(*trackITS)(1)<<"\n"; getchar();
 676
 677     if(outinters==-1) continue;
 678
 679     Int_t flaghit=0;
 680     if(outinters==0){
 681       TVector toucLad(9), toucDet(9);
 682       Int_t lycur=layerfin;
 683       ladp=ladinters+1;
 684       ladm=ladinters-1;
 685                 if(ladm <= 0) ladm=fNlad[layerfin-1];
 686                 if(ladp > fNlad[layerfin-1]) ladp=1;
 687       detp=detinters+1;
 688       detm=detinters-1;
 689       Int_t idetot=1;
 690       toucLad(0)=ladinters; toucLad(1)=ladm; toucLad(2)=ladp;
 691       toucLad(3)=ladinters; toucLad(4)=ladm; toucLad(5)=ladp;
 692       toucLad(6)=ladinters; toucLad(7)=ladm; toucLad(8)=ladp;
 693       toucDet(0)=detinters; toucDet(1)=detinters; toucDet(2)=detinters;
 694                 if(detm > 0 && detp <= fNdet[layerfin-1]) {
 695         idetot=9;
 696         toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm;
 697         toucDet(6)=detp; toucDet(7)=detp; toucDet(8)=detp;
 698       }
 699
 700                 if(detm > 0 && detp > fNdet[layerfin-1]) {
 701         idetot=6;
 702         toucDet(3)=detm; toucDet(4)=detm; toucDet(5)=detm;
 703       }
 704
 705                 if(detm <= 0 && detp <= fNdet[layerfin-1]) {
 706         idetot=6;
 707         toucDet(3)=detp; toucDet(4)=detp; toucDet(5)=detp;
 708       }
 709       Int_t iriv;
 710       for (iriv=0; iriv<idetot; iriv++) {  //for on detectors
 711         //AliITSgeom *g1 = aliITS->GetITSgeom();  //vvecchia
 712                   AliITSgeom *g1 = fITS->GetITSgeom();  //nnuova
 713         TVector ctf(9);
 714         g1->GetCenterThetaPhi(layerInit-1,(Int_t)toucLad(iriv),(Int_t)toucDet(iriv),ctf);
 715
 716         // cout<<" layer, ladder, det, xo, yo, zo = "<<layerInit-1<<" "<<(Int_t)toucLad(iriv)<<
 717         // " "<<(Int_t)toucDet(iriv)<<" "<<ctf(0)<<" "<<ctf(1)<<" "<<ctf(2)<< " "<<ctf(6)<<"\n"; getchar();
 718
 719         ////////////////////////////////////////////////////////////////////////////////////////////////
 720
 721         /*** Rec points sorted by module *****/
 722         /**************************************/
 723
 724         Int_t index;
 725         AliITSRecPoint *recp;
 726         //AliITSgeom *geom = aliITS->GetITSgeom();  //vvecchia
 727                   AliITSgeom *geom = fITS->GetITSgeom();  //nnuova
 728         index = geom->GetModuleIndex(lycur,toucLad(iriv),toucDet(iriv));
 729         Int_t lay,lad,det;
 730         geom->GetModuleId(index,lay,lad,det);
 731         //aliITS->ResetRecPoints();   //vvecchia
 732                   fITS->ResetRecPoints();   //nnuova
 733         //gAlice->TreeR()->GetEvent(index+1); //first entry in TreeR is empty
 734         gAlice->TreeR()->GetEvent(index); //first entry in TreeR is empty
 735
 736         Int_t npoints=frecPoints->GetEntries();
 737         Int_t *indlist=new Int_t[npoints+1];
 738         Int_t counter=0;
 739         Int_t ind;
 740         for (ind=0; ind<=npoints; ind++) {
 741           indlist[ind]=-1;
 742                if (*(fvettid[index]+ind)==0) {
 743               indlist[counter]=ind;
 744                    counter++;
 745               }
 746         }
 747
 748         ind=-1;
 749
 750         for(;;) {
 751           ind++;
 752           if(indlist[ind] < 0) recp=0;
 753                else recp = (AliITSRecPoint*)frecPoints->UncheckedAt(indlist[ind]);
 754
 755           if((!recp)  )  break;
 756           TVector cluster(3),vecclust(9);
 757           vecclust(6)=vecclust(7)=vecclust(8)=-1.;
 758           Double_t sigma[2];
 759           // set veclust in global
 760           Float_t global[3], local[3];
 761           local[0]=recp->GetX();
 762           local[1]=0.;
 763           local[2]= recp->GetZ();
 764           //AliITSgeom *g1 = aliITS->GetITSgeom();  //vvecchia
 765                          AliITSgeom *g1 = fITS->GetITSgeom();  //nnuova
 766           Int_t play = lycur;
 767           Int_t plad = TMath::Nint(toucLad(iriv));
 768           Int_t pdet = TMath::Nint(toucDet(iriv));
 769           g1->LtoG(play,plad,pdet,local,global);
 770
 771           vecclust(0)=global[0];
 772           vecclust(1)=global[1];
 773           vecclust(2)=global[2];
 774
 775
 776           vecclust(3) = (float)recp->fTracks[0];
 777           vecclust(4) = (float)indlist[ind];
 778           vecclust(5) = (float)index;
 779           vecclust(6) = (float)recp->fTracks[0];
 780           vecclust(7) = (float)recp->fTracks[1];
 781           vecclust(8) = (float)recp->fTracks[2];
 782
 783           sigma[0] = (Double_t)  recp->GetSigmaX2();
 784           sigma[1] = (Double_t) recp->GetSigmaZ2();
 785
 786          //now we are in r,phi,z in global
 787           cluster(0) = TMath::Sqrt(vecclust(0)*vecclust(0)+vecclust(1)*vecclust(1));//r hit
 788          // cluster(1) = PhiDef(vecclust(0),vecclust(1));    // phi hit  //vecchio
 789                          cluster(1) = TMath::ATan2(vecclust(1),vecclust(0)); if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi();  //nuovo
 790           cluster(2) = vecclust(2);                   // z hit
 791
 792          // cout<<" layer = "<<play<<"\n";
 793          // cout<<" cluster prima = "<<vecclust(0)<<" "<<vecclust(1)<<" "
 794          // <<vecclust(2)<<"\n"; getchar();
 795           //cluster(1)= cluster(1)-trackITS->Getalphaprov();  //provvisorio;
 796                          //if(cluster(1)<0.) cluster(1)+=2.*TMath::Pi(); //provvisorio
 797                          //cout<<" cluster(1) dopo = "<<cluster(1)<< " alphaprov = "<<trackITS->Getalphaprov()<<"\n";
 798           Float_t sigmatotphi, sigmatotz;
 799
 800           //Float_t epsphi=3.2, epsz=3.;
 801                Float_t epsphi=5.0, epsz=5.0;
 802           if(fPtref<0.2) {epsphi=3.; epsz=3.;}
 803
 804           Double_t rTrack=(*trackITS).Getrtrack();
 805           Double_t sigmaphi=sigma[0]/(rTrack*rTrack);
 806           sigmatotphi=epsphi*TMath::Sqrt(sigmaphi + (*trackITS).GetSigmaphi());
 807
 808           sigmatotz=epsz*TMath::Sqrt(sigma[1] + (*trackITS).GetSigmaZ());
 809   //cout<<"cluster e sigmatotphi e track = "<<cluster(0)<<" "<<cluster(1)<<" "<<sigmatotphi<<" "<<vecclust(3)<<"\n";
 810   //if(vecclust(3)==481) getchar();
 811                if(cluster(1)<6. && (*trackITS).Getphi()>6.) cluster(1)=cluster(1)+(2.*TMath::Pi());
 812                if(cluster(1)>6. && (*trackITS).Getphi()<6.) cluster(1)=cluster(1)-(2.*TMath::Pi());
 813           if(TMath::Abs(cluster(1)-(*trackITS).Getphi()) > sigmatotphi) continue;
 814                         // cout<<" supero sigmaphi \n";
 815           AliITSTrackV1 *newTrack = new AliITSTrackV1((*trackITS));
 816           (*newTrack).SetLayer((*trackITS).GetLayer()-1);
 817
 818                if (TMath::Abs(rTrack-cluster(0))/rTrack>1e-6)
 819                                                (*newTrack).Correct(Double_t(cluster(0)));
 820                         //cout<<" cluster(2) e (*newTrack).GetZ() = "<<cluster(2)<<" "<<        (*newTrack).GetZ()<<"\n";
 821           if(TMath::Abs(cluster(2)-(*newTrack).GetZ()) > sigmatotz){
 822              delete newTrack;
 823              continue;}
 824
 825           if(iriv == 0) flaghit=1;
 826
 827           (*newTrack).AddMS(frl);  // add the multiple scattering matrix to the covariance matrix
 828           (*newTrack).AddEL(frl,1.,0);
 829
 830           Double_t sigmanew[2];
 831           sigmanew[0]= sigmaphi;
 832           sigmanew[1]=sigma[1];
 833
 834           if(fflagvert)
 835             KalmanFilterVert(newTrack,cluster,sigmanew);
 836           else
 837             KalmanFilter(newTrack,cluster,sigmanew);
 838
 839
 840           (*newTrack).PutCluster(layernew, vecclust);
 841            newTrack->AddClustInTrack();
 842
 843            listoftrack.AddLast(newTrack);
 844
 845         }   // end of for(;;) on rec points
 846
 847         delete [] indlist;
 848
 849       }  // end of for on detectors
 850
 851     }//end if(outinters==0)
 852
 853     if(flaghit==0 || outinters==-2) {
 854       AliITSTrackV1 *newTrack = new AliITSTrackV1(*trackITS);
 855       (*newTrack).SetLayer((*trackITS).GetLayer()-1);
 856       (*newTrack).AddMS(frl);  // add the multiple scattering matrix to the covariance matrix
 857       (*newTrack).AddEL(frl,1.,0);
 858
 859       listoftrack.AddLast(newTrack);
 860     }
 861
 862
 863     //gObjectTable->Print();   // stampa memoria
 864
 865     RecursiveTracking(&listoftrack);
 866     listoftrack.Delete();
 867   } // end of for on tracks
 868
 869   //gObjectTable->Print();   // stampa memoria
 870
 871 }
 872
 873 Int_t AliITSTrackerV1::Intersection(AliITSTrackV1 &track, Double_t rk,Int_t layer, Int_t &ladder,
 874 Int_t &detector) {
 875 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 876 // Found the intersection and the detector
 877
 878   if(track.DoNotCross(rk)){ /*cout<< " Do not cross \n";*/ return -1;}
 879   track.Propagation(rk);
 880   Double_t zinters=track.GetZ();
 881   Double_t phinters=track.Getphi();
 882   //cout<<"zinters = "<<zinters<<"  phinters = "<<phinters<<"\n";
 883
 884   //////////////////////////////////      limits for Geometry 5      /////////////////////////////
 885
 886   //Int_t NLadder[]= {20, 40, 14, 22, 34, 38};
 887   //Int_t NDetector[]= {4,  4,   6,  8, 23, 26};
 888
 889   //Float_t Detx[]= {0.64, 0.64, 3.509, 3.509, 3.65, 3.65 };
 890   //Float_t Detz[]= {4.19, 4.19, 3.75 , 3.75 , 2   , 2    };
 891
 892   ////////////////////////////////////////////////////////////////////////////////////////////////
 893
 894   TVector det(9);
 895   TVector listDet(2);
 896   TVector distZCenter(2);
 897   AliITSgeom *g1 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
 898
 899   Int_t iz=0;
 900   Double_t epsz=1.2;
 901   Double_t epszpixel=0.05;
 902
 903   Int_t iD;
 904   for(iD = 1; iD<= fNdet[layer-1]; iD++) {
 905     g1->GetCenterThetaPhi(layer,1,iD,det);
 906          Double_t zmin=det(2)-fDetz[layer-1];
 907          if(iD==1) zmin=det(2)-(fDetz[layer-1])*epsz;
 908          Double_t zmax=det(2)+fDetz[layer-1];
 909          if(iD==fNdet[layer-1]) zmax=det(2)+(fDetz[layer-1])*epsz;
 910     //added to take into account problem on drift
 911     if(layer == 4 || layer==3) zmin=zmin-epszpixel; zmax=zmax+epszpixel;
 912     //cout<<"zmin zinters zmax det(2)= "<<zmin<<" "<<zinters<<" "<<zmax<<" "<<det(2)<<"\n";
 913     if(zinters > zmin && zinters <= zmax) {
 914       if(iz>1) {cout<< " Errore su iz in Intersection \n"; getchar();}
 915       else {
 916         listDet(iz)= iD; distZCenter(iz)=TMath::Abs(zinters-det(2)); iz++;
 917       }
 918     }
 919   }
 920
 921   if(iz==0) {/* cout<< " No detector along Z \n";*/ return -2;}
 922   detector=Int_t (listDet(0));
 923   if(iz>1 && (distZCenter(0)>distZCenter(1)))   detector=Int_t (listDet(1));
 924
 925   AliITSgeom *g2 = ((AliITS*)gAlice->GetDetector("ITS"))->GetITSgeom();
 926   Float_t global[3];
 927   Float_t local[3];
 928   TVector listLad(2);
 929   TVector distPhiCenter(2);
 930   Int_t ip=0;
 931   Double_t pigre=TMath::Pi();
 932
 933   Int_t iLd;
 934   for(iLd = 1; iLd<= fNlad[layer-1]; iLd++) {
 935           g1->GetCenterThetaPhi(layer,iLd,detector,det);
 936  // Double_t phidet=PhiDef(Double_t(det(0)),Double_t(det(1)));  //vecchio
 937     Double_t phidet= TMath::ATan2(Double_t(det(1)),Double_t(det(0))); if(phidet<0.) phidet+=2.*TMath::Pi(); //nuovo
 938   // cout<<" layer phidet e det(6) = "<< layer<<" "<<phidet<<" "<<det(6)<<"\n"; getchar();
 939   Double_t xmin,ymin,xmax,ymax;
 940  // Double_t phiconfr=0.0;
 941   //cout<<" phiconfr inizio =  "<<phiconfr <<"\n"; getchar();
 942   local[1]=local[2]=0.;
 943   local[0]= -(fDetx[layer-1]);
 944   if(layer==1)    local[0]= (fDetx[layer-1]);  //take into account different reference system
 945   g2->LtoG(layer,iLd,detector,local,global);
 946   xmax=global[0]; ymax=global[1];
 947   local[0]= (fDetx[layer-1]);
 948   if(layer==1)    local[0]= -(fDetx[layer-1]);  //take into account different reference system
 949   g2->LtoG(layer,iLd,detector,local,global);
 950   xmin=global[0]; ymin=global[1];
 951  // Double_t phimin=PhiDef(xmin,ymin);  //vecchio
 952      Double_t phimin= TMath::ATan2(ymin,xmin); if(phimin<0.) phimin+=2.*TMath::Pi();  //nuovo
 953  // Double_t phimax=PhiDef(xmax,ymax);
 954     Double_t phimax= TMath::ATan2(ymax,xmax); if(phimax<0.) phimax+=2.*TMath::Pi();  //nuovo
 955   //cout<<" xmin ymin = "<<xmin<<" "<<ymin<<"\n";
 956   // cout<<" xmax ymax = "<<xmax<<" "<<ymax<<"\n";
 957   // cout<<" iLd phimin phimax ="<<iLd<<" "<<phimin<<" "<<phimax<<"\n";
 958
 959   Double_t phiconfr=phinters;
 960  if(phimin>phimax ){
 961      if(phimin <5.5) {cout<<" Error in Intersection for phi \n"; getchar();}
 962     phimin=phimin-(2.*pigre);
 963     if(phinters>(1.5*pigre)) phiconfr=phinters-(2.*pigre);
 964     if(phidet>(1.5*pigre)) phidet=phidet-(2.*pigre);
 965   }
 966   //  cout<<" phiconfr finale = "<<phiconfr<<"\n"; getchar();
 967   if(phiconfr>phimin && phiconfr<= phimax) {
 968     if(ip>1) {
 969       cout<< " Errore su ip in Intersection \n"; getchar();
 970     }
 971       else  {
 972         listLad(ip)= iLd; distPhiCenter(ip)=TMath::Abs(phiconfr-phidet); ip++;
 973       }
 974     }
 975   }
 976   if(ip==0) { cout<< " No detector along phi \n"; getchar();}
 977   ladder=Int_t (listLad(0));
 978   if(ip>1 && (distPhiCenter(0)>distPhiCenter(1)))   ladder=Int_t (listLad(1));
 979
 980   return 0;
 981 }
 982
 983
 984 void AliITSTrackerV1::KalmanFilter(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]){
 985 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
 986 // Kalman filter without vertex constraint
 987
 988
 989   ////////////////////////////// Evaluation of the measurement vector /////////////////////////////////////
 990
 991   Double_t m[2];
 992   Double_t rk,phik,zk;
 993   rk=cluster(0);   phik=cluster(1);  zk=cluster(2);
 994   m[0]=phik;    m[1]=zk;
 995
 996   ///////////////////////////////////// Evaluation of the error matrix V  ///////////////////////////////
 997
 998   Double_t v00=sigma[0];
 999   Double_t v11=sigma[1];
1000
1001   ///////////////////////////////////////////////////////////////////////////////////////////
1002
1003
1004   Double_t cin00,cin10,cin20,cin30,cin40,cin11,cin21,cin31,cin41,cin22,cin32,cin42,cin33,cin43,cin44;
1005
1006   newTrack->GetCElements(cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,
1007                          cin41,cin42,cin43,cin44); //get C matrix
1008
1009   Double_t rold00=cin00+v00;
1010   Double_t rold10=cin10;
1011   Double_t rold11=cin11+v11;
1012
1013 //////////////////////////////////// R matrix inversion  ///////////////////////////////////////////////
1014
1015   Double_t det=rold00*rold11-rold10*rold10;
1016   Double_t r00=rold11/det;
1017   Double_t r10=-rold10/det;
1018   Double_t r11=rold00/det;
1019
1020 ////////////////////////////////////////////////////////////////////////////////////////////////////////
1021
1022   Double_t k00=cin00*r00+cin10*r10;
1023   Double_t k01=cin00*r10+cin10*r11;
1024   Double_t k10=cin10*r00+cin11*r10;
1025   Double_t k11=cin10*r10+cin11*r11;
1026   Double_t k20=cin20*r00+cin21*r10;
1027   Double_t k21=cin20*r10+cin21*r11;
1028   Double_t k30=cin30*r00+cin31*r10;
1029   Double_t k31=cin30*r10+cin31*r11;
1030   Double_t k40=cin40*r00+cin41*r10;
1031   Double_t k41=cin40*r10+cin41*r11;
1032
1033   Double_t x0,x1,x2,x3,x4;
1034   newTrack->GetXElements(x0,x1,x2,x3,x4);     // get the state vector
1035
1036   Double_t savex0=x0, savex1=x1;
1037
1038   x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1);
1039   x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1);
1040   x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1);
1041   x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1);
1042   x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1);
1043
1044   Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44;
1045
1046   c00=cin00-k00*cin00-k01*cin10;
1047   c10=cin10-k00*cin10-k01*cin11;
1048   c20=cin20-k00*cin20-k01*cin21;
1049   c30=cin30-k00*cin30-k01*cin31;
1050   c40=cin40-k00*cin40-k01*cin41;
1051
1052   c11=cin11-k10*cin10-k11*cin11;
1053   c21=cin21-k10*cin20-k11*cin21;
1054   c31=cin31-k10*cin30-k11*cin31;
1055   c41=cin41-k10*cin40-k11*cin41;
1056
1057   c22=cin22-k20*cin20-k21*cin21;
1058   c32=cin32-k20*cin30-k21*cin31;
1059   c42=cin42-k20*cin40-k21*cin41;
1060
1061   c33=cin33-k30*cin30-k31*cin31;
1062   c43=cin43-k30*cin40-k31*cin41;
1063
1064   c44=cin44-k40*cin40-k41*cin41;
1065
1066   newTrack->PutXElements(x0,x1,x2,x3,x4);               // put the new state vector
1067
1068   newTrack->PutCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); // put in track the
1069                                                                                        // new cov matrix
1070   Double_t vmcold00=v00-c00;
1071   Double_t vmcold10=-c10;
1072   Double_t vmcold11=v11-c11;
1073
1074 ///////////////////////////////////// Matrix vmc inversion  ////////////////////////////////////////////////
1075
1076   det=vmcold00*vmcold11-vmcold10*vmcold10;
1077   Double_t vmc00=vmcold11/det;
1078   Double_t vmc10=-vmcold10/det;
1079   Double_t vmc11=vmcold00/det;
1080
1081 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
1082
1083   Double_t chi2=(m[0]-x0)*( vmc00*(m[0]-x0) + 2.*vmc10*(m[1]-x1) ) +
1084                 (m[1]-x1)*vmc11*(m[1]-x1);
1085
1086   newTrack->SetChi2(newTrack->GetChi2()+chi2);
1087
1088 }
1089
1090
1091 void AliITSTrackerV1::KalmanFilterVert(AliITSTrackV1 *newTrack,TVector &cluster,Double_t sigma[2]){
1092 //Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it
1093 // Kalman filter with vertex constraint
1094
1095   ////////////////////////////// Evaluation of the measurement vector m ///////////////
1096
1097   Double_t m[4];
1098   Double_t rk,phik,zk;
1099   rk=cluster(0);   phik=cluster(1);  zk=cluster(2);
1100   m[0]=phik;    m[1]=zk;
1101
1102   Double_t cc=(*newTrack).GetC();
1103   Double_t zv=(*newTrack).GetZv();
1104   Double_t dv=(*newTrack).GetDv();
1105   Double_t cy=cc/2.;
1106   Double_t tgl= (zk-zv)*cy/TMath::ASin(cy*rk);
1107   m[2]=dv;    m[3]=tgl;
1108
1109   ///////////////////////////////////// Evaluation of the error matrix V  //////////////
1110   Int_t layer=newTrack->GetLayer();
1111   Double_t v00=sigma[0];
1112   Double_t v11=sigma[1];
1113   Double_t v31=sigma[1]/rk;
1114   Double_t sigmaDv=newTrack->GetsigmaDv();
1115   Double_t v22=sigmaDv*sigmaDv  + newTrack->Getd2(layer-1);
1116   Double_t v32=newTrack->Getdtgl(layer-1);
1117   Double_t sigmaZv=newTrack->GetsigmaZv();
1118   Double_t v33=(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(layer-1);
1119   ///////////////////////////////////////////////////////////////////////////////////////
1120
1121   Double_t cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,cin41,cin42,cin43,cin44;
1122
1123   newTrack->GetCElements(cin00,cin10,cin11,cin20,cin21,cin22,cin30,cin31,cin32,cin33,cin40,
1124                          cin41,cin42,cin43,cin44); //get C matrix
1125
1126   Double_t r[4][4];
1127   r[0][0]=cin00+v00;
1128   r[1][0]=cin10;
1129   r[2][0]=cin20;
1130   r[3][0]=cin30;
1131   r[1][1]=cin11+v11;
1132   r[2][1]=cin21;
1133   r[3][1]=cin31+sigma[1]/rk;
1134   r[2][2]=cin22+sigmaDv*sigmaDv+newTrack->Getd2(layer-1);
1135   r[3][2]=cin32+newTrack->Getdtgl(layer-1);
1136   r[3][3]=cin33+(sigma[1]+sigmaZv*sigmaZv)/(rk*rk) + newTrack->Gettgl2(layer-1);
1137
1138   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];
1139   r[2][3]=r[3][2];
1140
1141 /////////////////////  Matrix R inversion ////////////////////////////////////////////
1142
1143   const Int_t kn=4;
1144   Double_t big, hold;
1145   Double_t d=1.;
1146   Int_t ll[kn],mm[kn];
1147
1148   Int_t i,j,k;
1149
1150   for(k=0; k<kn; k++) {
1151     ll[k]=k;
1152     mm[k]=k;
1153     big=r[k][k];
1154     for(j=k; j<kn ; j++) {
1155       for (i=j; i<kn; i++) {
1156         if(TMath::Abs(big) < TMath::Abs(r[i][j]) ) { big=r[i][j]; ll[k]=i; mm[k]=j; }
1157       }
1158     }
1159 //
1160     j= ll[k];
1161     if(j > k) {
1162       for(i=0; i<kn; i++) { hold=-r[k][i]; r[k][i]=r[j][i]; r[j][i]=hold; }
1163
1164     }
1165 //
1166     i=mm[k];
1167     if(i > k ) {
1168       for(j=0; j<kn; j++) { hold=-r[j][k]; r[j][k]=r[j][i]; r[j][i]=hold; }
1169     }
1170 //
1171     if(!big) {
1172       d=0.;
1173       cout << "Singular matrix\n";
1174     }
1175     for(i=0; i<kn; i++) {
1176       if(i == k) { continue; }
1177       r[i][k]=r[i][k]/(-big);
1178     }
1179 //
1180     for(i=0; i<kn; i++) {
1181       hold=r[i][k];
1182       for(j=0; j<kn; j++) {
1183         if(i == k || j == k) { continue; }
1184         r[i][j]=hold*r[k][j]+r[i][j];
1185       }
1186     }
1187 //
1188     for(j=0; j<kn; j++) {
1189       if(j == k) { continue; }
1190       r[k][j]=r[k][j]/big;
1191     }
1192 //
1193     d=d*big;
1194 //
1195     r[k][k]=1./big;
1196   }
1197 //
1198   for(k=kn-1; k>=0; k--) {
1199     i=ll[k];
1200     if(i > k) {
1201       for (j=0; j<kn; j++) {hold=r[j][k]; r[j][k]=-r[j][i]; r[j][i]=hold;}
1202     }
1203     j=mm[k];
1204     if(j > k) {
1205       for (i=0; i<kn; i++) {hold=r[k][i]; r[k][i]=-r[j][i]; r[j][i]=hold;}
1206       }
1207   }
1208 //////////////////////////////////////////////////////////////////////////////////
1209
1210
1211   Double_t k00=cin00*r[0][0]+cin10*r[1][0]+cin20*r[2][0]+cin30*r[3][0];
1212   Double_t k01=cin00*r[1][0]+cin10*r[1][1]+cin20*r[2][1]+cin30*r[3][1];
1213   Double_t k02=cin00*r[2][0]+cin10*r[2][1]+cin20*r[2][2]+cin30*r[3][2];
1214   Double_t k03=cin00*r[3][0]+cin10*r[3][1]+cin20*r[3][2]+cin30*r[3][3];
1215   Double_t k10=cin10*r[0][0]+cin11*r[1][0]+cin21*r[2][0]+cin31*r[3][0];
1216   Double_t k11=cin10*r[1][0]+cin11*r[1][1]+cin21*r[2][1]+cin31*r[3][1];
1217   Double_t k12=cin10*r[2][0]+cin11*r[2][1]+cin21*r[2][2]+cin31*r[3][2];
1218   Double_t k13=cin10*r[3][0]+cin11*r[3][1]+cin21*r[3][2]+cin31*r[3][3];
1219   Double_t k20=cin20*r[0][0]+cin21*r[1][0]+cin22*r[2][0]+cin32*r[3][0];
1220   Double_t k21=cin20*r[1][0]+cin21*r[1][1]+cin22*r[2][1]+cin32*r[3][1];
1221   Double_t k22=cin20*r[2][0]+cin21*r[2][1]+cin22*r[2][2]+cin32*r[3][2];
1222   Double_t k23=cin20*r[3][0]+cin21*r[3][1]+cin22*r[3][2]+cin32*r[3][3];
1223   Double_t k30=cin30*r[0][0]+cin31*r[1][0]+cin32*r[2][0]+cin33*r[3][0];
1224   Double_t k31=cin30*r[1][0]+cin31*r[1][1]+cin32*r[2][1]+cin33*r[3][1];
1225   Double_t k32=cin30*r[2][0]+cin31*r[2][1]+cin32*r[2][2]+cin33*r[3][2];
1226   Double_t k33=cin30*r[3][0]+cin31*r[3][1]+cin32*r[3][2]+cin33*r[3][3];
1227   Double_t k40=cin40*r[0][0]+cin41*r[1][0]+cin42*r[2][0]+cin43*r[3][0];
1228   Double_t k41=cin40*r[1][0]+cin41*r[1][1]+cin42*r[2][1]+cin43*r[3][1];
1229   Double_t k42=cin40*r[2][0]+cin41*r[2][1]+cin42*r[2][2]+cin43*r[3][2];
1230   Double_t k43=cin40*r[3][0]+cin41*r[3][1]+cin42*r[3][2]+cin43*r[3][3];
1231
1232   Double_t x0,x1,x2,x3,x4;
1233   newTrack->GetXElements(x0,x1,x2,x3,x4);     // get the state vector
1234
1235   Double_t savex0=x0, savex1=x1, savex2=x2, savex3=x3;
1236
1237   x0+=k00*(m[0]-savex0)+k01*(m[1]-savex1)+k02*(m[2]-savex2)+
1238       k03*(m[3]-savex3);
1239   x1+=k10*(m[0]-savex0)+k11*(m[1]-savex1)+k12*(m[2]-savex2)+
1240       k13*(m[3]-savex3);
1241   x2+=k20*(m[0]-savex0)+k21*(m[1]-savex1)+k22*(m[2]-savex2)+
1242       k23*(m[3]-savex3);
1243   x3+=k30*(m[0]-savex0)+k31*(m[1]-savex1)+k32*(m[2]-savex2)+
1244       k33*(m[3]-savex3);
1245   x4+=k40*(m[0]-savex0)+k41*(m[1]-savex1)+k42*(m[2]-savex2)+
1246       k43*(m[3]-savex3);
1247
1248   Double_t c00,c10,c20,c30,c40,c11,c21,c31,c41,c22,c32,c42,c33,c43,c44;
1249
1250   c00=cin00-k00*cin00-k01*cin10-k02*cin20-k03*cin30;
1251   c10=cin10-k00*cin10-k01*cin11-k02*cin21-k03*cin31;
1252   c20=cin20-k00*cin20-k01*cin21-k02*cin22-k03*cin32;
1253   c30=cin30-k00*cin30-k01*cin31-k02*cin32-k03*cin33;
1254   c40=cin40-k00*cin40-k01*cin41-k02*cin42-k03*cin43;
1255
1256   c11=cin11-k10*cin10-k11*cin11-k12*cin21-k13*cin31;
1257   c21=cin21-k10*cin20-k11*cin21-k12*cin22-k13*cin32;
1258   c31=cin31-k10*cin30-k11*cin31-k12*cin32-k13*cin33;
1259   c41=cin41-k10*cin40-k11*cin41-k12*cin42-k13*cin43;
1260
1261   c22=cin22-k20*cin20-k21*cin21-k22*cin22-k23*cin32;
1262   c32=cin32-k20*cin30-k21*cin31-k22*cin32-k23*cin33;
1263   c42=cin42-k20*cin40-k21*cin41-k22*cin42-k23*cin43;
1264
1265   c33=cin33-k30*cin30-k31*cin31-k32*cin32-k33*cin33;
1266   c43=cin43-k30*cin40-k31*cin41-k32*cin42-k33*cin43;
1267
1268   c44=cin44-k40*cin40-k41*cin41-k42*cin42-k43*cin43;
1269
1270   newTrack->PutXElements(x0,x1,x2,x3,x4);               // put the new state vector
1271
1272   newTrack->PutCElements(c00,c10,c11,c20,c21,c22,c30,c31,c32,c33,c40,c41,c42,c43,c44); // put in track the
1273                                                                                        // new cov matrix
1274
1275   Double_t vmc[4][4];
1276
1277   vmc[0][0]=v00-c00; vmc[1][0]=-c10; vmc[2][0]=-c20; vmc[3][0]=-c30;
1278   vmc[1][1]=v11-c11; vmc[2][1]=-c21; vmc[3][1]=v31-c31;
1279   vmc[2][2]=v22-c22; vmc[3][2]=v32-c32;
1280   vmc[3][3]=v33-c33;
1281
1282   vmc[0][1]=vmc[1][0]; vmc[0][2]=vmc[2][0]; vmc[0][3]=vmc[3][0];
1283   vmc[1][2]=vmc[2][1]; vmc[1][3]=vmc[3][1];
1284   vmc[2][3]=vmc[3][2];
1285
1286
1287 /////////////////////// vmc matrix inversion ///////////////////////////////////
1288
1289   d=1.;
1290
1291   for(k=0; k<kn; k++) {
1292     ll[k]=k;
1293     mm[k]=k;
1294     big=vmc[k][k];
1295     for(j=k; j<kn ; j++) {
1296       for (i=j; i<kn; i++) {
1297         if(TMath::Abs(big) < TMath::Abs(vmc[i][j]) ) { big=vmc[i][j]; ll[k]=i; mm[k]=j; }
1298       }
1299     }
1300 //
1301     j= ll[k];
1302     if(j > k) {
1303       for(i=0; i<kn; i++) { hold=-vmc[k][i]; vmc[k][i]=vmc[j][i]; vmc[j][i]=hold; }
1304
1305     }
1306 //
1307     i=mm[k];
1308     if(i > k ) {
1309       for(j=0; j<kn; j++) { hold=-vmc[j][k]; vmc[j][k]=vmc[j][i]; vmc[j][i]=hold; }
1310     }
1311 //
1312     if(!big) {
1313       d=0.;
1314       cout << "Singular matrix\n";
1315     }
1316     for(i=0; i<kn; i++) {
1317       if(i == k) { continue; }
1318       vmc[i][k]=vmc[i][k]/(-big);
1319     }
1320 //
1321     for(i=0; i<kn; i++) {
1322       hold=vmc[i][k];
1323       for(j=0; j<kn; j++) {
1324         if(i == k || j == k) { continue; }
1325         vmc[i][j]=hold*vmc[k][j]+vmc[i][j];
1326       }
1327     }
1328 //
1329     for(j=0; j<kn; j++) {
1330       if(j == k) { continue; }
1331       vmc[k][j]=vmc[k][j]/big;
1332     }
1333 //
1334     d=d*big;
1335 //
1336     vmc[k][k]=1./big;
1337   }
1338 //
1339   for(k=kn-1; k>=0; k--) {
1340     i=ll[k];
1341     if(i > k) {
1342       for (j=0; j<kn; j++) {hold=vmc[j][k]; vmc[j][k]=-vmc[j][i]; vmc[j][i]=hold;}
1343     }
1344     j=mm[k];
1345     if(j > k) {
1346       for (i=0; i<kn; i++) {hold=vmc[k][i]; vmc[k][i]=-vmc[j][i]; vmc[j][i]=hold;}
1347       }
1348   }
1349
1350
1351 ////////////////////////////////////////////////////////////////////////////////
1352
1353   Double_t chi2=(m[0]-x0)*( vmc[0][0]*(m[0]-x0) + 2.*vmc[1][0]*(m[1]-x1) +
1354                    2.*vmc[2][0]*(m[2]-x2)+ 2.*vmc[3][0]*(m[3]-x3) ) +
1355                 (m[1]-x1)* ( vmc[1][1]*(m[1]-x1) + 2.*vmc[2][1]*(m[2]-x2)+
1356                    2.*vmc[3][1]*(m[3]-x3) ) +
1357                 (m[2]-x2)* ( vmc[2][2]*(m[2]-x2)+ 2.*vmc[3][2]*(m[3]-x3) ) +
1358                 (m[3]-x3)*vmc[3][3]*(m[3]-x3);
1359
1360   newTrack->SetChi2(newTrack->GetChi2()+chi2);
1361
1362 }
1363