]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - HBTAN/AliHBTReaderTPC.cxx
Comment added
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx
index 394baa75385f05373f3cbc1c92919d91fbf44219..cc17bc1f97cf747e2b7a9f8f222f882f3aa9c8fa 100644 (file)
 #include "AliHBTReaderTPC.h"
-
-#include <iostream.h>
-//#include <fstream.h>
+//______________________________________________
+//
+// class AliHBTReaderTPC
+//
+// reader for TPC tracks
+// needs galice.root
+// just to shut up coding conventions checker
+// 
+// more info: http://aliweb.cern.ch/people/skowron/analyzer/index.html
+// Piotr.Skowronski@cern.ch
+//
+///////////////////////////////////////////////////////////////////////////
 #include <TTree.h>
-#include <TFile.h>
 #include <TParticle.h>
+#include <TH1.h>
+
 
 #include <AliRun.h>
+#include <AliLoader.h>
+#include <AliStack.h>
 #include <AliMagF.h>
 #include <AliTPCtrack.h>
-#include <AliTPCParam.h>
-#include <AliTPCtracker.h>
+#include <AliTPCLoader.h>
 
-#include "AliHBTRun.h"
 #include "AliHBTEvent.h"
 #include "AliHBTParticle.h"
-#include "AliHBTParticleCut.h"
+#include "AliHBTTrackPoints.h"
+#include "AliHBTClusterMap.h"
 
 
 ClassImp(AliHBTReaderTPC)
-//reader for TPC tracking
-//needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc 
-//
-//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
-//Piotr.Skowronski@cern.ch
-
-AliHBTReaderTPC::
- AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
-                 const Char_t* galicefilename):
-                 fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
-                 fGAliceFileName(galicefilename)
+
+AliHBTReaderTPC::AliHBTReaderTPC():
+ fFileName("galice.root"),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE),
+ fNTrackPoints(0),
+ fdR(0.0),
+ fClusterMap(kFALSE),
+ fNClustMin(0),
+ fNClustMax(150),
+ fNChi2PerClustMin(0.0),
+ fNChi2PerClustMax(10e5),
+ fC00Min(0.0),
+ fC00Max(10e5),
+ fC11Min(0.0),
+ fC11Max(10e5),
+ fC22Min(0.0),
+ fC22Max(10e5),
+ fC33Min(0.0),
+ fC33Max(10e5),
+ fC44Min(0.0),
+ fC44Max(10e5)
+{
+  //constructor
+  
+}
+/********************************************************************/
+
+AliHBTReaderTPC::AliHBTReaderTPC(const Char_t* galicefilename):
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE),
+ fNTrackPoints(0),
+ fdR(0.0),
+ fNClustMin(0),
+ fNClustMax(150),
+ fNChi2PerClustMin(0.0),
+ fNChi2PerClustMax(10e5),
+ fC00Min(0.0),
+ fC00Max(10e5),
+ fC11Min(0.0),
+ fC11Max(10e5),
+ fC22Min(0.0),
+ fC22Max(10e5),
+ fC33Min(0.0),
+ fC33Max(10e5),
+ fC44Min(0.0),
+ fC44Max(10e5)
 {
   //constructor, 
   //Defaults:
-  //  trackfilename = "AliTPCtracks.root"
-  //  clusterfilename = "AliTPCclusters.root"
-  //  galicefilename = ""  - this means: Do not open gAlice file - 
-  //                         just leave the global pointer untached
-  
-  fParticles = new AliHBTRun();
-  fTracks    = new AliHBTRun();
-  fIsRead = kFALSE;
 }
 /********************************************************************/
-AliHBTReaderTPC::
-AliHBTReaderTPC(TObjArray* dirs,
-                  const Char_t* trackfilename = "AliTPCtracks.root",
-                  const Char_t* clusterfilename = "AliTPCclusters.root",
-                  const Char_t* galicefilename = "galice.root"):AliHBTReader(dirs),
-                  fTrackFileName(trackfilename),
-                 fClusterFileName(clusterfilename),fGAliceFileName(galicefilename)
 
