[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderTPC.cxx

#include "AliHBTReaderTPC.h"

#include <Riostream.h>
//#include <Riostream.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>

#include <AliRun.h>
#include <AliMagF.h>
#include <AliTPCtrack.h>
#include <AliTPCParam.h>
#include <AliTPCtracker.h>

#include "AliHBTRun.h"
#include "AliHBTEvent.h"
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.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)
{
  //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, const Char_t* clusterfilename,
                  const Char_t* galicefilename):
                  AliHBTReader(dirs), fTrackFileName(trackfilename),
                  fClusterFileName(clusterfilename),fGAliceFileName(galicefilename)

{
  //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()
 {
 //desctructor
   delete fParticles;
   delete fTracks;
 }
/********************************************************************/

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);
 }
/********************************************************************/
AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
 {
 //returns Nth event with reconstructed tracks
   if (!fIsRead) 
    if(Read(fParticles,fTracks))
     {
       Error("GetTrackEvent","Error in reading");
       return 0x0;
     }
   return fTracks->GetEvent(n);
 }
/********************************************************************/

Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
 {
 //returns number of events of particles
   if (!fIsRead) 
    if ( Read(fParticles,fTracks))
     {
       Error("GetNumberOfPartEvents","Error in reading");
       return 0;
     }
   return fParticles->GetNumberOfEvents();
 }

/********************************************************************/
Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
 {
 //returns number of events of tracks
  if (!fIsRead)
    if(Read(fParticles,fTracks))
     {
       Error("GetNumberOfTrackEvents","Error in reading");
       return 0;
     }
  return fTracks->GetNumberOfEvents();
 }
/********************************************************************/


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
   {
     Error("Read"," particles object must instatiated before passing it to the reader");
   }
  if (!tracks)  //check if an object is instatiated
   {
     Error("Read"," tracks object must instatiated before passing it to the reader");
   }
  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
   {
     Ndirs = fDirs->GetEntries(); //get the number if directories
   }
  else
   {
     Ndirs = 0; //if the array is not supplied read only from current directory
   }

  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
   {
    
    if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
     {
       Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
       currentdir++;
       continue;
     }
  
    
    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: B="<<gAlice->Field()->SolenoidField()<<"T"<<endl;
      AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
     }
    else
     {//if not return an error
       Error("Read","Can not find Header tree (TreeE) in gAlice");
       currentdir++;
       continue;
     }
  
    aClustersFile->cd();//set cluster file active 
    AliTPCParam *TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60_150x60");
    if (!TPCParam) 
      { 
       TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60");
       if (!TPCParam) 
        { 
          Error("Read","TPC parameters have not been found !\n");
          currentdir++;
          continue;
        }
      }

  
    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"); 
            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("Read","Can't get a branch with TPC tracks !\n"); 
            continue;
          }
         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"); 
            continue;
          }
         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(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);
            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);


  delete tarray;
  fIsRead = kTRUE;
  return 0;
 }

/********************************************************************/
Int_t AliHBTReaderTPC::OpenFiles
(TFile*& aTracksFile, TFile*& aClustersFile, TFile*& agAliceFile,Int_t event)
{
 //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;
     }
  
   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;
    }

   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;
    }

   return 0; 
}
/********************************************************************/

/********************************************************************/
  
void AliHBTReaderTPC::CloseFiles(TFile*& tracksFile, TFile*& clustersFile, TFile*& gAliceFile)
{
  //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;
   }
}

/********************************************************************/

/********************************************************************/
/********************************************************************/
/********************************************************************/
Commit	Line	Data
1b446896	1	#include "AliHBTReaderTPC.h"
1b446896	2
d0c23b58	3	#include <Riostream.h>
d0c23b58	4	//#include <Riostream.h>
1b446896	5	#include <TTree.h>
1b446896	6	#include <TFile.h>
16701d1b	7	#include <TParticle.h>
1b446896	8
16701d1b	9	#include <AliRun.h>
16701d1b	10	#include <AliMagF.h>
1b446896	11	#include <AliTPCtrack.h>
	12	#include <AliTPCParam.h>
	13	#include <AliTPCtracker.h>
	14
