Changes made to be consitant with changes in TPC tracking classes (B.Nilsen)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 May 2001 06:04:58 +0000 (06:04 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 30 May 2001 06:04:58 +0000 (06:04 +0000)
ITS/AliITS.cxx
ITS/AliTPCtracknew.C [new file with mode: 0644]

index 2fddab3..f31ed6e 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.51  2001/05/16 14:57:15  alibrary
+New files for folders and Stack
+
 Revision 1.50  2001/05/11 09:15:21  barbera
 Corrected to make fast point creation working with PPR geometry
 
@@ -247,7 +250,7 @@ the AliITS class.
 #include "AliITSRad.h"   
 #include "../TPC/AliTPC.h"
 #include "../TPC/AliTPCParam.h"
-
+#include "../TPC/AliTPCtracker.h"
 
 ClassImp(AliITS)
  
@@ -717,9 +720,6 @@ void AliITS::SetDefaults()
   } 
 
 }
-
-
-
 //_____________________________________________________________________________
 void AliITS::SetDefaultSimulation()
 {
@@ -1290,8 +1290,6 @@ void AliITS::DigitsToRecPoints(Int_t evNumber,Int_t lastentry,Option_t *opt)
    treeC->Write(hname,TObject::kOverwrite);
    treeC->Reset();
 }
-
-
 //____________________________________________________________________________
 void AliITS::HitsToFastRecPoints(Int_t evNumber,Int_t bgrev,Int_t size,
 Option_t *option,Option_t *opt,Text_t *filename)
@@ -1307,7 +1305,7 @@ Option_t *option,Option_t *opt,Text_t *filename)
    //if(ver!=5) return; 
 
    const char *all = strstr(opt,"All");
-   const char *det[3] = {strstr(opt,"SPD"),strstr(opt,"SDD"),strstr(opt,"SSD")};
+   const char *det[3] ={strstr(opt,"SPD"),strstr(opt,"SDD"),strstr(opt,"SSD")};
 
    Int_t nmodules;
    InitModules(size,nmodules);
@@ -1335,7 +1333,8 @@ Option_t *option,Option_t *opt,Text_t *filename)
        AliITSDetType *iDetType=DetType(id); 
        sim = (AliITSsimulation*)iDetType->GetSimulationModel();
        if (!sim) {
-           Error("HitsToFastPoints","The simulation class was not instantiated!");
+           Error("HitsToFastPoints",
+                "The simulation class was not instantiated!");
            exit(1);
           // or SetDefaultSimulation();
        }
@@ -1368,14 +1367,13 @@ Option_t *option,Option_t *opt,Text_t *filename)
    delete [] random;
 
 }
-
 //________________________________________________________________
