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

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