1b446896	15	#include "AliHBTRun.h"
	16	#include "AliHBTEvent.h"
	17	#include "AliHBTParticle.h"
	18	#include "AliHBTParticleCut.h"
	19
	20
	21	ClassImp(AliHBTReaderTPC)
	22	//reader for TPC tracking
	23	//needs galice.root, AliTPCtracks.root, AliTPCclusters.root, good_tracks_tpc
	24	//
	25	//more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html
	26	//Piotr.Skowronski@cern.ch
	27
	28	AliHBTReaderTPC::
	29	AliHBTReaderTPC(const Char_t* trackfilename,const Char_t* clusterfilename,
16701d1b	30	const Char_t* galicefilename):
1b446896	31	fTrackFileName(trackfilename),fClusterFileName(clusterfilename),
16701d1b	32	fGAliceFileName(galicefilename)
1b446896	33	{
16701d1b	34	//constructor,
1b446896	35	//Defaults:
	36	// trackfilename = "AliTPCtracks.root"
	37	// clusterfilename = "AliTPCclusters.root"
1b446896	38	// galicefilename = "" - this means: Do not open gAlice file -
	39	// just leave the global pointer untached
	40
	41	fParticles = new AliHBTRun();
	42	fTracks = new AliHBTRun();
16701d1b	43	fIsRead = kFALSE;
	44	}
	45	/********************************************************************/
	46	AliHBTReaderTPC::
	47	AliHBTReaderTPC(TObjArray* dirs,
df9a1768	48	const Char_t* trackfilename, const Char_t* clusterfilename,
	49	const Char_t* galicefilename):
	50	AliHBTReader(dirs), fTrackFileName(trackfilename),
	51	fClusterFileName(clusterfilename),fGAliceFileName(galicefilename)
