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