+AliHBTReaderTPC::AliHBTReaderTPC(TObjArray* dirs, const Char_t* galicefilename):
+ AliHBTReader(dirs), 
+ fFileName(galicefilename),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(0.0),
+ fUseMagFFromRun(kTRUE),
+ fNTrackPoints(0),
+ fdR(0.0),
+ fClusterMap(kFALSE),
+ fNClustMin(0),
+ fNClustMax(150),
+ fNChi2PerClustMin(0.0),
+ fNChi2PerClustMax(10e5),
+ fC00Min(0.0),
+ fC00Max(10e5),
+ fC11Min(0.0),
+ fC11Max(10e5),
+ fC22Min(0.0),
+ fC22Max(10e5),
+ fC33Min(0.0),
+ fC33Max(10e5),
+ fC44Min(0.0),
+ fC44Max(10e5)
 {
   //constructor, 
   //Defaults:
-  //  trackfilename = "AliTPCtracks.root"
-  //  clusterfilename = "AliTPCclusters.root"
   //  galicefilename = ""  - this means: Do not open gAlice file - 
   //                         just leave the global pointer untached
   
-  fParticles = new AliHBTRun();
-  fTracks    = new AliHBTRun();
-  fIsRead = kFALSE;
-  
+}
+/********************************************************************/
+AliHBTReaderTPC::AliHBTReaderTPC(const AliHBTReaderTPC& in):
+ AliHBTReader(in),
+ fFileName(in.fFileName),
+ fRunLoader(0x0),
+ fTPCLoader(0x0),
+ fMagneticField(in.fMagneticField),
+ fUseMagFFromRun(in.fUseMagFFromRun),
+ fNTrackPoints(in.fNTrackPoints),
+ fdR(in.fdR),
+ fNClustMin(in.fNClustMin),
+ fNClustMax(in.fNClustMax),
+ fNChi2PerClustMin(in.fNChi2PerClustMin),
+ fNChi2PerClustMax(in.fNChi2PerClustMax),
+ fC00Min(in.fC00Min),
+ fC00Max(in.fC00Max),
+ fC11Min(in.fC11Min),
+ fC11Max(in.fC11Max),
+ fC22Min(in.fC22Min),
+ fC22Max(in.fC22Max),
+ fC33Min(in.fC33Min),
+ fC33Max(in.fC33Max),
+ fC44Min(in.fC44Min),
+ fC44Max(in.fC44Max)
+{
+  //cpy constructor, 
 }
 /********************************************************************/
 
 AliHBTReaderTPC::~AliHBTReaderTPC()
- {
+{
  //desctructor
-   delete fParticles;
-   delete fTracks;
- }
+   if (AliHBTParticle::GetDebug()) 
+    {
+      Info("~AliHBTReaderTPC","deleting Run Loader");
+      AliLoader::SetDebug(AliHBTParticle::GetDebug());
+    }
+   
+   delete fRunLoader;
+   
+   if (AliHBTParticle::GetDebug()) 
+    {
+      Info("~AliHBTReaderTPC","deleting Run Loader Done");
+    }
+}
 /********************************************************************/
 
-AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
- {
- //returns Nth event with simulated particles
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
-     {
-       Error("GetParticleEvent","Error in reading");
-       return 0x0;
-     }
-   return fParticles->GetEvent(n);
- }
+AliHBTReaderTPC& AliHBTReaderTPC::operator=(const AliHBTReaderTPC& in)
+{
+//Assigment operator
+
+ delete fRunLoader;
+
+ fFileName = in.fFileName;
+ fRunLoader = 0x0;
+ fTPCLoader = 0x0;
+ fMagneticField = in.fMagneticField;
+ fUseMagFFromRun = in.fUseMagFFromRun;
+ fNTrackPoints = in.fNTrackPoints;
+ fdR = in.fdR;
+ fNClustMin = in.fNClustMin;
+ fNClustMax = in.fNClustMax;
+ fNChi2PerClustMin = in.fNChi2PerClustMin;
+ fNChi2PerClustMax = in.fNChi2PerClustMax;
+ fC00Min = in.fC00Min;
+ fC00Max = in.fC00Max;
+ fC11Min = in.fC11Min;
+ fC11Max = in.fC11Max;
+ fC22Min = in.fC22Min;
+ fC22Max = in.fC22Max;
+ fC33Min = in.fC33Min;
+ fC33Max = in.fC33Max;
+ fC44Min = in.fC44Min;
+ fC44Max = in.fC44Max;
+ return *this; 
+} 
+/********************************************************************/
+
+void AliHBTReaderTPC::Rewind()
+{
+//Rewind reading to the beginning
+  delete fRunLoader;
+  fRunLoader = 0x0;
+  fCurrentDir = 0;
+  fNEventsRead= 0;
+  if (fTrackCounter) fTrackCounter->Reset();
+}
 /********************************************************************/
-AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
+
+Int_t AliHBTReaderTPC::ReadNext()
  {
- //returns Nth event with reconstructed tracks
-   if (!fIsRead) 
-    if(Read(fParticles,fTracks))
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+  Info("Read","");
+  
+  TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
+  tarray->SetOwner(); //set the ownership of the objects it contains
+                      //when array is is deleted or cleared all objects 
+                      //that it contains are deleted
+
+  if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
+  if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
+
+  fParticlesEvent->Reset();
+  fTracksEvent->Reset();
+
+  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
+   {
+      
+    if (fRunLoader == 0x0) 
+      if (OpenNextSession()) continue;//directory counter is increased inside in case of error
+
+    if (fCurrentEvent == fRunLoader->GetNumberOfEvents())
      {
-       Error("GetTrackEvent","Error in reading");
-       return 0x0;
+       //read next directory
+       delete fRunLoader;//close current session
+       fRunLoader = 0x0;//assure pointer is null
+       fCurrentDir++;//go to next dir
+       continue;//directory counter is increased inside in case of error
      }
-   return fTracks->GetEvent(n);
- }
-/********************************************************************/
+     
+    Info("ReadNext","Reading Event %d",fCurrentEvent);
+    
+    fRunLoader->GetEvent(fCurrentEvent);
+    TTree *tracktree = fTPCLoader->TreeT();//get the tree 
+    if (!tracktree) //check if we got the tree
+      {//if not return with error
+         Error("ReadNext","Can't get a tree with TPC tracks !\n"); 
+         fCurrentEvent++;//go to next dir
+         continue;
+      }
+   
+    TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
+    if (!trackbranch) ////check if we got the branch
+      {//if not return with error
+        Error("ReadNext","Can't get a branch with TPC tracks !\n"); 
+        fCurrentEvent++;//go to next dir
+        continue;
+      }
+    Int_t ntpctracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
+    Info("ReadNext","Found %d TPC tracks.",ntpctracks);
+    //Copy tracks to array
 
-Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
- {
- //returns number of events of particles
-   if (!fIsRead) 
-    if ( Read(fParticles,fTracks))
+    AliTPCtrack *iotrack=0;
+
+    for (Int_t i=0; i<ntpctracks; i++) //loop over all tpc tracks
      {
-       Error("GetNumberOfPartEvents","Error in reading");
-       return 0;
+       iotrack=new AliTPCtrack;   //create new tracks
+       trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
+       tracktree->GetEvent(i); //stream track i to the iotrack
+       tarray->AddLast(iotrack); //put the track in the array
      }
-   return fParticles->GetNumberOfEvents();
- }
 
-/********************************************************************/
-Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
- {
- //returns number of events of tracks
-  if (!fIsRead)
-    if(Read(fParticles,fTracks))
+    Double_t xk;
+    Double_t par[5];
+    Float_t phi, lam, pt;//angles and transverse momentum
+    Int_t label; //label of the current track
+
+    AliStack* stack = fRunLoader->Stack();
+    if (stack == 0x0)
      {
-       Error("GetNumberOfTrackEvents","Error in reading");
-       return 0;
+       Error("ReadNext","Can not get stack for current event",fCurrentEvent);
+       fCurrentEvent++;
+       continue;
      }
-  return fTracks->GetNumberOfEvents();
+    stack->Particles();
+
+    for (Int_t i=0; i<ntpctracks; i++) //loop over all good tracks
+     { 
+       iotrack = (AliTPCtrack*)tarray->At(i);
+       label = iotrack->GetLabel();
+
+       if (label < 0) continue;
+       
+       if (CheckTrack(iotrack)) continue;//Checks the cuts on track parameters cov. mtx etc
+       
+       TParticle *p = (TParticle*)stack->Particle(label);
+
+       if(p == 0x0) continue; //if returned pointer is NULL
+       if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
+
+       if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
+                                   //if not take next partilce
+
+       AliHBTParticle* part = new AliHBTParticle(*p,i);
+       if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
+                                               //if it does not delete it and take next good track
+
+//       iotrack->PropagateTo(3.,0.0028,65.19);
+//       iotrack->PropagateToVertex(36.66,1.2e-3);//orig
+       iotrack->PropagateToVertex(50.,0.0353);
+       
+       iotrack->GetExternalParameters(xk,par);     //get properties of the track
+       
+       phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
+       if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
+       if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
+       lam=par[3]; 
+       pt=1.0/TMath::Abs(par[4]);
+
+       Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
+       Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
+       Double_t tpz = pt * lam; //track z coordinate of momentum
+
+       Double_t mass = p->GetMass();
+       Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
+
+       AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(),i, tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
+       if(Pass(track))//check if meets all criteria of any of our cuts
+                    //if it does not delete it and take next good track
+        { 
+          delete track;
+          delete part;
+          continue;
+        }
+        
+       if (fNTrackPoints > 0) 
+        {
+          AliHBTTrackPoints* tpts = new AliHBTTrackPoints(fNTrackPoints,iotrack,fdR);
+          track->SetTrackPoints(tpts);
+        }
+       if (  fClusterMap ) 
+        {
+          AliHBTClusterMap* cmap = new AliHBTClusterMap(iotrack);
+          track->SetClusterMap(cmap);
+        }
+    
+       fParticlesEvent->AddParticle(part);
+       fTracksEvent->AddParticle(track);
+      }
+      
+    Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
+            fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
+            fNEventsRead,fCurrentEvent,fCurrentDir);
+    fTrackCounter->Fill(fTracksEvent->GetNumberOfParticles());
+    fCurrentEvent++;
+    fNEventsRead++;
+    delete tarray;
+    return 0;
+   }while(fCurrentDir < GetNumberOfDirs());
+
+  delete tarray;
+  return 1;
  }
 /********************************************************************/
 