1b446896	52
16701d1b	53	{
	54	//constructor,
	55	//Defaults:
	56	// trackfilename = "AliTPCtracks.root"
	57	// clusterfilename = "AliTPCclusters.root"
	58	// galicefilename = "" - this means: Do not open gAlice file -
	59	// just leave the global pointer untached
	60
16701d1b	61	fParticles = new AliHBTRun();
16701d1b	62	fTracks = new AliHBTRun();
1b446896	63	fIsRead = kFALSE;
16701d1b	64
1b446896	65	}
	66	/********************************************************************/
	67
	68	AliHBTReaderTPC::~AliHBTReaderTPC()
	69	{
	70	//desctructor
	71	delete fParticles;
	72	delete fTracks;
	73	}
	74	/********************************************************************/
	75
	76	AliHBTEvent* AliHBTReaderTPC::GetParticleEvent(Int_t n)
	77	{
	78	//returns Nth event with simulated particles
16701d1b	79	if (!fIsRead)
	80	if(Read(fParticles,fTracks))
	81	{
	82	Error("GetParticleEvent","Error in reading");
	83	return 0x0;
	84	}
1b446896	85	return fParticles->GetEvent(n);
	86	}
	87	/********************************************************************/
	88	AliHBTEvent* AliHBTReaderTPC::GetTrackEvent(Int_t n)
	89	{
	90	//returns Nth event with reconstructed tracks
16701d1b	91	if (!fIsRead)
	92	if(Read(fParticles,fTracks))
	93	{
	94	Error("GetTrackEvent","Error in reading");
	95	return 0x0;
	96	}
1b446896	97	return fTracks->GetEvent(n);
	98	}
	99	/********************************************************************/
	100
	101	Int_t AliHBTReaderTPC::GetNumberOfPartEvents()
	102	{
	103	//returns number of events of particles
16701d1b	104	if (!fIsRead)
	105	if ( Read(fParticles,fTracks))
	106	{
	107	Error("GetNumberOfPartEvents","Error in reading");
	108	return 0;
	109	}
1b446896	110	return fParticles->GetNumberOfEvents();
	111	}
	112
	113	/********************************************************************/
	114	Int_t AliHBTReaderTPC::GetNumberOfTrackEvents()
	115	{
	116	//returns number of events of tracks
16701d1b	117	if (!fIsRead)
	118	if(Read(fParticles,fTracks))
	119	{
	120	Error("GetNumberOfTrackEvents","Error in reading");
	121	return 0;
	122	}
1b446896	123	return fTracks->GetNumberOfEvents();
	124	}
	125	/********************************************************************/
	126
	127
	128	Int_t AliHBTReaderTPC::Read(AliHBTRun* particles, AliHBTRun *tracks)
	129	{
	130	//reads data and puts put to the particles and tracks objects
	131	//reurns 0 if everything is OK
	132	//
b71a5edb	133	cout<<"AliHBTReaderTPC::Read()"<<endl;
1b446896	134	Int_t i; //iterator and some temprary values
16701d1b	135	Int_t Nevents = 0;
	136	Int_t totalNevents = 0;
	137	TFile *aTracksFile;//file with tracks
	138	TFile *aClustersFile;//file with clusters
	139	TFile *aGAliceFile;//!ile name with galice
	140
1b446896	141	if (!particles) //check if an object is instatiated
1b446896	142	{
7da9bdd2	143	Error("Read"," particles object must instatiated before passing it to the reader");
1b446896	144	}
	145	if (!tracks) //check if an object is instatiated
	146	{
7da9bdd2	147	Error("Read"," tracks object must instatiated before passing it to the reader");
1b446896	148	}
	149	particles->Reset();//clear runs == delete all old events
	150	tracks->Reset();
16701d1b	151
	152	TObjArray *tarray = new TObjArray(5000); //cotainer for tpc tracks
	153	tarray->SetOwner(); //set the ownership of the objects it contains
	154	//when array is is deleted or cleared all objects
	155	//that it contains are deleted
	156	Int_t currentdir = 0;
b71a5edb	157
	158	Int_t Ndirs;
	159	if (fDirs) //if array with directories is supplied by user
	160	{
	161	Ndirs = fDirs->GetEntries(); //get the number if directories
	162	}
	163	else
	164	{
	165	Ndirs = 0; //if the array is not supplied read only from current directory
	166	}
	167
16701d1b	168	do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
1b446896	169	{
16701d1b	170
	171	if( (i=OpenFiles(aTracksFile,aClustersFile,aGAliceFile,currentdir)) )
	172	{
7da9bdd2	173	Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
5ee8b891	174	currentdir++;
5ee8b891	175	continue;
16701d1b	176	}
1b446896	177
1b446896	178
16701d1b	179	if (gAlice->TreeE())//check if tree E exists
	180	{
	181	Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
	182	cout<<"________________________________________________________\n";
	183	cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
5ee8b891	184	cout<<"Setting Magnetic Field: B="<<gAlice->Field()->SolenoidField()<<"T"<<endl;
5ee8b891	185	AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
16701d1b	186	}
	187	else
	188	{//if not return an error
7da9bdd2	189	Error("Read","Can not find Header tree (TreeE) in gAlice");
5ee8b891	190	currentdir++;
5ee8b891	191	continue;
16701d1b	192	}
1b446896	193
16701d1b	194	aClustersFile->cd();//set cluster file active
2f1287d2	195	AliTPCParam TPCParam= (AliTPCParam)aClustersFile->Get("75x40_100x60_150x60");
16701d1b	196	if (!TPCParam)
16701d1b	197	{
2f1287d2	198	TPCParam= (AliTPCParam*)aClustersFile->Get("75x40_100x60");
	199	if (!TPCParam)
	200	{
	201	Error("Read","TPC parameters have not been found !\n");
	202	currentdir++;
	203	continue;
	204	}
16701d1b	205	}
1b446896	206
1b446896	207
16701d1b	208	for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
	209	{
	210	cout<<"Reading Event "<<currentEvent<<endl;
	211	/**************************************/
	212	/**************************************/
	213	/**************************************/
1b446896	214
16701d1b	215	aTracksFile->cd();//set track file active
16701d1b	216
1b446896	217	Char_t treename[100];
	218	sprintf(treename,"TreeT_TPC_%d",currentEvent);//prepare name of the tree
	219
	220	TTree *tracktree=0;
	221
16701d1b	222	tracktree=(TTree*)aTracksFile->Get(treename);//get the tree
1b446896	223	if (!tracktree) //check if we got the tree
1b446896	224	{//if not return with error
7da9bdd2	225	Error("Read","Can't get a tree with TPC tracks !\n");
5ee8b891	226	continue;
1b446896	227	}
	228
	229	TBranch *trackbranch=tracktree->GetBranch("tracks");//get the branch with tracks
	230	if (!trackbranch) ////check if we got the branch
	231	{//if not return with error
7da9bdd2	232	Error("Read","Can't get a branch with TPC tracks !\n");
5ee8b891	233	continue;
1b446896	234	}
	235	Int_t NTPCtracks=(Int_t)tracktree->GetEntries();//get number of TPC tracks
	236	cout<<"Found "<<NTPCtracks<<" TPC tracks.\n";
	237	//Copy tracks to array
	238
	239	AliTPCtrack *iotrack=0;
	240
16701d1b	241	aClustersFile->cd();//set cluster file active
1b446896	242	AliTPCtracker *tracker = new AliTPCtracker(TPCParam,currentEvent);//create the tacker for this event
	243	if (!tracker) //check if it has created succeffuly
	244	{//if not return with error
7da9bdd2	245	Error("Read","Can't get a tracker !\n");
5ee8b891	246	continue;
1b446896	247	}
	248	tracker->LoadInnerSectors();
	249	tracker->LoadOuterSectors();
	250
	251	for (i=0; i<NTPCtracks; i++) //loop over all tpc tracks
	252	{
	253	iotrack=new AliTPCtrack; //create new tracks
	254	trackbranch->SetAddress(&iotrack); //tell the branch ehere to put track data from tree(file)
	255	tracktree->GetEvent(i); //stream track i to the iotrack
	256	tracker->CookLabel(iotrack,0.1); //calculate (cook) the label of the tpc track
	257	//which is the label of corresponding simulated particle
	258	tarray->AddLast(iotrack); //put the track in the array
	259	}
	260
16701d1b	261	aTracksFile->Delete(treename);//delete tree from memmory (and leave untached on disk)- we do not need it any more
16701d1b	262	aTracksFile->Delete("tracks");//delete branch from memmory
1b446896	263	delete tracker; //delete tracker
	264
	265	tracker = 0x0;
	266	trackbranch = 0x0;
	267	tracktree = 0x0;
1b446896	268
	269	Double_t xk;
	270	Double_t par[5];
	271	Float_t phi, lam, pt;//angles and transverse momentum
	272	Int_t label; //label of the current track
16701d1b	273
	274	aGAliceFile->cd();
	275	gAlice->GetEvent(currentEvent);
	276
	277	gAlice->Particles();
	278
	279	for (i=0; i<NTPCtracks; i++) //loop over all good tracks
1b446896	280	{
16701d1b	281	iotrack = (AliTPCtrack*)tarray->At(i);
	282	label = iotrack->GetLabel();
	283
	284	if (label < 0) continue;
1b446896	285
16701d1b	286	TParticle p = (TParticle)gAlice->Particle(label);
1b446896	287
5ee8b891	288	if(p == 0x0) continue; //if returned pointer is NULL
5ee8b891	289	if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
f1e2da22	290
16701d1b	291	if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
16701d1b	292	//if not take next partilce
1b446896	293
16701d1b	294	AliHBTParticle* part = new AliHBTParticle(*p);
1b446896	295	if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
	296	//if it does not delete it and take next good track
	297
16701d1b	298	iotrack->PropagateToVertex();
7da9bdd2	299
1b446896	300	iotrack->GetExternalParameters(xk,par); //get properties of the track
	301	phi=TMath::ASin(par[2]) + iotrack->GetAlpha();
	302	if (phi<-TMath::Pi()) phi+=2*TMath::Pi();
	303	if (phi>=TMath::Pi()) phi-=2*TMath::Pi();
	304	lam=par[3];
	305	pt=1.0/TMath::Abs(par[4]);
	306
	307	Double_t tpx = pt * TMath::Cos(phi); //track x coordinate of momentum
	308	Double_t tpy = pt * TMath::Sin(phi); //track y coordinate of momentum
	309	Double_t tpz = pt * lam; //track z coordinate of momentum
	310
16701d1b	311	Double_t mass = p->GetMass();
1b446896	312	Double_t tEtot = TMath::Sqrt( tpxtpx + tpytpy + tpztpz + massmass);//total energy of the track
1b446896	313
16701d1b	314	AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), tpx, tpy , tpz, tEtot, 0., 0., 0., 0.);
5ee8b891	315	if(Pass(track))//check if meets all criteria of any of our cuts
7da9bdd2	316	//if it does not delete it and take next good track
	317	{
	318	delete track;
	319	delete part;
	320	continue;
	321	}
16701d1b	322	particles->AddParticle(totalNevents,part);//put track and particle on the run
16701d1b	323	tracks->AddParticle(totalNevents,track);
1b446896	324
	325	}
	326	tarray->Clear(); //clear the array
16701d1b	327
16701d1b	328	/**************************************/
1b446896	329	/**************************************/
16701d1b	330	/**************************************/
	331	totalNevents++;
	332	}
