Comparison working with ESD (Yu.Belikov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Nov 2003 07:58:43 +0000 (07:58 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 21 Nov 2003 07:58:43 +0000 (07:58 +0000)
ITS/AliCascadeComparison.C
ITS/AliITSComparisonV2.C
ITS/AliV0Comparison.C
TPC/AliTPCComparison.C

index 0d1c04b..d21aa88 100644 (file)
   #include <fstream.h>
 
   #include "AliRun.h"
+  #include "AliMC.h"
   #include "AliHeader.h"
   #include "AliRunLoader.h"
   #include "AliITSLoader.h"
 
   #include "TH1.h"
   #include "TFile.h"
-  #include "TTree.h"
+  #include "TKey.h"
   #include "TObjArray.h"
   #include "TStyle.h"
   #include "TCanvas.h"
@@ -29,7 +30,7 @@
   #include "TPDGCode.h"
 
   #include "AliRun.h"
-  #include "AliPDG.h"
+  #include "AliESD.h"
   #include "AliCascadeVertex.h"
 #endif
 
@@ -49,11 +50,57 @@ Int_t AliCascadeComparison(Int_t code=3312) {
   //code=-3312; //kXiPlusBar
   //code= 3334; //kOmegaMinus
   //code=-3334; //kOmegaPlusBar
-
    cerr<<"Doing comparison...\n";
-
    TStopwatch timer;
 
+   /*** Check if the file with the "good" cascades exists ***/
+   GoodCascade gc[100];
+   Int_t ngood=0;
+   ifstream in("good_cascades");
+   if (in) {
+      cerr<<"Reading good cascades...\n";
+      while (in>>gc[ngood].nlab>>gc[ngood].plab>>
+                gc[ngood].blab>>gc[ngood].code>>
+                 gc[ngood].px>>gc[ngood].py>>gc[ngood].pz>>
+                 gc[ngood].x >>gc[ngood].y >>gc[ngood].z) {
+         ngood++;
+         cerr<<ngood<<'\r';
+         if (ngood==100) {
+            cerr<<"Too many good cascades !\n";
+            break;
+         }
+      }
+      if (!in.eof()) cerr<<"Read error (good_cascades) !\n";
+   } else {
+     /*** generate a file with the "good" cascades ***/
+      cerr<<"Marking good cascades (this will take a while)...\n";
+      ngood=good_cascades(gc,100);
+      ofstream out("good_cascades");
+      if (out) {
+         for (Int_t ngd=0; ngd<ngood; ngd++)            
+            out<<gc[ngd].nlab<<' '<<gc[ngd].plab<<' '<<
+                gc[ngd].blab<<' '<<gc[ngd].code<<' '<<
+                 gc[ngd].px<<' '<<gc[ngd].py<<' '<<gc[ngd].pz<<' '<<
+                 gc[ngd].x <<' '<<gc[ngd].y <<' '<<gc[ngd].z <<endl;
+      } else cerr<<"Can not open file (good_cascades) !\n";
+      out.close();
+   }
+
+   AliESD *event=0;
+   { /*** Load reconstructed vertices ***/
+   TFile *ef=TFile::Open("AliESDcas.root");
+   if ((!ef)||(!ef->IsOpen())) {
+     cerr<<"AliCascadeComparison.C: Can't open AliESDcas.root !\n";
+     return 1;
+   }
+   TKey *key=0;
+   TIter next(ef->GetListOfKeys());
+   if ((key=(TKey*)next())!=0) event=(AliESD*)key->ReadObj();
+   ef->Close();
+   }
+   Int_t nentr=event->GetNumberOfCascades();
+
+
    const Double_t cascadeWindow=0.05, cascadeWidth=0.015; 
    Double_t ptncut=0.12, ptpcut=0.33, kine0cut=0.003;
    Double_t ptbcut=0.11, kinecut=0.002;
@@ -109,77 +156,23 @@ Int_t AliCascadeComparison(Int_t code=3312) {
    TH1F *csf =new TH1F("csf","Fake Cascade Effective Mass",40, mmin, mmax);
    csf->SetXTitle("(GeV)"); csf->SetFillColor(6);
 
-   if (gAlice) {
-      delete gAlice->GetRunLoader();
-      delete gAlice; 
-      gAlice=0;
-   }   
-   AliRunLoader *rl = AliRunLoader::Open("galice.root");
-   if (!rl) {
-       cerr<<"AliV0Comparison.C :Can't start sesion !\n";
-       return 1;
-   }
-   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
-   if (itsl == 0x0) {
-       cerr<<"AliV0Comparison.C : Can not find the ITSLoader\n";
-       delete rl;
-       return 2;
-   }
-
-   /*** Load reconstructed cascades ***/
-   TObjArray carray(1000);
-   itsl->LoadCascades();
-   TTree *xTree=itsl->TreeX();
-   TBranch *branch=xTree->GetBranch("cascades");
-   Int_t nentr=(Int_t)xTree->GetEntries();
-   for (Int_t i=0; i<nentr; i++) {
-       AliCascadeVertex *iovertex=new AliCascadeVertex; 
-       branch->SetAddress(&iovertex);
-       xTree->GetEvent(i);
-       carray.AddLast(iovertex);
-   }
-
-   /*** Check if the file with the "good" cascades exists ***/
-   GoodCascade gc[100];
-   Int_t ngood=0;
-   ifstream in("good_cascades");
-   if (in) {
-      cerr<<"Reading good cascades...\n";
-      while (in>>gc[ngood].nlab>>gc[ngood].plab>>
-                gc[ngood].blab>>gc[ngood].code>>
-                 gc[ngood].px>>gc[ngood].py>>gc[ngood].pz>>
-                 gc[ngood].x >>gc[ngood].y >>gc[ngood].z) {
-         ngood++;
-         cerr<<ngood<<'\r';
-         if (ngood==100) {
-            cerr<<"Too many good cascades !\n";
-            break;
-         }
-      }
-      if (!in.eof()) cerr<<"Read error (good_cascades) !\n";
-   } else {
-     /*** generate a file with the "good" cascades ***/
-      cerr<<"Marking good cascades (this will take a while)...\n";
-      ngood=good_cascades(gc,100);
-      ofstream out("good_cascades");
-      if (out) {
-         for (Int_t ngd=0; ngd<ngood; ngd++)            
-            out<<gc[ngd].nlab<<' '<<gc[ngd].plab<<' '<<
-                gc[ngd].blab<<' '<<gc[ngd].code<<' '<<
-                 gc[ngd].px<<' '<<gc[ngd].py<<' '<<gc[ngd].pz<<' '<<
-                 gc[ngd].x <<' '<<gc[ngd].y <<' '<<gc[ngd].z <<endl;
-      } else cerr<<"Can not open file (good_cascades) !\n";
-      out.close();
-   }
-
    Double_t pxg=0.,pyg=0.,ptg=0.;
    Int_t nlab=-1, plab=-1, blab=-1;
    Int_t i;
    for (i=0; i<nentr; i++) {
-       AliCascadeVertex *cascade=(AliCascadeVertex*)carray.UncheckedAt(i);
-       nlab=TMath::Abs(cascade->GetNindex()); 
-       plab=TMath::Abs(cascade->GetPindex());
-       blab=TMath::Abs(cascade->GetBindex());
+       AliESDcascade *cascade=event->GetCascade(i);
+
+       Int_t nidx=TMath::Abs(cascade->GetNindex());
+       Int_t pidx=TMath::Abs(cascade->GetPindex());
+       Int_t bidx=TMath::Abs(cascade->GetBindex());
+
+       AliESDtrack *ntrack=event->GetTrack(nidx);
+       AliESDtrack *ptrack=event->GetTrack(pidx);
+       AliESDtrack *btrack=event->GetTrack(bidx);
+
+       nlab=TMath::Abs(ntrack->GetLabel()); 
+       plab=TMath::Abs(ptrack->GetLabel());
+       blab=TMath::Abs(btrack->GetLabel());
 
        /** Kinematical cuts **/
        Double_t pxn,pyn,pzn; cascade->GetNPxPyPz(pxn,pyn,pzn); 
@@ -245,7 +238,7 @@ Int_t AliCascadeComparison(Int_t code=3312) {
      cerr<<"Cascade ("<<nlab<<','<<plab<<","<<blab<<") has not been found !\n";
    }
 
-   carray.Delete();
+   delete event;
 
    Stat_t ng=hgood->GetEntries();
    Stat_t nf=hfound->GetEntries();
@@ -331,8 +324,6 @@ Int_t AliCascadeComparison(Int_t code=3312) {
 
    timer.Stop(); timer.Print();
 
-   delete rl;
-
    return 0;
 }
 
@@ -363,18 +354,17 @@ Int_t good_cascades(GoodCascade *gc, Int_t max) {
    }
 
    /*** Get an access to the kinematics ***/
-   AliRunLoader *rl =
-        AliRunLoader::GetRunLoader(AliConfig::fgkDefaultEventFolderName);
-   if (rl == 0x0) {
-  ::Fatal("AliCascadeComparison.C::good_cascades","Can not find Run Loader !");
-   }
-
-   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
-   if (itsl == 0x0) {
-       cerr<<"AliITSComparisonV2.C : Can not find TPCLoader\n";
-       delete rl;
+   if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice; 
+      gAlice=0;
+   }   
+   AliRunLoader *rl = AliRunLoader::Open("galice.root");
+   if (!rl) {
+       cerr<<"AliCascadeComparison.C, good_cascades :Can't start session !\n";
        return 1;
    }
+
    rl->LoadgAlice();
    rl->LoadHeader();
    rl->LoadKinematics();
index 952d88b..c754483 100644 (file)
@@ -1,4 +1,10 @@
 /****************************************************************************
+ *           Very important, delicate and rather obscure macro.             *
+ *                                                                          *
+ *               Creates list of "trackable" tracks,                        *
+ *             calculates efficiency, resolutions etc.                      *
+ *  The ESD tracks must be in an appropriate state (kITSin or kITSrefit)    *
+ *                                                                          *
  *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
  ****************************************************************************/
 
   #include "AliITStrackV2.h"
   #include "AliITSclusterV2.h"
   #include "AliMagF.h"
+  #include "AliESD.h"
 
   #include "TFile.h"
+  #include "TKey.h"
   #include "TTree.h"
   #include "TH1.h"
   #include "TH2.h"
@@ -38,6 +46,9 @@ struct GoodTrackITS {
   Float_t x,y,z;
 };
 
+Int_t good_tracks_its(GoodTrackITS *gt, const Int_t max, 
+const char* evfoldname = AliConfig::fgkDefaultEventFolderName);
+
 extern AliRun *gAlice;
 
 Int_t AliITSComparisonV2() {
@@ -55,28 +66,18 @@ Int_t AliITSComparisonV2() {
    }
    rl->LoadgAlice();
    if (rl->GetAliRun())
-   AliKalmanTrack::
-   SetConvConst(1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField());
+   AliKalmanTrack::SetConvConst(
+     1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
+   );
    else {
       cerr<<"AliITSComparisonV2.C :Can't get AliRun !\n";
       return 1;
    }
    //rl->UnloadgAlice();
-    
-   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
-   if (itsl == 0x0) {
-       cerr<<"AliITSComparisonV2.C : Can not find ITSLoader\n";
-       delete rl;
-       return 3;
-   }
    
-   const Int_t MAX=15000;
-   Int_t good_tracks_its(
-     GoodTrackITS *gt, const Int_t max, 
-     const char* evfoldname = AliConfig::fgkDefaultEventFolderName
-   );//declaration only
 
    /* Generate a list of "good" tracks */
+   const Int_t MAX=15000;
    GoodTrackITS gt[MAX];
    Int_t ngood=0;
    ifstream in("good_tracks_its");
@@ -106,22 +107,28 @@ Int_t AliITSComparisonV2() {
       out.close();
    }
 
-   Int_t nentr=0; TObjArray tarray(2000);
-   {/* Load tracks */ 
-     itsl->LoadTracks();
-     TTree *tracktree=itsl->TreeT();
-     if (!tracktree) {cerr<<"Can't get a tree with ITS tracks !\n"; return 4;}
-     TBranch *tbranch=tracktree->GetBranch("tracks");
-     nentr=(Int_t)tracktree->GetEntries();
-
-     for (Int_t i=0; i<nentr; i++) {
-        AliITStrackV2 *iotrack=new AliITStrackV2;
-        tbranch->SetAddress(&iotrack);
-        tracktree->GetEvent(i);
+   TObjArray tarray(2000);
+   { /*Load tracks*/
+   TFile *ef=TFile::Open("AliESDits.root");
+   if ((!ef)||(!ef->IsOpen())) {
+      ::Fatal("AliITSComparisonV2.C","Can't open AliESDits.root !");
+   }
+   TKey *key=0;
+   TIter next(ef->GetListOfKeys());
+   if ((key=(TKey*)next())!=0) {
+     AliESD *event=(AliESD*)key->ReadObj();
+     Int_t ntrk=event->GetNumberOfTracks();
+     for (Int_t i=0; i<ntrk; i++) {
+        AliESDtrack *t=event->GetTrack(i);
+        if ((t->GetStatus()&AliESDtrack::kITSin)==0) continue;
+        AliITStrackV2 *iotrack=new AliITStrackV2(*t);
         tarray.AddLast(iotrack);
      }
-     itsl->UnloadTracks();
+     delete event;
+   }
+   ef->Close();
    }
+   Int_t nentr=tarray.GetEntriesFast();
 
    TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
    TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
index 5a27e8e..c6f31f7 100644 (file)
   #include <fstream.h>
 
   #include "AliRun.h"
+  #include "AliMC.h"
   #include "AliHeader.h"
   #include "AliRunLoader.h"
   #include "AliITSLoader.h"
 
   #include "TH1.h"
   #include "TFile.h"
+  #include "TKey.h"
   #include "TTree.h"
   #include "TObjArray.h"
   #include "TStyle.h"
@@ -28,7 +30,7 @@
   #include "TStopwatch.h"
 
   #include "AliRun.h"
-  #include "AliPDG.h"
+  #include "AliESD.h"
   #include "AliV0vertex.h"
 #endif
 
@@ -44,9 +46,55 @@ extern AliRun *gAlice;
 
 Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
    cerr<<"Doing comparison...\n";
-
    TStopwatch timer;
 
+   /*** Check if the file with the "good" vertices exists ***/
+   GoodVertex gv[1000];
+   Int_t ngood=0;
+   ifstream in("good_vertices");
+   if (in) {
+      cerr<<"Reading good vertices...\n";
+      while (in>>gv[ngood].nlab>>gv[ngood].plab>>gv[ngood].code>>
+                 gv[ngood].px>>gv[ngood].py>>gv[ngood].pz>>
+                 gv[ngood].x >>gv[ngood].y >>gv[ngood].z) {
+         ngood++;
+         cerr<<ngood<<'\r';
+         if (ngood==1000) {
+            cerr<<"Too many good vertices !\n";
+            break;
+         }
+      }
+      if (!in.eof()) cerr<<"Read error (good_vertices) !\n";
+   } else {
+     /*** generate a file with the "good" vertices ***/
+      cerr<<"Marking good vertices (this will take a while)...\n";
+      ngood=good_vertices(gv,1000);
+      ofstream out("good_vertices");
+      if (out) {
+         for (Int_t ngd=0; ngd<ngood; ngd++)            
+           out<<gv[ngd].nlab<<' '<<gv[ngd].plab<<' '<<gv[ngd].code<<' '<<
+                 gv[ngd].px<<' '<<gv[ngd].py<<' '<<gv[ngd].pz<<' '<<
+                 gv[ngd].x <<' '<<gv[ngd].y <<' '<<gv[ngd].z <<endl;
+      } else cerr<<"Can not open file (good_vertices) !\n";
+      out.close();
+   }
+
+   
+   AliESD *event=0;
+   { /*** Load reconstructed vertices ***/
+   TFile *ef=TFile::Open("AliESDv0.root");
+   if ((!ef)||(!ef->IsOpen())) {
+     cerr<<"AliV0Comparison.C: Can't open AliESDv0.root !\n";
+     return 1;
+   }
+   TKey *key=0;
+   TIter next(ef->GetListOfKeys());
+   if ((key=(TKey*)next())!=0) event=(AliESD*)key->ReadObj();
+   ef->Close();
+   }
+   Int_t nentr=event->GetNumberOfV0s();
+
+
    const Double_t V0window=0.05;
    Double_t ptncut=0.13, ptpcut=0.13, kinecut=0.03; 
    Double_t V0mass=0.497672, V0width=0.020;
@@ -94,74 +142,20 @@ Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
    v0sf->SetXTitle("(GeV)"); v0sf->SetFillColor(6);
 
 
-   if (gAlice) {
-      delete gAlice->GetRunLoader();
-      delete gAlice; 
-      gAlice=0;
-   }   
-   AliRunLoader *rl = AliRunLoader::Open("galice.root");
-   if (!rl) {
-       cerr<<"AliV0Comparison.C :Can't start sesion !\n";
-       return 1;
-   }
-   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
-   if (itsl == 0x0) {
-       cerr<<"AliV0Comparison.C : Can not find the ITSLoader\n";
-       delete rl;
-       return 2;
-   }
-
-   /*** Load reconstructed vertices ***/
-   TObjArray varray(1000);
-   itsl->LoadV0s();
-   TTree *vTree=itsl->TreeV0();
-   TBranch *branch=vTree->GetBranch("vertices");
-   Int_t nentr=(Int_t)vTree->GetEntries();
-   for (Int_t i=0; i<nentr; i++) {
-       AliV0vertex *iovertex=new AliV0vertex; branch->SetAddress(&iovertex);
-       vTree->GetEvent(i);
-       varray.AddLast(iovertex);
-   }
-
-   /*** Check if the file with the "good" vertices exists ***/
-   GoodVertex gv[1000];
-   Int_t ngood=0;
-   ifstream in("good_vertices");
-   if (in) {
-      cerr<<"Reading good vertices...\n";
-      while (in>>gv[ngood].nlab>>gv[ngood].plab>>gv[ngood].code>>
-                 gv[ngood].px>>gv[ngood].py>>gv[ngood].pz>>
-                 gv[ngood].x >>gv[ngood].y >>gv[ngood].z) {
-         ngood++;
-         cerr<<ngood<<'\r';
-         if (ngood==1000) {
-            cerr<<"Too many good vertices !\n";
-            break;
-         }
-      }
-      if (!in.eof()) cerr<<"Read error (good_vertices) !\n";
-   } else {
-     /*** generate a file with the "good" vertices ***/
-      cerr<<"Marking good vertices (this will take a while)...\n";
-      ngood=good_vertices(gv,1000);
-      ofstream out("good_vertices");
-      if (out) {
-         for (Int_t ngd=0; ngd<ngood; ngd++)            
-           out<<gv[ngd].nlab<<' '<<gv[ngd].plab<<' '<<gv[ngd].code<<' '<<
-                 gv[ngd].px<<' '<<gv[ngd].py<<' '<<gv[ngd].pz<<' '<<
-                 gv[ngd].x <<' '<<gv[ngd].y <<' '<<gv[ngd].z <<endl;
-      } else cerr<<"Can not open file (good_vertices) !\n";
-      out.close();
-   }
-
-
    Double_t pxg=0.,pyg=0.,ptg=0.;
    Int_t nlab=-1, plab=-1;
    Int_t i;
    for (i=0; i<nentr; i++) {
-       AliV0vertex *vertex=(AliV0vertex*)varray.UncheckedAt(i);
-       nlab=TMath::Abs(vertex->GetNindex());
-       plab=TMath::Abs(vertex->GetPindex());
+       AliESDv0 *vertex=event->GetV0(i);
+
+       Int_t nidx=TMath::Abs(vertex->GetNindex());
+       Int_t pidx=TMath::Abs(vertex->GetPindex());
+
+       AliESDtrack *ntrack=event->GetTrack(nidx);
+       AliESDtrack *ptrack=event->GetTrack(pidx);
+
+       nlab=TMath::Abs(ntrack->GetLabel());
+       plab=TMath::Abs(ptrack->GetLabel());
 
        /** Kinematical cuts **/
        Double_t pxn,pyn,pzn; vertex->GetNPxPyPz(pxn,pyn,pzn); 
@@ -221,7 +215,7 @@ Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
       cerr<<"Vertex ("<<nlab<<','<<plab<<") has not been found !\n";
    }
 
-   varray.Delete();
+   delete event;
 
    Stat_t ng=hgood->GetEntries();
    Stat_t nf=hfound->GetEntries();
@@ -305,8 +299,6 @@ Int_t AliV0Comparison(Int_t code=310) { //Lambda=3122, LambdaBar=-3122
 
    timer.Stop(); timer.Print();
 
-   delete rl;
-
    return 0;
 }
 
@@ -336,18 +328,17 @@ Int_t good_vertices(GoodVertex *gv, Int_t max) {
    }
 
    /*** Get an access to the kinematics ***/
-   AliRunLoader *rl = 
-        AliRunLoader::GetRunLoader(AliConfig::fgkDefaultEventFolderName);
-   if (rl == 0x0) {
-     ::Fatal("AliV0Comparison.C::good_vertices","Can not find Run Loader !");
-   }
-
-   AliITSLoader* itsl = (AliITSLoader*)rl->GetLoader("ITSLoader");
-   if (itsl == 0x0) {
-       cerr<<"AliITSComparisonV2.C : Can not find TPCLoader\n";
-       delete rl;
+   if (gAlice) {
+      delete gAlice->GetRunLoader();
+      delete gAlice; 
+      gAlice=0;
+   }   
+   AliRunLoader *rl = AliRunLoader::Open("galice.root");
+   if (!rl) {
+       cerr<<"AliV0Comparison.C, good_vertices :Can't start session !\n";
        return 1;
    }
+
    rl->LoadgAlice();
    rl->LoadHeader();
    rl->LoadKinematics();
index 632b575..283ea93 100644 (file)
@@ -3,18 +3,43 @@
  *                                                                          *
  *               Creates list of "trackable" tracks,                        *
  *             calculates efficiency, resolutions etc.                      *
+ *  The ESD tracks must be in an appropriate state (kTPCin or kTPCrefit)    *
  *                                                                          *
  *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
  * with several nice improvements by: M.Ivanov, GSI, m.ivanov@gsi.de        *
  ****************************************************************************/
 
 #if !defined(__CINT__) || defined(__MAKECINT__)
-  #include "alles.h"
-  #include "AliTPCtracker.h"
+  #include <TMath.h>
+  #include <TError.h>
+  #include <Riostream.h>
+  #include <TH1F.h>
+  #include <TH2F.h>
+  #include <TTree.h>
+  #include <TParticle.h>
+  #include <TPad.h>
+  #include <TCanvas.h>
+  #include <TLine.h>
+  #include <TText.h>
+  #include <TBenchmark.h>
+  #include <TStyle.h>
+  #include <TKey.h>
+
+  #include "AliStack.h"
+  #include "AliHeader.h"
+  #include "AliTrackReference.h"
   #include "AliRunLoader.h"
-  #include "AliTPCLoader.h"
   #include "AliMagF.h"
   #include "AliRun.h"
+  #include "AliESD.h"
+
+  #include "AliSimDigits.h"
+  #include "AliTPC.h"
+  #include "AliTPCParamSR.h"
+  #include "AliTPCClustersArray.h"
+  #include "AliTPCClustersRow.h"
+  #include "AliTPCtracker.h"
+  #include "AliTPCLoader.h"
 #endif
 
 struct GoodTrackTPC {
@@ -24,12 +49,16 @@ struct GoodTrackTPC {
   Float_t x,y,z;
 };
 
-Int_t 
-good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname);
+Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, 
+const char* evfoldname=AliConfig::fgkDefaultEventFolderName);
 
 extern AliRun *gAlice;
+extern TBenchmark *gBenchmark;
 
-Int_t AliTPCComparison(Bool_t thcut=1.0) {
+Int_t AliTPCComparison() {
+   gBenchmark->Start("AliTPCComparison");
+
+   ::Info("AliTPCComparison.C","Doing comparison...");
 
    if (gAlice) { 
        delete gAlice->GetRunLoader();
@@ -37,20 +66,9 @@ Int_t AliTPCComparison(Bool_t thcut=1.0) {
        gAlice = 0x0;
     }
 
-   cerr<<"Doing comparison...\n";
-
-   const Int_t MAX=20000;
-   Int_t good_tracks_tpc(
-      GoodTrackTPC *gt, const Int_t max, 
-      const char* evfoldname = AliConfig::fgkDefaultEventFolderName
-   );//declaration only
-
-   gBenchmark->Start("AliTPCComparison");
-
    AliRunLoader *rl = AliRunLoader::Open("galice.root","COMPARISON");
    if (!rl) {
-       cerr<<"Can't start sesion !\n";
-       return 1;
+       ::Fatal("AliTPCComparison.C","Can't start session !");
    }
 
    rl->LoadgAlice();
@@ -59,49 +77,32 @@ Int_t AliTPCComparison(Bool_t thcut=1.0) {
      1000/0.299792458/rl->GetAliRun()->Field()->SolenoidField()
    );
    else {
-       cerr<<"AliTPCComparison.C :Can't get AliRun !\n";
-       return 1;
+       ::Fatal("AliTPCComparison.C","Can't get AliRun !");
    }
    //rl->UnloadgAlice();
   
-   AliTPCLoader * tpcl = (AliTPCLoader *)rl->GetLoader("TPCLoader");
-   if (tpcl == 0x0) {
-       cerr<<"AliTPCComparison.C : Can not find TPCLoader\n";
-       delete rl;
-       return 3;
-   }
 
    /* Generate a list of "good" tracks */
-
+   const Int_t MAX=20000;
    GoodTrackTPC gt[MAX];
    Int_t ngood=0;
    ifstream in("good_tracks_tpc");
    if (in) {
-      cerr<<"Reading good tracks...\n";
+      ::Info("AliTPCComparison.C","Reading good tracks...");
       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) {
-       Double_t rin = TMath::Sqrt( gt[ngood].x*gt[ngood].x+gt[ngood].y*gt[ngood].y);
-       if (rin<1) continue;
-       Double_t theta = gt[ngood].z/rin;
-       //theta =0;
-       if (TMath::Abs(theta)<thcut){ 
          ngood++;
-         cerr<<ngood<<'\r';
          if (ngood==MAX) {
-           cerr<<"Too many good tracks !\n";
+           ::Warning("AliTPCComparison.C","Too many good tracks !");
            break;
          }     
-       //ngood++;
-        // cerr<<ngood<<'\r';
-        // if (ngood==MAX) {
-        //    cerr<<"Too many good tracks !\n";
-        //    break;
-       }
       }
-      if (!in.eof()) cerr<<"Read error (good_tracks_tpc) !\n";
+      if (!in.eof()) 
+        ::Fatal("AliTPCComparison.C","Read error (good_tracks_tpc) !");
    } else {
-     cerr<<"Marking good tracks (this will take a while)...\n";
+      ::Info
+      ("AliTPCComparison","Marking good tracks (this will take a while)...");
       ngood=good_tracks_tpc(gt,MAX,"COMPARISON");
       ofstream out("good_tracks_tpc");
       if (out) {
@@ -109,29 +110,33 @@ Int_t AliTPCComparison(Bool_t thcut=1.0) {
            out<<gt[ngd].lab<<' '<<gt[ngd].code<<' '<<
                  gt[ngd].px<<' '<<gt[ngd].py<<' '<<gt[ngd].pz<<' '<<
                  gt[ngd].x <<' '<<gt[ngd].y <<' '<<gt[ngd].z <<endl;
-      } else cerr<<"Can not open file (good_tracks_tpc) !\n";
+      } else 
+         ::Fatal("AliTPCComparison.C","Can not open file (good_tracks_tpc) !");
       out.close();
    }
 
-   Int_t nentr=0,i=0; TObjArray tarray(MAX);
+   Int_t i=0; TObjArray tarray(MAX);
    { /*Load tracks*/
-   
-     tpcl->LoadTracks();
-     
-     TTree *tracktree=tpcl->TreeT();
-     if (!tracktree) {cerr<<"Can't get a tree with TPC tracks !\n"; return 4;}
-
-     TBranch *tbranch=tracktree->GetBranch("tracks");
-     nentr=(Int_t)tracktree->GetEntries();
-     AliTPCtrack *iotrack=0;
-     for (i=0; i<nentr; i++) {
-       iotrack=new AliTPCtrack;
-       tbranch->SetAddress(&iotrack);
-       tracktree->GetEvent(i);
-       tarray.AddLast(iotrack);
-     }   
-     tpcl->UnloadTracks();
+   TFile *ef=TFile::Open("AliESDtpc.root");
+   if ((!ef)||(!ef->IsOpen())) {
+      ::Fatal("AliTPCComparison.C","Can't open AliESDtpc.root !");
+   }
+   TKey *key=0;
+   TIter next(ef->GetListOfKeys());
+   if ((key=(TKey*)next())!=0) {
+     AliESD *event=(AliESD*)key->ReadObj();
+     Int_t ntrk=event->GetNumberOfTracks();
+     for (Int_t i=0; i<ntrk; i++) {
+        AliESDtrack *t=event->GetTrack(i);
+        if ((t->GetStatus()&AliESDtrack::kTPCin)==0) continue;
+        AliTPCtrack *iotrack=new AliTPCtrack(*t);
+        tarray.AddLast(iotrack);
+     }
+     delete event;
+   }
+   ef->Close();
    }
+   Int_t nentr=tarray.GetEntriesFast();
 
    TH1F *hp=new TH1F("hp","PHI resolution",50,-20.,20.); hp->SetFillColor(4);
    TH1F *hl=new TH1F("hl","LAMBDA resolution",50,-20,20);hl->SetFillColor(4);
@@ -140,12 +145,12 @@ Int_t AliTPCComparison(Bool_t thcut=1.0) {
    TH1F *hmpt=new TH1F("hmpt","Relative Pt resolution (pt>4GeV/c)",30,-60,60); 
    hmpt->SetFillColor(6);
 
-   TH1F *hgood=new TH1F("hgood","Good tracks",30,0.1,6.1);    
-   TH1F *hfound=new TH1F("hfound","Found tracks",30,0.1,6.1);
-   TH1F *hfake=new TH1F("hfake","Fake tracks",30,0.1,6.1);
-   TH1F *hg=new TH1F("hg","",30,0.1,6.1); //efficiency for good tracks
+   TH1F *hgood=new TH1F("hgood","Good tracks",30,0.2,6.1);    
+   TH1F *hfound=new TH1F("hfound","Found tracks",30,0.2,6.1);
+   TH1F *hfake=new TH1F("hfake","Fake tracks",30,0.2,6.1);
+   TH1F *hg=new TH1F("hg","Efficiency for good tracks",30,0.2,6.1);
    hg->SetLineColor(4); hg->SetLineWidth(2);
-   TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,0.1,6.1);
+   TH1F *hf=new TH1F("hf","Efficiency for fake tracks",30,0.2,6.1);
    hf->SetFillColor(1); hf->SetFillStyle(3013); hf->SetLineWidth(2);
 
    TH1F *he =new TH1F("he","dE/dX for pions with 0.4<p<0.5 GeV/c",50,0.,100.);
@@ -337,9 +342,7 @@ good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname) {
    }
    AliTPCLoader *tpcl = (AliTPCLoader *)rl->GetLoader("TPCLoader");
    if (tpcl == 0x0) {
-       cerr<<"AliTPCHits2Digits.C : Can not find TPCLoader\n";
-       delete rl;
-       return 0;
+       ::Fatal("AliTPCHits2Digits.C","Can not find TPCLoader !");
    }
    
    rl->LoadgAlice();
@@ -350,7 +353,9 @@ good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname) {
 
    rl->CdGAFile();
    AliTPCParamSR *digp=(AliTPCParamSR*)gDirectory->Get("75x40_100x60_150x60");
-   if (!digp) { cerr<<"TPC parameters have not been found !\n"; exit(6); }
+   if (!digp) { 
+     ::Fatal("AliTPCHits2Digits.C","TPC parameters have not been found !"); 
+   }
    TPC->SetParam(digp);
 
    rl->LoadHeader();
@@ -458,9 +463,7 @@ good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname) {
      }
       break;
    default:
-      cerr<<"Invalid TPC version !\n";
-      delete rl;
-      exit(7);
+      ::Fatal("AliTPCComparison.C","Invalid TPC version !");
    }
 
    rl->LoadKinematics();
@@ -510,42 +513,38 @@ good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const char* evfoldname) {
 
    /** check if there is also information at the entrance of the TPC **/
    
-   tpcl->LoadHits();
-   TTree *TH=tpcl->TreeH();
-   TPC->SetTreeAddress();
-   np=(Int_t)TH->GetEntries();
+   rl->LoadTrackRefs();
+   TTree *TR=rl->TreeTR();
+   TBranch *branch=TR->GetBranch("TPC");
+   if (branch==0) {
+      ::Fatal("AliTPCComparison.C::good_tracks_tpc","No track references !");
+   }
+   TClonesArray *references=new TClonesArray("AliTrackReference",1000);
+   branch->SetAddress(&references);
+   np=(Int_t)TR->GetEntries();
    for (i=0; i<np; i++) {
-      TPC->ResetHits();
-      TH->GetEvent(i);
-      AliTPChit *phit = (AliTPChit*)TPC->FirstHit(-1);
-      for ( ; phit; phit=(AliTPChit*)TPC->NextHit() ) {
-        if (phit->fQ !=0. ) continue;
-
-         Double_t px=phit->X(), py=phit->Y(), pz=phit->Z();
-
-         if ((phit=(AliTPChit*)TPC->NextHit())==0) break;
-         if (phit->fQ != 0.) continue;
-
-         Double_t x=phit->X(), y=phit->Y(), z=phit->Z();
-        if (TMath::Sqrt(x*x+y*y)>90.) continue;
-
-         Int_t j, lab=phit->Track();
-         for (j=0; j<nt; j++) {if (gt[j].lab==lab) break;}
-         if (j==nt) continue;         
-
-        // (px,py,pz) - in global coordinate system, (x,y,z) - in local !
-         gt[j].px=px; gt[j].py=py; gt[j].pz=pz;
-         Float_t cs,sn; digp->AdjustCosSin(phit->fSector,cs,sn);
-         gt[j].x = x*cs + y*sn;
-         gt[j].y =-x*sn + y*cs;
-         gt[j].z = z;
-      }
-      cerr<<np-i<<"                \r";
+      references->Clear();
+      TR->GetEvent(i);
+      Int_t nref=references->GetEntriesFast();
+      if (nref==0) continue;
+      AliTrackReference *ref=(AliTrackReference*)references->UncheckedAt(0);
+
+      Int_t j, lab=ref->Label();
+      for (j=0; j<nt; j++) {if (gt[j].lab==lab) break;}
+      if (j==nt) continue;         
+
+      // (px,py,pz) - in global coordinate system, (x,y,z) - in local !
+      gt[j].px=ref->Px(); gt[j].py=ref->Py(); gt[j].pz=ref->Pz();
+      gt[j].x = ref->LocalX();
+      gt[j].y = ref->LocalY();
+      gt[j].z = ref->Z();
    }
 
+   delete references;
    delete[] good;
    
    tpcl->UnloadHits();
+   rl->UnloadTrackRefs();
    rl->UnloadgAlice();
 
    gBenchmark->Stop("AliTPCComparison");