-
-Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
- {
- //reads data and puts put to the particles and tracks objects
- //reurns 0 if everything is OK
- //
-  cout<<"AliHBTReaderTPC::Read()"<<endl;
-  Int_t i; //iterator and some temprary values
-  Int_t Nevents = 0;
-  Int_t totalNevents = 0;
-  TFile *aTracksFile;//file with tracks
-  TFile *aClustersFile;//file with clusters
-  TFile *aGAliceFile;//!ile name with galice
-
-  if (!particles) //check if an object is instatiated
+Int_t AliHBTReaderTPC::OpenNextSession()
+{
+//Opens session thats from fCurrentDir 
+  TString filename = GetDirName(fCurrentDir);
+  if (filename.IsNull())
    {
-     Error("Read"," particles object must instatiated before passing it to the reader");
+     DoOpenError("Can not get directory name");
+     return 1;
    }
-  if (!tracks)  //check if an object is instatiated
+  filename = filename +"/"+ fFileName;
+  fRunLoader = AliRunLoader::Open(filename,AliConfig::fgkDefaultEventFolderName);
+  if( fRunLoader == 0x0)
    {
-     Error("Read"," tracks object must instatiated before passing it to the reader");
+     DoOpenError("Can not open session.");
+     return 1;
    }
-  particles->Reset();//clear runs == delete all old events
-  tracks->Reset();
-  TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
-  tarray->SetOwner(); //set the ownership of the objects it contains
-                      //when array is is deleted or cleared all objects 
-                      //that it contains are deleted
-  Int_t currentdir = 0;
 
-  Int_t Ndirs;
-  if (fDirs) //if array with directories is supplied by user
+  fRunLoader->LoadHeader();
+  if (fRunLoader->GetNumberOfEvents() <= 0)
    {
-     Ndirs = fDirs->GetEntries(); //get the number if directories
+     DoOpenError("There is no events in this directory.");
+     return 1;
    }
-  else
+
+  if (fRunLoader->LoadKinematics())
    {
-     Ndirs = 0; //if the array is not supplied read only from current directory
+     DoOpenError("Error occured while loading kinematics.");
+     return 1;
    }
-
-  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
+  fTPCLoader = (AliTPCLoader*)fRunLoader->GetLoader("TPCLoader");
+  if ( fTPCLoader == 0x0)
    {
-    
-    if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
-     {
-       Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
-       return i;
-     }
-  
-    
-    if (gAlice->TreeE())//check if tree E exists
-     {
-      Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
-      cout<<"________________________________________________________\n";
-      cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
-      cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
-      AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());  
-     }
-    else
-     {//if not return an error
-       Error("Read","Can not find Header tree (TreeE) in gAlice");
-       return 1;
-     }
-  
-    aClustersFile->cd();//set cluster file active 
-    AliTPCParam *TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60");
-    if (!TPCParam) 
-      { 
-       Error("Read","TPC parameters have not been found !\n");
-       return 1;
+     DoOpenError("Exiting due to problems with opening files.");
+     return 1;
+   }
+  Info("OpenNextSession","________________________________________________________");
+  Info("OpenNextSession","Found %d event(s) in directory %s",
+        fRunLoader->GetNumberOfEvents(),GetDirName(fCurrentDir).Data());
+  Float_t mf;
+  if (fUseMagFFromRun)
+   {
+     fRunLoader->LoadgAlice();
+     mf = fRunLoader->GetAliRun()->Field()->SolenoidField();
+     Info("OpenNextSession","Setting Magnetic Field from run: B=%fT",mf/10.);
+     fRunLoader->UnloadgAlice();
+   }
+  else
+   {
+     Info("OpenNextSession","Setting Own Magnetic Field: B=%fT",fMagneticField);
+     if (fMagneticField == 0x0)
+      {
+        Fatal("OpenNextSession","Magnetic field can not be 0.");
+        return 1;//pro forma
       }
+     mf = fMagneticField*10.;
+   }
+  AliKalmanTrack::SetConvConst(1000/0.299792458/mf);
 
+
+  if (fTPCLoader->LoadTracks())
+   {
+     DoOpenError("Error occured while loading TPC tracks.");
+     return 1;
+   }
   
-    for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
-     {
-       cout<<"Reading Event "<<currentEvent<<endl;
-       /**************************************/
-        /**************************************/
-         /**************************************/ 
-         
-         aTracksFile->cd();//set track file active
-          
-         Char_t  treename[100];
-         sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
-   
-         TTree *tracktree=0;
-         
-         tracktree=(TTree*)aTracksFile->Get(treename);//get the tree 
-         if (!tracktree) //check if we got the tree
-          {//if not return with error
-            Error("Read","Can't get a tree with TPC tracks !\n"); 
-            
-            return 1;
-          }
-   
-         TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
-         if (!trackbranch) ////check if we got the branch
-          {//if not return with error
-            Error("Read","Can't get a branch with TPC tracks !\n"); 
-            return 2;
-          }
-         Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks 
-         cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
-         //Copy tracks to array
-         
-         AliTPCtrack *iotrack=0;
-         
-         aClustersFile->cd();//set cluster file active 
-         AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
-         if (!tracker) //check if it has created succeffuly
-          {//if not return with error
-            Error("Read","Can't get a tracker !\n"); 
-            return 3;
-          }
-         tracker->LoadInnerSectors();
-         tracker->LoadOuterSectors();
-   
-         for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
-          {
-            iotrack=new AliTPCtrack;   //create new tracks
-            trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
-            tracktree->GetEvent(i); //stream track i to the iotrack
-            tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
-                                             //which is the label of corresponding simulated particle 
-            tarray->AddLast(iotrack); //put the track in the array
-          }
-         
-         aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
-         aTracksFile->Delete("tracks");//delete branch from memmory
-         delete tracker; //delete tracker
-         
-         tracker = 0x0;
-         trackbranch = 0x0;
-         tracktree = 0x0;
-   
-         Double_t xk;
-         Double_t par[5];
-         Float_t phi, lam, pt;//angles and transverse momentum
-         Int_t label; //label of the current track
-         
-         aGAliceFile->cd();
-         gAlice->GetEvent(currentEvent); 
-
-         gAlice->Particles();
-         
-         for (i=0; i<NTPCtracks; i++) //loop over all good tracks
-          { 
-            iotrack = (AliTPCtrack*)tarray->At(i);
-            label = iotrack->GetLabel();
-
-            if (label < 0) continue;
-            
-            TParticle *p = (TParticle*)gAlice->Particle(label);
-            
-            if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type 
-                                        //if not take next partilce
-            
-            AliHBTParticle* part = new AliHBTParticle(*p);
-            if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
-                                                    //if it does not delete it and take next good track
-         
-         
-            iotrack->PropagateToVertex();
-
-            iotrack->GetExternalParameters(xk,par);     //get properties of the track
-            phi=TMath::ASin(par[2]) + iotrack->GetAlpha(); 
-            if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
-            if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
-            lam=par[3]; 
-            pt=1.0/TMath::Abs(par[4]);
-            
-            Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
-            Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
-            Double_t tpz = pt * lam; //track z coordinate of momentum
-            
-            Double_t mass = p->GetMass();
-            Double_t tEtot = TMath::Sqrt( tpx*tpx + tpy*tpy + tpz*tpz + mass*mass);//total energy of the track
-            
-            AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
-           if(Pass(track))//check if meets all criteria of any of our cuts
-                         //if it does not delete it and take next good track
-             { 
-               delete track;
-               delete part;
-               continue;
-             }
-            particles->AddParticle(totalNevents,part);//put track and particle on the run
-            tracks->AddParticle(totalNevents,track);
-
-          }
-         tarray->Clear(); //clear the array
-         
-        /**************************************/
-       /**************************************/
-      /**************************************/  
-     totalNevents++;
-    }
-  
-    //save environment (resouces) --
-    //clean your place after the work
-    CloseFiles(aTracksFile,aClustersFile,aGAliceFile); 
-    currentdir++;
-   }while(currentdir < Ndirs);
+  fCurrentEvent = 0;
+  return 0;
+}
+/********************************************************************/
 