-AliITStrack  AliITS::Tracking(AliITStrack &track, AliITStrack *reference,TObjArray *fastpoints, Int_t
-**vettid, Bool_t flagvert,  AliITSRad *rl  ) { 
-                                                                               
-//Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it, Giuseppe.S.Pappalardo@ct.infn.it 
+AliITStrack  AliITS::Tracking(AliITStrack &track, AliITStrack *reference,
+                             TObjArray *fastpoints, Int_t **vettid,
+                             Bool_t flagvert,  AliITSRad *rl  ) {
+// Origin  A. Badala' and G.S. Pappalardo:  e-mail Angela.Badala@ct.infn.it,
+// Giuseppe.S.Pappalardo@ct.infn.it
 
-                                                                                   
   TList *list= new TList();   
 
   AliITStrack tr(track);
@@ -1385,7 +1383,8 @@ AliITStrack  AliITS::Tracking(AliITStrack &track, AliITStrack *reference,TObjArr
   Double_t Pt=(tr).GetPt();
 //  cout << "\n Pt = " << Pt <<"\n";  //stampa
 
-  AliITStracking obj(list, reference, this, fastpoints,TMath::Abs(Pt),vettid, flagvert, rl);
+  AliITStracking obj(list, reference, this, fastpoints,
+                    TMath::Abs(Pt),vettid, flagvert, rl);
   list->Delete();
   delete list;
 
@@ -1395,17 +1394,18 @@ AliITStrack  AliITS::Tracking(AliITStrack &track, AliITStrack *reference,TObjArr
   for(lay=5; lay>=0; lay--) {
     TVector VecLabref(3); 
     VecLabref=(*reference).GetLabTrack(lay);
-    Float_t ClustZ=(*reference).GetZclusterTrack( lay);   //aggiunta il 5-3-2001
-    for(k=0; k<3; k++){ //{itot++; VecTotLabref(itot)=VecLabref(k);} //cambiata il 5-3-2001 
-               Int_t lpp=(Int_t)VecLabref(k);
-               if(lpp>=0) {
-                 TParticle *p=(TParticle*) gAlice->Particle(lpp);
-                 Int_t pcode=p->GetPdgCode();
-                 if(pcode==11) VecLabref(k)=p->GetFirstMother();
-               }    
-    itot++; VecTotLabref(itot)=VecLabref(k);
-    if(VecLabref(k)==0. && ClustZ == 0.) VecTotLabref(itot) =-3.; }  
-  }
+    Float_t ClustZ=(*reference).GetZclusterTrack( lay); //aggiunta il 5-3-2001
+    for(k=0; k<3; k++){
+       //{itot++; VecTotLabref(itot)=VecLabref(k);} //cambiata il 5-3-2001 
+       Int_t lpp=(Int_t)VecLabref(k);
+       if(lpp>=0) {
+           TParticle *p=(TParticle*) gAlice->Particle(lpp);
+           Int_t pcode=p->GetPdgCode();
+           if(pcode==11) VecLabref(k)=p->GetFirstMother();
+       } // end if lpp>=0    
+       itot++; VecTotLabref(itot)=VecLabref(k);
+       if(VecLabref(k)==0. && ClustZ == 0.) VecTotLabref(itot) =-3.; }  
+  } // end for lay
   Long_t labref;
   Int_t freq;  
   (*reference).Search(VecTotLabref, labref, freq);
@@ -1418,145 +1418,144 @@ AliITStrack  AliITS::Tracking(AliITStrack &track, AliITStrack *reference,TObjArr
   return *reference; 
 
 }
+//________________________________________________________________------
+void AliITS::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t,
+                       TFile *file, Bool_t flagvert) {
+    //   ex macro for tracking ITS
 
+    printf("begin DoTracking - file %p\n",file);
 
+    //const char *pname="75x40_100x60";
+  
+    Int_t imax=200,jmax=450;
+    AliITSRad *rl = new AliITSRad(imax,jmax);
+    //cout<<" dopo costruttore AliITSRad\n"; getchar();
 
-//________________________________________________________________
+    struct GoodTrack {
+       Int_t lab,code;
+       Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
+       Bool_t flag;
+    };
 
+    gAlice->GetEvent(0);
 
+    AliKalmanTrack *kkprov;
+    kkprov->SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());  
 
-void AliITS::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, Bool_t flagvert) {
+    /*  //modificato il 26-4-2001
+       AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
+       AliTPCParam *digp = (AliTPCParam*)file->Get(pname);
+       if (digp!=0) TPC->SetParam(digp);
+    */
+    TFile *cf=TFile::Open("AliTPCclusters.root");  
+    AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
+    if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
+
+    AliTPCtracker *tracker = new AliTPCtracker(digp); //aggiunto il 23-5 
+    //aggiunto il 23-5 
+    // Load clusters
+    tracker->LoadInnerSectors();
+    tracker->LoadOuterSectors();
+
+    GoodTrack gt[15000];
+    Int_t ngood=0;
+    ifstream in("itsgood_tracks");
+
+    cerr<<"Reading itsgood tracks...\n";
+    while (in>>gt[ngood].lab>>gt[ngood].code
+          >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
+          >>gt[ngood].x  >>gt[ngood].y >>gt[ngood].z
+          >>gt[ngood].pxg  >>gt[ngood].pyg >>gt[ngood].pzg
+          >>gt[ngood].ptg >>gt[ngood].flag) {
+       ngood++;
+       cerr<<ngood<<'\r';
+       if (ngood==15000) {
+           cerr<<"Too many good tracks !\n";
+           break;
+       } // end if ngood==1500
+    } // end while
+    if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
+
+    // Load tracks
+    // TFile *tf=TFile::Open("tpctracks.root"); //commentato 26-4-2001 
+    TFile *tf=TFile::Open("AliTPCtracks.root"); 
+    if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
+    TObjArray tracks(200000);
+    //TTree *tracktree=(TTree*)tf->Get("TreeT");
+    TTree *tracktree=(TTree*)tf->Get("TPCf"); //aggiunto il 23-5
+    if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n";} 
+    TBranch *tbranch=tracktree->GetBranch("tracks");
+    Int_t nentr=(Int_t)tracktree->GetEntries();
+    Int_t kk;
+    /*  commentato il 26-4-2001
+       for (kk=0; kk<nentr; kk++) {
+       AliTPCtrack *iotrack=new AliTPCtrack;
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(kk);
+       tracks.AddLast(iotrack);
+       } 
+    */
+    AliTPCtrack *iotracktpc=0;    
+    for (kk=0; kk<nentr; kk++) {
+       iotracktpc=new AliTPCtrack; 
+       tbranch->SetAddress(&iotracktpc);
+       tracktree->GetEvent(kk);
+       tracker->CookLabel(iotracktpc,0.1); //aggiunto 23-5
+       tracks.AddLast(iotracktpc);    
+    } // end for kk
+    delete tracker;
+    tf->Close();
 
-//   ex macro for tracking ITS
+    Int_t nt = tracks.GetEntriesFast();
+    cerr<<"Number of found tracks "<<nt<<endl;
 
-  printf("begin DoTracking - file %p\n",file);
+    TVector DataOut(9);
+    Int_t kkk=0;
 
-  //const char *pname="75x40_100x60";
-  
- Int_t imax=200,jmax=450;
-  AliITSRad *rl = new AliITSRad(imax,jmax);
-  //cout<<" dopo costruttore AliITSRad\n"; getchar();
-    
-  struct GoodTrack {
-    Int_t lab,code;
-    Float_t px,py,pz,x,y,z,pxg,pyg,pzg,ptg;
-    Bool_t flag;
-  };
-  
+    Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
 
-  gAlice->GetEvent(0);
-  //  AliKalmanTrack *kkprov;
-  //  kkprov->SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());  
-
-  /*  //modificato il 26-4-2001
-  AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
-  AliTPCParam *digp = (AliTPCParam*)file->Get(pname);
-  if (digp!=0) TPC->SetParam(digp);
-  */
-  TFile *cf=TFile::Open("AliTPCclusters.root");  
-  AliTPCParam *digp= (AliTPCParam*)cf->Get("75x40_100x60");
-  if (!digp) { cerr<<"TPC parameters have not been found !\n"; getchar();}
-     
-  
-  GoodTrack gt[15000];
-  Int_t ngood=0;
-  ifstream in("itsgood_tracks");
-
-  cerr<<"Reading itsgood tracks...\n";
-  while (in>>gt[ngood].lab>>gt[ngood].code
-         >>gt[ngood].px >>gt[ngood].py>>gt[ngood].pz
-         >>gt[ngood].x  >>gt[ngood].y >>gt[ngood].z
-         >>gt[ngood].pxg  >>gt[ngood].pyg >>gt[ngood].pzg
-         >>gt[ngood].ptg >>gt[ngood].flag) {
-    ngood++;
-    cerr<<ngood<<'\r';
-    if (ngood==15000) {
-      cerr<<"Too many good tracks !\n";
-      break;
-    }
-  }
-  if (!in.eof()) cerr<<"Read error (itsgood_tracks) !\n";
-  
-  
-// Load tracks
- // TFile *tf=TFile::Open("tpctracks.root"); //commentato 26-4-2001 
-  TFile *tf=TFile::Open("AliTPCtracks.root"); 
-  if (!tf->IsOpen()) {cerr<<"Can't open AliTPCtracks.root !\n"; return ;}
-  TObjArray tracks(200000);
-  TTree *tracktree=(TTree*)tf->Get("TreeT");
-  TBranch *tbranch=tracktree->GetBranch("tracks");
-  Int_t nentr=(Int_t)tracktree->GetEntries();
-  Int_t kk;
-  /*  commentato il 26-4-2001
-  for (kk=0; kk<nentr; kk++) {
-    AliTPCtrack *iotrack=new AliTPCtrack;
-    tbranch->SetAddress(&iotrack);
-    tracktree->GetEvent(kk);
-    tracks.AddLast(iotrack);
-  } 
-  */
-   AliTPCtrack *iotracktpc=0;    
-  for (kk=0; kk<nentr; kk++) {
-    iotracktpc=new AliTPCtrack; 
-    tbranch->SetAddress(&iotracktpc);
-    tracktree->GetEvent(kk);
-    tracks.AddLast(iotracktpc);    
-  }      
-  tf->Close();
-
-
-  Int_t nt = tracks.GetEntriesFast();
-  cerr<<"Number of found tracks "<<nt<<endl;
-  
-  TVector DataOut(9);
-  Int_t kkk=0;
-  
-  Double_t ptg=0.,pxg=0.,pyg=0.,pzg=0.;
-
-  //////////////////////////////  good tracks definition in TPC  ////////////////////////////////
-      
-  ofstream out1 ("AliITSTrag.out");
-  Int_t i;
-  for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
-  out1.close();
-
-
-  TVector vec(5);
-  TTree *TR=gAlice->TreeR();
-  Int_t nent=(Int_t)TR->GetEntries();
-  TClonesArray  *recPoints = RecPoints();
-  Int_t numbpoints;
-  Int_t totalpoints=0;
-  Int_t *np = new Int_t[nent];
-  Int_t **vettid = new Int_t* [nent];
-  Int_t mod;
-  for (mod=0; mod<nent; mod++) {
-    vettid[mod]=0;
-    this->ResetRecPoints();
-    //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
-    gAlice->TreeR()->GetEvent(mod); //first entry in TreeR is empty
-    numbpoints = recPoints->GetEntries();
-    totalpoints+=numbpoints;
-    np[mod] = numbpoints;
-  //cout<<" mod = "<<mod<<"   numbpoints = "<<numbpoints<<"\n"; getchar();
-    vettid[mod] = new Int_t[numbpoints];
-    Int_t ii;
-    for (ii=0;ii<numbpoints; ii++) *(vettid[mod]+ii)=0;
-  }
-
-  AliTPCtrack *track=0;
+    //////////////  good tracks definition in TPC  ///////////////
 
-     
-  if(min_t < 0) {min_t = 0; max_t = nt-1;}   
+    ofstream out1 ("AliITSTrag.out");
+    Int_t i;
+    for (i=0; i<ngood; i++) out1 << gt[i].ptg << "\n";
+    out1.close();
+
+
+    TVector vec(5);
+    TTree *TR=gAlice->TreeR();
+    Int_t nent=(Int_t)TR->GetEntries();
+    TClonesArray  *recPoints = RecPoints();
+    Int_t numbpoints;
+    Int_t totalpoints=0;
+    Int_t *np = new Int_t[nent];
+    Int_t **vettid = new Int_t* [nent];
+    Int_t mod;
+    for (mod=0; mod<nent; mod++) {
+       vettid[mod]=0;
+       this->ResetRecPoints();
+       //gAlice->TreeR()->GetEvent(mod+1); //first entry in TreeR is empty
+       gAlice->TreeR()->GetEvent(mod); //first entry in TreeR is empty
+       numbpoints = recPoints->GetEntries();
+       totalpoints+=numbpoints;
+       np[mod] = numbpoints;
+       //cout<<" mod = "<<mod<<"   numbpoints = "<<numbpoints<<"\n";getchar();
+       vettid[mod] = new Int_t[numbpoints];
+       Int_t ii;
+       for (ii=0;ii<numbpoints; ii++) *(vettid[mod]+ii)=0;
+    } // end for mod
+
+    AliTPCtrack *track=0;
+    if(min_t < 0) {min_t = 0; max_t = nt-1;}   
 
 /*
-  ///////////////////////////////// Definition of vertex end its error ////////////////////////////
-  ////////////////////////// In the future it will be given by a method ///////////////////////////
+  //////////////////// Definition of vertex end its error ////////////////////
+  ////////////// In the future it will be given by a method //////////////////
   Double_t Vx=0.;
   Double_t Vy=0.;
   Double_t Vz=0.;
-  
+
   Float_t sigmavx=0.0050;      // 50  microns
   Float_t sigmavy=0.0050;      // 50  microns
   Float_t sigmavz=0.010;       // 100 microns
@@ -1565,275 +1564,274 @@ void AliITS::DoTracking(Int_t evNumber, Int_t min_t, Int_t max_t, TFile *file, B
   TVector vertex(3), ervertex(3)
   vertex(0)=Vx; vertex(1)=Vy; vertex(2)=Vz;
   ervertex(0)=sigmavx;  ervertex(1)=sigmavy;  ervertex(2)=sigmavz;
-  /////////////////////////////////////////////////////////////////////////////////////////////////
-*/      
-
-  TTree tracktree1("TreeT","Tree with ITS tracks");
-  AliITSiotrack *iotrack=0;
-  tracktree1.Branch("ITStracks","AliITSiotrack",&iotrack,32000,0);
+  //////////////////////////////////////////////////////////////////////////
+*/
 
-  ofstream out ("AliITSTra.out");
-  //ofstream outprova ("AliITSprova.out"); //commentato il 26-4-2001 
-  
-  Int_t j;       
-  for (j=min_t; j<=max_t; j++) {     
-    track=(AliTPCtrack*)tracks.UncheckedAt(j);
-    Int_t flaglab=0;
-    if (!track) continue;
-    ////// elimination of not good tracks ////////////  
-    Int_t ilab=TMath::Abs(track->GetLabel());
-    Int_t iii;
-    for (iii=0;iii<ngood;iii++) {
-        //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n"; getchar();
-      if (ilab==gt[iii].lab) { 
-       flaglab=1;
-       ptg=gt[iii].ptg; 
-       pxg=gt[iii].pxg;
-       pyg=gt[iii].pyg;
-       pzg=gt[iii].pzg;        
-       break;
-      }
-    }
-        //cout<<" j flaglab =  " <<j<<" "<<flaglab<<"\n";  getchar();
-    if (!flaglab) continue;  
-        //cout<<" j =  " <<j<<"\n";  getchar();
-  /*           
-    ////// old propagation to the end of TPC //////////////       
-    Double_t xk=76.;
-    track->PropagateTo(xk);
-    xk-=0.11;
-    track->PropagateTo(xk,42.7,2.27); //C
-    xk-=2.6;
-    track->PropagateTo(xk,36.2,1.98e-3); //C02
-    xk-=0.051;
-    track->PropagateTo(xk,42.7,2.27); //C 
-    /////////////////////////////////////////////////// 
-        */
-                
-        ////// new propagation to the end of TPC //////////////
-    Double_t xk=77.415;
-    track->PropagateTo(xk, 28.94, 1.204e-3);    //Ne
-        xk -=0.01;
-    track->PropagateTo(xk, 44.77, 1.71);        //Tedlar
-        xk -=0.04;
-    track->PropagateTo(xk, 44.86, 1.45);        //Kevlar
-        xk -=2.0;
-    track->PropagateTo(xk, 41.28, 0.029);       //Nomex         
-    xk-=16;
-    track->PropagateTo(xk,36.2,1.98e-3); //C02
-        xk -=0.01;
-    track->PropagateTo(xk, 24.01, 2.7);         //Al    
-        xk -=0.01;
-    track->PropagateTo(xk, 44.77, 1.71);        //Tedlar
-        xk -=0.04;
-    track->PropagateTo(xk, 44.86, 1.45);        //Kevlar
-        xk -=0.5;
-    track->PropagateTo(xk, 41.28, 0.029);       //Nomex                                                 
-                    
-       ///////////////////////////////////////////////////////////////                          
-  
-   ///////////////////////////////////////////////////////////////
-    AliITStrack trackITS(*track);
-    AliITStrack result(*track);
-    AliITStrack primarytrack(*track); 
+    TTree tracktree1("TreeT","Tree with ITS tracks");
+    AliITSiotrack *iotrack=0;
+    tracktree1.Branch("ITStracks","AliITSiotrack",&iotrack,32000,0);
+
+    ofstream out ("AliITSTra.out");
+    //ofstream outprova ("AliITSprova.out"); //commentato il 26-4-2001
+
+    Int_t j;       
+    for (j=min_t; j<=max_t; j++) {     
+       track=(AliTPCtrack*)tracks.UncheckedAt(j);
+       Int_t flaglab=0;
+       if (!track) continue;
+       ////// elimination of not good tracks ////////////       
+       Int_t ilab=TMath::Abs(track->GetLabel());
+       Int_t iii;
+       for (iii=0;iii<ngood;iii++) {
+           //cout<<" ilab, gt[iii].lab = "<<ilab<<" "<<gt[iii].lab<<"\n";getchar();
+           if (ilab==gt[iii].lab) { 
+               flaglab=1;
+               ptg=gt[iii].ptg; 
+               pxg=gt[iii].pxg;
+               pyg=gt[iii].pyg;
+               pzg=gt[iii].pzg;        
+               break;
+           } // end if ilab==
+       } // end for iii
+       //cout<<" j flaglab =  " <<j<<" "<<flaglab<<"\n";  getchar();
+       if (!flaglab) continue;  
+       //cout<<" j =  " <<j<<"\n";  getchar();
+       /*
+         ////// old propagation to the end of TPC //////////////       
+         Double_t xk=76.;
+         track->PropagateTo(xk);
+         xk-=0.11;
+         track->PropagateTo(xk,42.7,2.27); //C
+         xk-=2.6;
+         track->PropagateTo(xk,36.2,1.98e-3); //C02
+         xk-=0.051;
+         track->PropagateTo(xk,42.7,2.27); //C 
+         /////////////////////////////////////////////////// 
+       */ 
+       ////// new propagation to the end of TPC //////////////
+       Double_t xk=77.415;
+       track->PropagateTo(xk, 28.94, 1.204e-3);         //Ne
+       xk -=0.01;
+       track->PropagateTo(xk, 44.77, 1.71);     //Tedlar
+       xk -=0.04;
+       track->PropagateTo(xk, 44.86, 1.45);     //Kevlar
+       xk -=2.0;
+       track->PropagateTo(xk, 41.28, 0.029);    //Nomex         
+       xk-=16;
+       track->PropagateTo(xk,36.2,1.98e-3); //C02
+       xk -=0.01;
+       track->PropagateTo(xk, 24.01, 2.7);      //Al    
+       xk -=0.01;
+       track->PropagateTo(xk, 44.77, 1.71);     //Tedlar
+       xk -=0.04;
+       track->PropagateTo(xk, 44.86, 1.45);     //Kevlar
+       xk -=0.5;
+       track->PropagateTo(xk, 41.28, 0.029);    //Nomex
+       ///////////////////////////////////////////////////////////////
+
+       ///////////////////////////////////////////////////////////////
+       AliITStrack trackITS(*track);
+       AliITStrack result(*track);
+       AliITStrack primarytrack(*track); 
     
-///////////////////////////////////////////////////////////////////////////////////////////////
-        TVector Vgeant(3);
-        Vgeant=result.GetVertex(); 
-                         
-  // Definition of Dv and Zv for vertex constraint     
-     Double_t sigmaDv=0.0050;  Double_t sigmaZv=0.010; 
-    //Double_t sigmaDv=0.0015;  Double_t sigmaZv=0.0015;                                 
+       /////////////////////////////////////////////////////////////////
+       TVector Vgeant(3);
+       Vgeant=result.GetVertex();
+
+       // Definition of Dv and Zv for vertex constraint        
+       Double_t sigmaDv=0.0050;  Double_t sigmaZv=0.010;       
+       //Double_t sigmaDv=0.0015;  Double_t sigmaZv=0.0015;
        Double_t uniform= gRandom->Uniform();
        Double_t signdv;
        if(uniform<=0.5) signdv=-1.;
-          else
-                signdv=1.;
-        
-       Double_t Vr=TMath::Sqrt(Vgeant(0)*Vgeant(0)+ Vgeant(1)*Vgeant(1));
-         Double_t Dv=gRandom->Gaus(signdv*Vr,(Float_t)sigmaDv); 
-    Double_t Zv=gRandom->Gaus(Vgeant(2),(Float_t)sigmaZv);
-                               
-  //cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";                            
-    trackITS.SetDv(Dv);  trackITS.SetZv(Zv);
-    trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv); 
-    result.SetDv(Dv);  result.SetZv(Zv);
-    result.SetsigmaDv(sigmaDv); result.SetsigmaZv(sigmaZv);
-    primarytrack.SetDv(Dv);  primarytrack.SetZv(Zv);
-    primarytrack.SetsigmaDv(sigmaDv); primarytrack.SetsigmaZv(sigmaZv);                                                                
-
-/////////////////////////////////////////////////////////////////////////////////////////////////       
-               
-    primarytrack.PrimaryTrack(rl);
-    TVector  d2=primarytrack.Getd2();
-    TVector  tgl2=primarytrack.Gettgl2();
-    TVector  dtgl=primarytrack.Getdtgl();
-    trackITS.Setd2(d2); trackITS.Settgl2(tgl2);  trackITS.Setdtgl(dtgl); 
-    result.Setd2(d2); result.Settgl2(tgl2);  result.Setdtgl(dtgl);        
-        /*                      
-    trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
-    result.SetVertex(vertex);   result.SetErrorVertex(ervertex);   
-    */
-                         
-    Tracking(trackITS,&result,recPoints,vettid, flagvert,rl);  
-            
-    // cout<<" progressive track number = "<<j<<"\r";
-   // cout<<j<<"\r";
-    Int_t NumofCluster=result.GetNumClust();  
-//    cout<<" progressive track number = "<<j<<"\n";    // stampa
-    Long_t labITS=result.GetLabel();
-//    cout << " ITS track label = " << labITS << "\n";         // stampa           
-    int lab=track->GetLabel();             
-//    cout << " TPC track label = " << lab <<"\n";      // stampa
-        
-            
-//propagation to vertex
-       
-    Double_t rbeam=3.;
-     
-    result.Propagation(rbeam);
-                
-        Double_t C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44;
-        result.GetCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
-                
-    Double_t pt=TMath::Abs(result.GetPt());
-    Double_t Dr=result.GetD();
-    Double_t Z=result.GetZ();
-    Double_t tgl=result.GetTgl();
-    Double_t C=result.GetC();
-    Double_t Cy=C/2.;
-    Double_t Dz=Z-(tgl/Cy)*TMath::ASin(result.arga(rbeam));
-        Dz-=Vgeant(2);
-         
-        // cout<<" Dr e dz alla fine = "<<Dr<<" "<<Dz<<"\n"; getchar();
-    Double_t phi=result.Getphi();
-    Double_t phivertex = phi - TMath::ASin(result.argA(rbeam));
-    Double_t duepi=2.*TMath::Pi();      
-    if(phivertex>duepi) phivertex-=duepi;
-    if(phivertex<0.) phivertex+=duepi;
-    Double_t Dtot=TMath::Sqrt(Dr*Dr+Dz*Dz);
-        
-//////////////////////////////////////////////////////////////////////////////////////////     
-  
-    Int_t idmodule,idpoint;
-        if(NumofCluster >=5)  {            // cinque - sei
-        //if(NumofCluster ==6)  {            // cinque - sei 
-
-
-      AliITSiotrack outtrack;
-
-      iotrack=&outtrack;
-
-      iotrack->SetStatePhi(phi);
-      iotrack->SetStateZ(Z);
-      iotrack->SetStateD(Dr);
-      iotrack->SetStateTgl(tgl);
-      iotrack->SetStateC(C);
-               Double_t radius=result.Getrtrack();
-               iotrack->SetRadius(radius);
-               Int_t charge;
-               if(C>0.) charge=-1;  else charge=1;
-               iotrack->SetCharge(charge);
-               
-
+       else
+           signdv=1.;
 
-      iotrack->SetCovMatrix(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);  
-
-      Double_t px=pt*TMath::Cos(phivertex);
-      Double_t py=pt*TMath::Sin(phivertex);
-      Double_t pz=pt*tgl;
-               
-      Double_t xtrack=Dr*TMath::Sin(phivertex);
-      Double_t ytrack=Dr*TMath::Cos(phivertex);
-      Double_t ztrack=Dz+Vgeant(2);
-
-
-      iotrack->SetPx(px);
-      iotrack->SetPy(py);
-      iotrack->SetPz(pz);
-      iotrack->SetX(xtrack);
-      iotrack->SetY(ytrack);
-      iotrack->SetZ(ztrack);
-      iotrack->SetLabel(labITS);
-               
-      Int_t il;                
-               for(il=0;il<6; il++){
-                 iotrack->SetIdPoint(il,result.GetIdPoint(il));
-                 iotrack->SetIdModule(il,result.GetIdModule(il));              
-               }
-      tracktree1.Fill();
-
-   //cout<<" labITS = "<<labITS<<"\n";
-       //cout<<" phi z Dr tgl C = "<<phi<<" "<<Z<<" "<<Dr<<" "<<tgl<<" "<<C<<"\n";  getchar();    
-
-     DataOut(kkk) = ptg; kkk++; DataOut(kkk)=labITS; kkk++; DataOut(kkk)=lab; kkk++;           
-
-      for (il=0;il<6;il++) {
-        idpoint=result.GetIdPoint(il);
-        idmodule=result.GetIdModule(il);
-       *(vettid[idmodule]+idpoint)=1; 
-       iotrack->SetIdPoint(il,idpoint);
-        iotrack->SetIdModule(il,idmodule);
-      }
-      
-    //  cout<<"  +++++++++++++  pt e ptg = "<<pt<<" "<<ptg<<"  ++++++++++\n";
-      Double_t difpt= (pt-ptg)/ptg*100.;                                       
-      DataOut(kkk)=difpt; kkk++;                                             
-      Double_t lambdag=TMath::ATan(pzg/ptg);
-      Double_t   lam=TMath::ATan(tgl);      
-      Double_t diflam = (lam - lambdag)*1000.;
-      DataOut(kkk) = diflam; kkk++;                                        
-      Double_t phig=TMath::ATan2(pyg,pxg);  if(phig<0) phig=2.*TMath::Pi()+phig;       
-      Double_t phi=phivertex;
-        
-      Double_t difphi = (phi - phig)*1000.;
-      DataOut(kkk)=difphi; kkk++;
-      DataOut(kkk)=Dtot*1.e4; kkk++;
-      DataOut(kkk)=Dr*1.e4; kkk++;
-      DataOut(kkk)=Dz*1.e4; kkk++; 
-      Int_t r;
-      for (r=0; r<9; r++) { out<<DataOut(r)<<" ";}
-      out<<"\n";
-      kkk=0;  
+       Double_t Vr=TMath::Sqrt(Vgeant(0)*Vgeant(0)+ Vgeant(1)*Vgeant(1));
+       Double_t Dv=gRandom->Gaus(signdv*Vr,(Float_t)sigmaDv); 
+       Double_t Zv=gRandom->Gaus(Vgeant(2),(Float_t)sigmaZv);
                
-           
-    } // end if on NumofCluster
-  //gObjectTable->Print();    // stampa memoria     
-  }  //  end for (int j=min_t; j<=max_t; j++)
-  
-  out.close();  
-  //outprova.close();   
+       //cout<<" Dv e Zv = "<<Dv<<" "<<Zv<<"\n";
+       trackITS.SetDv(Dv);  trackITS.SetZv(Zv);
+       trackITS.SetsigmaDv(sigmaDv); trackITS.SetsigmaZv(sigmaZv); 
+       result.SetDv(Dv);  result.SetZv(Zv);
+       result.SetsigmaDv(sigmaDv); result.SetsigmaZv(sigmaZv);
+       primarytrack.SetDv(Dv);  primarytrack.SetZv(Zv);
+       primarytrack.SetsigmaDv(sigmaDv); primarytrack.SetsigmaZv(sigmaZv);
+
+       ////////////////////////////////////////////////////////////////
+
+       primarytrack.PrimaryTrack(rl);
+       TVector  d2=primarytrack.Getd2();
+       TVector  tgl2=primarytrack.Gettgl2();
+       TVector  dtgl=primarytrack.Getdtgl();
+       trackITS.Setd2(d2); trackITS.Settgl2(tgl2);  trackITS.Setdtgl(dtgl); 
+       result.Setd2(d2); result.Settgl2(tgl2);  result.Setdtgl(dtgl);     
+       /*                       
+          trackITS.SetVertex(vertex); trackITS.SetErrorVertex(ervertex);
+          result.SetVertex(vertex);   result.SetErrorVertex(ervertex);   
+       */
+
+
+       Tracking(trackITS,&result,recPoints,vettid, flagvert,rl);
+       // cout<<" progressive track number = "<<j<<"\r";
+       // cout<<j<<"\r";
+       Int_t NumofCluster=result.GetNumClust();  
+       // cout<<" progressive track number = "<<j<<"\n";    // stampa
+       Long_t labITS=result.GetLabel();
+       //  cout << " ITS track label = " << labITS << "\n";    // stampa
+       int lab=track->GetLabel();                  
+       // cout << " TPC track label = " << lab <<"\n";      // stampa
+
+       //propagation to vertex
+
+       Double_t rbeam=3.;
+
+       result.Propagation(rbeam);
+
+       Double_t C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44;
+       result.GetCElements(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);
+       Double_t pt=TMath::Abs(result.GetPt());
+       Double_t Dr=result.GetD();
+       Double_t Z=result.GetZ();
+       Double_t tgl=result.GetTgl();
+       Double_t C=result.GetC();
+       Double_t Cy=C/2.;
+       Double_t Dz=Z-(tgl/Cy)*TMath::ASin(result.arga(rbeam));
+       Dz-=Vgeant(2);
+
+       // cout<<" Dr e dz alla fine = "<<Dr<<" "<<Dz<<"\n"; getchar();
+       Double_t phi=result.Getphi();
+       Double_t phivertex = phi - TMath::ASin(result.argA(rbeam));
+       Double_t duepi=2.*TMath::Pi();   
+       if(phivertex>duepi) phivertex-=duepi;
+       if(phivertex<0.) phivertex+=duepi;
+       Double_t Dtot=TMath::Sqrt(Dr*Dr+Dz*Dz);
+
+       /////////////////////////////////////////////////////////////////////
+
+       Int_t idmodule,idpoint;
+       if(NumofCluster >=5)  {            // cinque - sei
+           //if(NumofCluster ==6)  {            // cinque - sei
+
+           AliITSiotrack outtrack;
+
+           iotrack=&outtrack;
+
+           iotrack->SetStatePhi(phi);
+           iotrack->SetStateZ(Z);
+           iotrack->SetStateD(Dr);
+           iotrack->SetStateTgl(tgl);
+           iotrack->SetStateC(C);
+           Double_t radius=result.Getrtrack();
+           iotrack->SetRadius(radius);
+           Int_t charge;
+           if(C>0.) charge=-1;  else charge=1;
+           iotrack->SetCharge(charge);
+
+
+           iotrack->SetCovMatrix(C00,C10,C11,C20,C21,C22,C30,C31,C32,C33,C40,C41,C42,C43,C44);  
+
+           Double_t px=pt*TMath::Cos(phivertex);
+           Double_t py=pt*TMath::Sin(phivertex);
+           Double_t pz=pt*tgl;
+
+           Double_t xtrack=Dr*TMath::Sin(phivertex);
+           Double_t ytrack=Dr*TMath::Cos(phivertex);
+           Double_t ztrack=Dz+Vgeant(2);
+
+
+           iotrack->SetPx(px);
+           iotrack->SetPy(py);
+           iotrack->SetPz(pz);
+           iotrack->SetX(xtrack);
+           iotrack->SetY(ytrack);
+           iotrack->SetZ(ztrack);
+           iotrack->SetLabel(labITS);
+
+           Int_t il;           
+           for(il=0;il<6; il++){
+               iotrack->SetIdPoint(il,result.GetIdPoint(il));
+               iotrack->SetIdModule(il,result.GetIdModule(il));
+           } // end for il
+           tracktree1.Fill();
+
+           //cout<<" labITS = "<<labITS<<"\n";
+           //cout<<" phi z Dr tgl C = "<<phi<<" "<<Z<<" "<<Dr<<" "<<tgl<<" "<<C<<"\n";  getchar();        
+
+           DataOut(kkk) = ptg; kkk++; DataOut(kkk)=labITS; kkk++; DataOut(kkk)=lab; kkk++;             
+
+           for (il=0;il<6;il++) {
+               idpoint=result.GetIdPoint(il);
+               idmodule=result.GetIdModule(il);
+               *(vettid[idmodule]+idpoint)=1; 
+               iotrack->SetIdPoint(il,idpoint);
+               iotrack->SetIdModule(il,idmodule);
+           } // end for il
+
+           //  cout<<"  +++++++++++++  pt e ptg = "<<pt<<" "<<ptg<<"  ++++++++++\n";
+           /* ///  provvisorio il 23-5-2001
+              Double_t pg=TMath::Sqrt(pxg*pxg+pyg*pyg+pzg*pzg);
+              Double_t cosl=TMath::Sqrt(1./(1.+tgl*tgl));
+              Double_t ptot=pt/cosl;
+              cout<<"ptot e pg = "<<ptot<<" "<<pg<<"\n";
+              Double_t difpt= (ptot-pg)/pg*100.;
+           */
+           ///////////////////////////////
+           Double_t difpt= (pt-ptg)/ptg*100.;  //coomentato prov il 23-5-2001
+           DataOut(kkk)=difpt; kkk++;
+           Double_t lambdag=TMath::ATan(pzg/ptg);
+           Double_t   lam=TMath::ATan(tgl);      
+           Double_t diflam = (lam - lambdag)*1000.;
+           DataOut(kkk) = diflam; kkk++;
+           Double_t phig=TMath::ATan2(pyg,pxg);
+           if(phig<0) phig=2.*TMath::Pi()+phig;       
+           Double_t phi=phivertex;
+
+           Double_t difphi = (phi - phig)*1000.;
+           DataOut(kkk)=difphi; kkk++;
+           DataOut(kkk)=Dtot*1.e4; kkk++;
+           DataOut(kkk)=Dr*1.e4; kkk++;
+           DataOut(kkk)=Dz*1.e4; kkk++; 
+           Int_t r;
+           for (r=0; r<9; r++) { out<<DataOut(r)<<" ";}
+           out<<"\n";
+           kkk=0;
+
+       } // end if on NumofCluster
+       //gObjectTable->Print();    // stampa memoria     
+    }  //  end for (int j=min_t; j<=max_t; j++)
+
+    out.close();  
+    //outprova.close();   
  
-  static Bool_t first=kTRUE;
-  static TFile *tfile;
+    static Bool_t first=kTRUE;
+    static TFile *tfile;
 
-       if(first) {
-           tfile=new TFile("itstracks.root","RECREATE");
-           //cout<<"I have opened itstracks.root file "<<endl;
-       }           
-       first=kFALSE;
-       tfile->cd();
-       tfile->ls();
+    if(first) {
+       tfile=new TFile("itstracks.root","RECREATE");
+       //cout<<"I have opened itstracks.root file "<<endl;
+    } // end if first      
+    first=kFALSE;
+    tfile->cd();
+    tfile->ls();
 
-   char hname[30];
-   sprintf(hname,"TreeT%d",evNumber);
+    char hname[30];
+    sprintf(hname,"TreeT%d",evNumber);
 
-  tracktree1.Write(hname);
+    tracktree1.Write(hname);
 
 
   
-           TTree *fAli=gAlice->TreeK();
-            TFile *fileAli=0;
+    TTree *fAli=gAlice->TreeK();
+    TFile *fileAli=0;
            
-           if (fAli) fileAli =fAli->GetCurrentFile();
-           fileAli->cd();
+    if (fAli) fileAli =fAli->GetCurrentFile();
+    fileAli->cd();
      
-  ////////////////////////////////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////
 
-  printf("delete vectors\n");
-  if(np) delete [] np;
-  if(vettid) delete [] vettid;
-  
+    printf("delete vectors\n");
+    if(np) delete [] np;
+    if(vettid) delete [] vettid;
 }