1b446896	333
16701d1b	334	//save environment (resouces) --
	335	//clean your place after the work
	336	CloseFiles(aTracksFile,aClustersFile,aGAliceFile);
	337	currentdir++;
b71a5edb	338	}while(currentdir < Ndirs);
16701d1b	339
16701d1b	340
1b446896	341	delete tarray;
1b446896	342	fIsRead = kTRUE;
	343	return 0;
	344	}
	345
	346	/********************************************************************/
16701d1b	347	Int_t AliHBTReaderTPC::OpenFiles
16701d1b	348	(TFile& aTracksFile, TFile& aClustersFile, TFile*& agAliceFile,Int_t event)
1b446896	349	{
1b446896	350	//opens all the files
16701d1b	351
	352
	353	const TString& dirname = GetDirName(event);
	354	if (dirname == "")
	355	{
	356	Error("OpenFiles","Can not get directory name");
	357	return 4;
	358	}
	359
	360	TString filename = dirname +"/"+ fTrackFileName;
	361	aTracksFile = TFile::Open(filename.Data());
	362	if ( aTracksFile == 0x0 )
1b446896	363	{
16701d1b	364	Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
1b446896	365	return 1;
1b446896	366	}
16701d1b	367	if (!aTracksFile->IsOpen())
	368	{
	369	Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
	370	return 1;
	371	}
	372
	373	filename = dirname +"/"+ fClusterFileName;
	374	aClustersFile = TFile::Open(filename.Data());
	375	if ( aClustersFile == 0x0 )
1b446896	376	{
16701d1b	377	Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
1b446896	378	return 2;
1b446896	379	}
16701d1b	380	if (!aClustersFile->IsOpen())
	381	{
	382	Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
	383	return 2;
	384	}
	385
	386	filename = dirname +"/"+ fGAliceFileName;
	387	agAliceFile = TFile::Open(filename.Data());
	388	if ( agAliceFile== 0x0)
	389	{
	390	Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
	391	return 3;
	392	}
	393	if (!agAliceFile->IsOpen())
	394	{
	395	Error("OpenFiles","Can't open file with TPC clusters named %s",filename.Data());
	396	return 3;
	397	}
	398
	399	if (!(gAlice=(AliRun*)agAliceFile->Get("gAlice")))
	400	{
	401	Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
	402	return 5;
	403	}
1b446896	404
16701d1b	405	return 0;
1b446896	406	}
16701d1b	407	/********************************************************************/
1b446896	408
	409	/********************************************************************/
	410
16701d1b	411	void AliHBTReaderTPC::CloseFiles(TFile& tracksFile, TFile& clustersFile, TFile*& gAliceFile)
1b446896	412	{
1b446896	413	//closes the files
16701d1b	414	tracksFile->Close();
7da9bdd2	415	delete tracksFile;
16701d1b	416	tracksFile = 0x0;
16701d1b	417	clustersFile->Close();
7da9bdd2	418	delete clustersFile;
16701d1b	419	clustersFile = 0x0;
7da9bdd2	420
	421	delete gAlice;
	422	gAlice = 0;
	423
	424	if (gAliceFile)
	425	{
	426	gAliceFile->Close();
	427	delete gAliceFile;
	428	gAliceFile = 0x0;
	429	}
1b446896	430	}
	431
	432	/********************************************************************/
	433
	434	/********************************************************************/
	435	/********************************************************************/
	436	/********************************************************************/
	437