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