-
diff --git a/ITS/AliTPCtracknew.C b/ITS/AliTPCtracknew.C
new file mode 100644 (file)
index 0000000..b63909a
--- /dev/null
@@ -0,0 +1,370 @@
+/****************************************************************************
+ * This macro is supposed to do reconstruction in the barrel ALICE trackers *
+ * (Make sure you have TPC digits and ITS hits before using this macro !!!) *
+ * It does the following steps (April 12, 2001):                            *
+ *                   1) TPC cluster finding                                 *
+ *                   2) TPC track finding                                   *
+ *                   3) ITS cluster finding V2 (via fast points !)          *
+ *                   4) ITS track finding V2                                *
+ *                   5) ITS back track propagation V2                       *
+ *                   6) TPC back track propagation                          *
+ *                (Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch)          *
+ ****************************************************************************/
+
+#ifndef __CINT__
+  #include "alles.h"
+  #include "AliMagF.h"
+  #include "AliTPCtracker.h"
+
+  #include "AliITS.h"
+  #include "AliITSgeom.h"
+  #include "AliITSRecPoint.h"
+  #include "AliITSclusterV2.h"
+  #include "AliITSsimulationFastPoints.h"
+  #include "AliITStrackerV2.h"
+#endif
+
+Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname);
+Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname);
+Int_t TPCSortTracks(const Char_t *inname, const Char_t *outname);
+Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname);
+
+Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname);
+Int_t ITSFindTracks(const Char_t *inname, const Char_t *outname);
+Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname);
+
+
+Int_t AliTPCtracknew() {
+   const Char_t *TPCdigName="galice.root";
+   const Char_t *TPCclsName="AliTPCclusters.root";
+   const Char_t *TPCtrkName="AliTPCtracks.root";
+
+   const Char_t *ITSdigName="galice.root";
+   const Char_t *ITSclsName="AliITSclustersV2.root";
+   const Char_t *ITStrkName="AliITStracksV2.root";
+
+   AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
+  // AliKalmanTrack::SetConvConst(100/0.299792458/0.2/2);
+// ********** Find TPC clusters *********** //
+   if (TPCFindClusters(TPCdigName,TPCclsName)) {
+      cerr<<"Failed to get TPC clusters !\n";
+      return 1;
+   }      
+
+// ********** Find TPC tracks *********** //
+   if (TPCFindTracks(TPCclsName,TPCtrkName)) {
+      cerr<<"Failed to get TPC tracks !\n";
+      return 1;
+   }
+
+// ********** Sort and label TPC tracks *********** //
+   if (TPCSortTracks(TPCclsName,TPCtrkName)) {
+      cerr<<"Failed to sort TPC tracks !\n";
+      return 1;
+   }
+/*
+// ********** Find ITS clusters *********** //
+   if (ITSFindClusters(ITSdigName,ITSclsName)) {
+      cerr<<"Failed to get ITS clusters !\n";
+      return 1;
+   }
+
+// ********** Find ITS tracks *********** //
+   {TFile *clsFile=TFile::Open(ITSclsName);
+   if (ITSFindTracks(TPCtrkName,ITStrkName)) {
+      cerr<<"Failed to get ITS tracks !\n";
+      return 1;
+   }
+   clsFile->Close();}
+
+// ********** Back propagation of the ITS tracks *********** //
+   {TFile *clsFile=TFile::Open(ITSclsName);
+   if (ITSPropagateBack(ITStrkName,TPCtrkName)) {
+      cerr<<"Failed to propagate back the ITS tracks !\n";
+      return 1;
+   }
+   clsFile->Close();}
+
+
+// ********** Back propagation of the TPC tracks *********** //
+   {TFile *clsFile=TFile::Open(TPCclsName);
+   if (TPCPropagateBack(TPCtrkName,TPCtrkName)) {
+      cerr<<"Failed to propagate back the TPC tracks !\n";
+      return 1;
+   }
+   clsFile->Close();}
+*/
+   return 0;
+}
+
+
+Int_t TPCFindClusters(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="TPCFindClusters";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname);
+
+   AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+   AliTPCv2 tpc;
+   tpc.SetParam(param);
+
+   tpc.Digits2Clusters(out);
+
+   in->Close();
+   out->Close();
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t TPCFindTracks(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="TPCFindTracks";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname);
+   AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+
+   AliTPCtracker *tracker=new AliTPCtracker(param);
+   rc=tracker->Clusters2Tracks(0,out);
+   delete tracker;
+   in->Close();
+   out->Close();
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t TPCSortTracks(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="TPCSortTracks";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   TFile *out=TFile::Open(outname);
+   TFile *in =TFile::Open(inname);
+   AliTPCParam *param=(AliTPCParam *)in->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+
+   TObjArray tarray(10000);
+   AliTPCtrack *iotrack=0;
+   Int_t i;
+
+   out->cd();
+   TTree *tracktree=(TTree*)out->Get("TPCf");
+   TBranch *tbranch=tracktree->GetBranch("tracks");
+   Int_t nentr=(Int_t)tracktree->GetEntries();
+   for (i=0; i<nentr; i++) {
+       iotrack=new AliTPCtrack;
+       tbranch->SetAddress(&iotrack);
+       tracktree->GetEvent(i);
+       tarray.AddLast(iotrack);
+   }   
+   tarray.Sort();
+   out->Close();
+   
+   //assign thacks GEANT labels
+   in->cd();
+   AliTPCtracker *tracker = new AliTPCtracker(param);
+   tracker->LoadInnerSectors();
+   tracker->LoadOuterSectors();
+   for (i=0; i<nentr; i++) {
+       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
+       tracker->CookLabel(iotrack,0.1);
+   }   
+   delete tracker;
+   in->Close();
+   //end of GEANT label assignment
+
+   out=TFile::Open(outname,"recreate");
+   tracktree=new TTree("TPCf","Tree with TPC tracks");
+   tracktree->Branch("tracks","AliTPCtrack",&iotrack,32000,0);
+   for (i=0; i<nentr; i++) {
+       iotrack=(AliTPCtrack*)tarray.UncheckedAt(i);
+       tracktree->Fill();
+   }
+   tracktree->Write();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t ITSFindClusters(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="ITSFindClusters";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname,"update");
+
+   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+      cerr<<"Can't get gAlice !\n";
+      return 1;
+   }
+   Int_t ev=0;
+   gAlice->GetEvent(ev);
+   AliITS *ITS  = (AliITS*)gAlice->GetModule("ITS");
+   if (!ITS) { cerr<<"Can't get the ITS !\n"; return 1;}
+
+   gAlice->MakeTree("R"); ITS->MakeBranch("R",0);
+   //////////////// Taken from ITSHitsToFastPoints.C ///////////////////////
+   AliITSsimulationFastPoints *sim = new AliITSsimulationFastPoints();
+   for (Int_t i=0;i<3;i++) { ITS->SetSimulationModel(i,sim); }
+   Int_t nsignal=25;
+   Int_t size=-1;
+   Int_t bgr_ev=Int_t(ev/nsignal);
+   ITS->HitsToFastRecPoints(ev,bgr_ev,size," ","All"," ");
+   //////////////////////////////////////////////////////////////////////////
+
+   gAlice->GetEvent(ev);
+
+   AliITSgeom *geom=ITS->GetITSgeom();
+   out->cd();
+   geom->Write();
+
+   TClonesArray *clusters=new TClonesArray("AliITSclusterV2",10000);
+   TTree *cTree=new TTree("cTree","ITS clusters");
+   cTree->Branch("Clusters",&clusters);
+
+   TTree *pTree=gAlice->TreeR();
+   if (!pTree) { cerr<<"Can't get TreeR !\n"; return 1; }
+   TBranch *branch=pTree->GetBranch("ITSRecPoints");
+   if (!branch) { cerr<<"Can't get ITSRecPoints branch !\n"; return 1;}
+   TClonesArray *points=new TClonesArray("AliITSRecPoint",10000);
+   branch->SetAddress(&points);
+
+   TClonesArray &cl=*clusters;
+   Int_t nclusters=0;
+   Int_t nentr=(Int_t)pTree->GetEntries();
+   for (Int_t i=0; i<nentr; i++) {
+       if (!pTree->GetEvent(i)) {cTree->Fill(); continue;}
+       Int_t lay,lad,det; geom->GetModuleId(i,lay,lad,det);
+       Float_t x,y,zshift; geom->GetTrans(lay,lad,det,x,y,zshift); 
+       Double_t rot[9];    geom->GetRotMatrix(lay,lad,det,rot);
+       Double_t yshift = x*rot[0] + y*rot[1];
+       Int_t ndet=(lad-1)*geom->GetNdetectors(lay) + (det-1);
+       Int_t ncl=points->GetEntriesFast();
+       nclusters+=ncl;
+       for (Int_t j=0; j<ncl; j++) {
+          AliITSRecPoint *p=(AliITSRecPoint*)points->UncheckedAt(j);
+          Float_t lp[5];
+          lp[0]=-p->GetX()-yshift; if (lay==1) lp[0]=-lp[0];
+          lp[1]=p->GetZ()+zshift;
+          lp[2]=p->GetSigmaX2();
+          lp[3]=p->GetSigmaZ2();
+          lp[4]=p->GetQ();
+          Int_t lab[6]; 
+          lab[0]=p->GetLabel(0);lab[1]=p->GetLabel(1);lab[2]=p->GetLabel(2);
+          lab[3]=ndet;
+
+          Int_t label=lab[0];
+          TParticle *part=(TParticle*)gAlice->Particle(label);
+          label=-3;
+          while (part->P() < 0.005) {
+             Int_t m=part->GetFirstMother();
+             if (m<0) {cerr<<"Primary momentum: "<<part->P()<<endl; break;}
+             label=m;
+             part=(TParticle*)gAlice->Particle(label);
+          }
+          if      (lab[1]<0) lab[1]=label;
+          else if (lab[2]<0) lab[2]=label;
+          else cerr<<"No empty labels !\n";
+
+          new(cl[j]) AliITSclusterV2(lab,lp);
+       }
+       cTree->Fill(); clusters->Delete();
+       points->Delete();
+   }
+   cTree->Write();
+   cerr<<"Number of clusters: "<<nclusters<<endl;
+
+   delete cTree; delete clusters; delete points;
+
+   delete gAlice; gAlice=0;
+   in->Close();
+   out->Close();
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t ITSFindTracks(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="ITSFindTracks";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+   AliITStrackerV2 tracker(geom);
+
+   TFile *out=TFile::Open(outname,"recreate");
+   TFile *in =TFile::Open(inname);
+   rc=tracker.Clusters2Tracks(in,out);
+   in->Close();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t ITSPropagateBack(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="ITSPropagateBack";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   AliITSgeom *geom=(AliITSgeom*)gFile->Get("AliITSgeom");
+   if (!geom) { cerr<<"can't get ITS geometry !\n"; return 1;}
+   AliITStrackerV2 tracker(geom);
+
+   TFile *out=TFile::Open(outname,"update");
+   TFile *in =TFile::Open(inname);
+   rc=tracker.PropagateBack(in,out);
+   in->Close();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+
+Int_t TPCPropagateBack(const Char_t *inname, const Char_t *outname) {
+   Int_t rc=0;
+   const Char_t *name="TPCPropagateBack";
+   cerr<<'\n'<<name<<"...\n";
+   gBenchmark->Start(name);
+
+   AliTPCParam *param=(AliTPCParam *)gFile->Get("75x40_100x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 1;}
+   AliTPCtracker *tracker=new AliTPCtracker(param);
+
+   TFile *out=TFile::Open(outname,"update");
+   TFile *in =TFile::Open(inname);
+   rc=tracker->PropagateBack(in,out);
+   delete tracker;
+   in->Close();
+   out->Close();
+
+   gBenchmark->Stop(name);
+   gBenchmark->Show(name);
+
+   return rc;
+}
+