+void AliHBTReaderTPC::DoOpenError( const char *va_(fmt), ...)
+{
+  // Does error display and clean-up in case error caught on Open Next Session
 
-  delete tarray;
-  fIsRead = kTRUE;
-  return 0;
- }
+   va_list ap;
+   va_start(ap,va_(fmt));
+   Error("OpenNextSession", va_(fmt), ap);
+   va_end(ap);
+   
+   delete fRunLoader;
+   fRunLoader = 0x0;
+   fTPCLoader = 0x0;
+   fCurrentDir++;
+}
 
 /********************************************************************/
-Int_t AliHBTReaderTPC::OpenFiles
-(TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
+Bool_t AliHBTReaderTPC::CheckTrack(AliTPCtrack* t) const
 {
- //opens all the files
-   
-   
-   const TString& dirname = GetDirName(event); 
-   if (dirname == "")
-    {
-      Error("OpenFiles","Can not get directory name");
-      return 4;
-    }
-   
-   TString filename = dirname +"/"+ fTrackFileName;
-   aTracksFile = TFile::Open(filename.Data());
-   if ( aTracksFile  == 0x0 ) 
-     {
-       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
-       return 1;
-     }
-   if (!aTracksFile->IsOpen())
-     {
-       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
-       return 1;
-     }
+  //Performs check of the track
+  if ( (t->GetNumberOfClusters() > fNClustMax) || (t->GetNumberOfClusters() < fNClustMin) ) return kTRUE;
+
+  Float_t chisqpercl = t->GetChi2()/((Double_t)t->GetNumberOfClusters());
+  if ( (chisqpercl < fNChi2PerClustMin) || (chisqpercl > fNChi2PerClustMax) ) return kTRUE;
   
-   filename = dirname +"/"+ fClusterFileName;
-   aClustersFile = TFile::Open(filename.Data());
-   if ( aClustersFile == 0x0 )
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 2;
-    }
-   if (!aClustersFile->IsOpen())
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 2;
-    }
+  Double_t cc[15];
+  t->GetExternalCovariance(cc);
+
+  if ( (cc[0]  < fC00Min) || (cc[0]  > fC00Max) ) return kTRUE;
+  if ( (cc[2]  < fC11Min) || (cc[2]  > fC11Max) ) return kTRUE;
+  if ( (cc[5]  < fC22Min) || (cc[5]  > fC22Max) ) return kTRUE;
+  if ( (cc[9]  < fC33Min) || (cc[9]  > fC33Max) ) return kTRUE;
+  if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
+  
+  
+  return kFALSE;
+  
+}
+/********************************************************************/
 
-   filename = dirname +"/"+ fGAliceFileName;
-   agAliceFile = TFile::Open(filename.Data());
-   if ( agAliceFile== 0x0)
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 3;
-    }
-   if (!agAliceFile->IsOpen())
-    {
-      Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
-      return 3;
-    } 
-   
-   if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice"))) 
-    {
-      Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
-      return 5;
-    }
+void AliHBTReaderTPC::SetNClustersRange(Int_t min,Int_t max)
+{
+ //sets range of Number Of Clusters that tracks have to have
+ fNClustMin = min;
+ fNClustMax = max;
+}
+/********************************************************************/
+
+void AliHBTReaderTPC::SetChi2PerCluserRange(Float_t min, Float_t max)
+{
+  //sets range of Chi2 per Cluster
+  fNChi2PerClustMin = min;
+  fNChi2PerClustMax = max;
+}
+/********************************************************************/
 
-   return 0; 
+void AliHBTReaderTPC::SetC00Range(Float_t min, Float_t max)
+{
+ //Sets range of C00 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC00Min = min;
+  fC00Max = max;
 }
 /********************************************************************/
 
+void AliHBTReaderTPC::SetC11Range(Float_t min, Float_t max)
+{
+ //Sets range of C11 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC11Min = min;
+  fC11Max = max;
+}
 /********************************************************************/
-  
-void AliHBTReaderTPC::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
+
+void AliHBTReaderTPC::SetC22Range(Float_t min, Float_t max)
 {
-  //closes the files
-  tracksFile->Close();
-  delete tracksFile;
-  tracksFile = 0x0;
-  clustersFile->Close();
-  delete clustersFile;
-  clustersFile = 0x0;
-
-  delete gAlice;
-  gAlice = 0;
-
-  if (gAliceFile) 
-   {
-     gAliceFile->Close();
-     delete gAliceFile;
-     gAliceFile = 0x0;
-   }
+ //Sets range of C22 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC22Min = min;
+  fC22Max = max;
 }
+/********************************************************************/
 
+void AliHBTReaderTPC::SetC33Range(Float_t min, Float_t max)
+{
+ //Sets range of C33 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC33Min = min;
+  fC33Max = max;
+}
 /********************************************************************/
 
+void AliHBTReaderTPC::SetC44Range(Float_t min, Float_t max)
+{
+ //Sets range of C44 parameter of covariance matrix of the track
+ //it defines uncertainty of the momentum
+  fC44Min = min;
+  fC44Max = max;
+}
+
 /********************************************************************/
 /********************************************************************/
 /********************************************************************/
 
+/*
+void (AliTPCtrack* iotrack, Double_t curv)
+{
+  Double_t x[5];
+  iotrack->GetExternalPara
+  //xk is a 
+  Double_t fP4=iotrack->GetC();
+  Double_t fP2=iotrack->GetEta();
+  
+  Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fP0, z1=fP1;
+  Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
+  Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
+  
+  fP0 += dx*(c1+c2)/(r1+r2);
+  fP1 += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
+
+}
+*/
+