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


#include "AliHBTReaderITSv1.h"
#include "AliHBTEvent.h"
#include "AliHBTRun.h"
#include "AliHBTParticle.h"
#include "AliHBTParticleCut.h"

#include <Riostream.h>

#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TBranch.h>
#include <TObjArray.h>
#include <TParticle.h>
#include <TString.h>
#include <TObjString.h>

#include <AliRun.h>
#include <AliMagF.h>
#include <AliKalmanTrack.h>
#include <AliITSIOTrack.h>

ClassImp(AliHBTReaderITSv1)
/********************************************************************/

AliHBTReaderITSv1::
AliHBTReaderITSv1(const Char_t* tracksfilename,const Char_t* galicefilename):
                 fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
 {
     fParticles = new AliHBTRun();
     fTracks    = new AliHBTRun();
     fIsRead = kFALSE;
 }
/********************************************************************/

AliHBTReaderITSv1::
AliHBTReaderITSv1(TObjArray* dirs, const Char_t* tracksfilename,const Char_t* galicefilename):
                 AliHBTReader(dirs),
                 fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
 {
   fParticles = new AliHBTRun();
   fTracks    = new AliHBTRun();
   fIsRead    = kFALSE;
 }
/********************************************************************/

AliHBTReaderITSv1::~AliHBTReaderITSv1()
{
   delete fParticles;
   delete fTracks;
}
/********************************************************************/

AliHBTEvent* AliHBTReaderITSv1::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* AliHBTReaderITSv1::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 AliHBTReaderITSv1::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 AliHBTReaderITSv1::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 AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
{
 cout<<"AliHBTReaderITSv1::Read()"<<endl;
 Int_t Nevents = 0;
 AliITSIOTrack *iotrack=new AliITSIOTrack;
 Int_t currentdir = 0;
 Int_t Ndirs;
 Int_t totalNevents = 0;
 
 if (fDirs)
  {
    Ndirs = fDirs->GetEntries();
  }
 else
  {
    Ndirs = 0;
  }
 
 do //do while is good even if 
  {  
   TFile* gAliceFile = OpenGAliceFile(currentdir);
   if(gAliceFile == 0x0)
    {
       Error("Read","Can not open the file with gAlice");
       delete iotrack;
       return 1;
    }
   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");
       delete iotrack;
       return 4;
     }

   TFile *file = OpenTrackFile(currentdir);
   if(file == 0x0)
    {
       Error("Read","Can not open the file with ITS tracks V1");
       delete iotrack;
       return 2;
    }
    
   Int_t naccepted = 0;
   char tname[30];
   
   for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
    { 
      cout<<"Reading Event "<<currentEvent;
      
      sprintf(tname,"TreeT%d",currentEvent);
      file->cd(); 
      TTree *tracktree=(TTree*)file->Get(tname);
      TBranch *tbranch=tracktree->GetBranch("ITStracks");
      tbranch->SetAddress(&iotrack);
      
      gAliceFile->cd();
      gAlice->GetEvent(currentEvent);
      gAlice->Particles();

      Int_t nentr=(Int_t)tracktree->GetEntries();
      
      cout<<".  Found "<<nentr<<" tracks.";
      fflush(0);
      
      for (Int_t i=0; i<nentr; i++) 
       {

        tracktree->GetEvent(i);
        if(!iotrack) continue;       
        Int_t label = iotrack->GetLabel();
        if (label < 0) 
         {
           continue;
         }

        TParticle *p = (TParticle*)gAlice->Particle(label);
        if(!p)
         {
           Warning("Read","Can not get particle with label &d",label);
           continue;
         }
        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
        
        Double_t px=iotrack->GetPx();
        Double_t py=iotrack->GetPy();
        Double_t pz=iotrack->GetPz();
        Double_t mass = p->GetMass();
        Double_t tEtot = TMath::Sqrt(px*px + py*py + pz*pz + mass*mass);//total energy of the track
 
        Double_t x= iotrack->GetX();
        Double_t y= iotrack->GetY();
        Double_t z= iotrack->GetZ();
        
        AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, px, py , pz, tEtot, x, y, z, 0.);
        if(Pass(track)) { delete  track;continue;}//check if meets all criteria of any of our cuts
                                                  //if it does not delete it and take next good track

        particles->AddParticle(totalNevents,part);//put track and particle on the run
        tracks->AddParticle(totalNevents,track);
        naccepted++;
       }//end loop over tracks in the event

       totalNevents++;
       cout<<"  Accepted "<<naccepted<<" tracks"<<endl;
     }//end of loop over events in current directory
    
    gAliceFile->Close();
    delete gAliceFile;
    gAliceFile = 0;
    
    file->Close(); 
    delete file;
    file = 0;
    currentdir++;
   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array


  delete iotrack;
  fIsRead = kTRUE;
  return 0;
 
 }
/********************************************************************/

TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
{
//opens files to be read for given directoru nomber in fDirs Array
   const TString& dirname = GetDirName(ndir); 
   if (dirname == "")
    {
      Error("OpenGAliceFile","Can not get directory name");
      return 0x0;
    }
   TString filename = dirname + "/" + fITSTracksFileName;

   TFile *file = TFile::Open(filename.Data());   
   if (!file)
    {
      Error("Read","Can not open file %s",filename.Data());
      return 0x0;
    }
   if (!file->IsOpen())
    {
      Error("Read","Can not open file %s",filename.Data());
      return 0x0;
    }
   
   return file;
}


/********************************************************************/
TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
{
  const TString& dirname = GetDirName(ndir); 
   if (dirname == "")
    {
      Error("OpenGAliceFile","Can not get directory name");
      return 0x0;
    }
  
  TString filename = dirname + "/" + fGAliceFileName;

  TFile* gAliceFile = TFile::Open(filename.Data());
  if ( gAliceFile== 0x0)
   {
     Error("OpenFiles","Can't open file named %s",filename.Data());
     return 0x0;
   }
  if (!gAliceFile->IsOpen())
   {
     Error("OpenFiles","Can't open file named %s",filename.Data());
     return 0x0;
   }

  if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
   {
     Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
     gAliceFile->Close();
     delete gAliceFile;
     return 0x0;
   }
  
  return gAliceFile;
}

/********************************************************************/
/********************************************************************/
Commit	Line	Data
	1
	2	#include "AliHBTReaderITSv1.h"
	3	#include "AliHBTEvent.h"
	4	#include "AliHBTRun.h"
	5	#include "AliHBTParticle.h"
	6	#include "AliHBTParticleCut.h"
	7
	8	#include <Riostream.h>
	9
	10	#include <TROOT.h>
	11	#include <TFile.h>
	12	#include <TTree.h>
	13	#include <TBranch.h>
	14	#include <TObjArray.h>
	15	#include <TParticle.h>
	16	#include <TString.h>
	17	#include <TObjString.h>
	18
	19	#include <AliRun.h>
	20	#include <AliMagF.h>
	21	#include <AliKalmanTrack.h>
	22	#include <AliITSIOTrack.h>
	23
	24	ClassImp(AliHBTReaderITSv1)
	25	/********************************************************************/
	26
	27	AliHBTReaderITSv1::
	28	AliHBTReaderITSv1(const Char_t* tracksfilename,const Char_t* galicefilename):
	29	fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
	30	{
	31	fParticles = new AliHBTRun();
	32	fTracks = new AliHBTRun();
	33	fIsRead = kFALSE;
	34	}
	35	/********************************************************************/
	36
	37	AliHBTReaderITSv1::
	38	AliHBTReaderITSv1(TObjArray* dirs, const Char_t* tracksfilename,const Char_t* galicefilename):
	39	AliHBTReader(dirs),
	40	fITSTracksFileName(tracksfilename),fGAliceFileName(galicefilename)
	41	{
	42	fParticles = new AliHBTRun();
	43	fTracks = new AliHBTRun();
	44	fIsRead = kFALSE;
	45	}
	46	/********************************************************************/
	47
	48	AliHBTReaderITSv1::~AliHBTReaderITSv1()
	49	{
	50	delete fParticles;
	51	delete fTracks;
	52	}
	53	/********************************************************************/
	54
	55	AliHBTEvent* AliHBTReaderITSv1::GetParticleEvent(Int_t n)
	56	{
	57	//returns Nth event with simulated particles
	58	if (!fIsRead)
	59	if(Read(fParticles,fTracks))
	60	{
	61	Error("GetParticleEvent","Error in reading");
	62	return 0x0;
	63	}
	64
	65	return fParticles->GetEvent(n);
	66	}
	67	/********************************************************************/
	68
	69	AliHBTEvent* AliHBTReaderITSv1::GetTrackEvent(Int_t n)
	70	{
	71	//returns Nth event with reconstructed tracks
	72	if (!fIsRead)
	73	if(Read(fParticles,fTracks))
	74	{
	75	Error("GetTrackEvent","Error in reading");
	76	return 0x0;
	77	}
	78	return fTracks->GetEvent(n);
	79	}
	80	/********************************************************************/
	81
	82	Int_t AliHBTReaderITSv1::GetNumberOfPartEvents()
	83	{
	84	//returns number of events of particles
	85	if (!fIsRead)
	86	if(Read(fParticles,fTracks))
	87	{
	88	Error("GetNumberOfPartEvents","Error in reading");
	89	return 0;
	90	}
	91	return fParticles->GetNumberOfEvents();
	92	}
	93
	94	/********************************************************************/
	95	Int_t AliHBTReaderITSv1::GetNumberOfTrackEvents()
	96	{
	97	//returns number of events of tracks
	98	if (!fIsRead)
	99	if(Read(fParticles,fTracks))
	100	{
	101	Error("GetNumberOfTrackEvents","Error in reading");
	102	return 0;
	103	}
	104	return fTracks->GetNumberOfEvents();
	105	}
	106	/********************************************************************/
	107
	108	Int_t AliHBTReaderITSv1::Read(AliHBTRun* particles, AliHBTRun *tracks)
	109	{
	110	cout<<"AliHBTReaderITSv1::Read()"<<endl;
	111	Int_t Nevents = 0;
	112	AliITSIOTrack *iotrack=new AliITSIOTrack;
	113	Int_t currentdir = 0;
	114	Int_t Ndirs;
	115	Int_t totalNevents = 0;
	116
	117	if (fDirs)
	118	{
	119	Ndirs = fDirs->GetEntries();
	120	}
	121	else
	122	{
	123	Ndirs = 0;
	124	}
	125
	126	do //do while is good even if
	127	{
	128	TFile* gAliceFile = OpenGAliceFile(currentdir);
	129	if(gAliceFile == 0x0)
	130	{
	131	Error("Read","Can not open the file with gAlice");
	132	delete iotrack;
	133	return 1;
	134	}
	135	if (gAlice->TreeE())//check if tree E exists
	136	{
	137	Nevents = (Int_t)gAlice->TreeE()->GetEntries();//if yes get number of events in gAlice
	138	cout<<"________________________________________________________\n";
	139	cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
	140	cout<<"Setting Magnetic Field. Factor is "<<gAlice->Field()->Factor()<<endl;
	141	AliKalmanTrack::SetConvConst(100/0.299792458/0.2/gAlice->Field()->Factor());
	142	}
	143	else
	144	{//if not return an error
	145	Error("Read","Can not find Header tree (TreeE) in gAlice");
	146	delete iotrack;
	147	return 4;
	148	}
	149
	150	TFile *file = OpenTrackFile(currentdir);
	151	if(file == 0x0)
	152	{
	153	Error("Read","Can not open the file with ITS tracks V1");
	154	delete iotrack;
	155	return 2;
	156	}
	157
	158	Int_t naccepted = 0;
	159	char tname[30];
	160
	161	for (Int_t currentEvent = 0; currentEvent < Nevents; currentEvent++)
	162	{
	163	cout<<"Reading Event "<<currentEvent;
	164
	165	sprintf(tname,"TreeT%d",currentEvent);
	166	file->cd();
	167	TTree tracktree=(TTree)file->Get(tname);
	168	TBranch *tbranch=tracktree->GetBranch("ITStracks");
	169	tbranch->SetAddress(&iotrack);
	170
	171	gAliceFile->cd();
	172	gAlice->GetEvent(currentEvent);
	173	gAlice->Particles();
	174
	175	Int_t nentr=(Int_t)tracktree->GetEntries();
	176
	177	cout<<". Found "<<nentr<<" tracks.";
	178	fflush(0);
	179
	180	for (Int_t i=0; i<nentr; i++)
	181	{
	182
	183	tracktree->GetEvent(i);
	184	if(!iotrack) continue;
	185	Int_t label = iotrack->GetLabel();
	186	if (label < 0)
	187	{
	188	continue;
	189	}
	190
	191	TParticle p = (TParticle)gAlice->Particle(label);
	192	if(!p)
	193	{
	194	Warning("Read","Can not get particle with label &d",label);
	195	continue;
	196	}
	197	if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
	198	//if not take next partilce
	199
	200	AliHBTParticle* part = new AliHBTParticle(*p,i);
	201	if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
	202	//if it does not delete it and take next good track
	203
	204	Double_t px=iotrack->GetPx();
	205	Double_t py=iotrack->GetPy();
	206	Double_t pz=iotrack->GetPz();
	207	Double_t mass = p->GetMass();
	208	Double_t tEtot = TMath::Sqrt(pxpx + pypy + pzpz + massmass);//total energy of the track
	209
	210	Double_t x= iotrack->GetX();
	211	Double_t y= iotrack->GetY();
	212	Double_t z= iotrack->GetZ();
	213
	214	AliHBTParticle* track = new AliHBTParticle(p->GetPdgCode(), i, px, py , pz, tEtot, x, y, z, 0.);
	215	if(Pass(track)) { delete track;continue;}//check if meets all criteria of any of our cuts
	216	//if it does not delete it and take next good track
	217
	218	particles->AddParticle(totalNevents,part);//put track and particle on the run
	219	tracks->AddParticle(totalNevents,track);
	220	naccepted++;
	221	}//end loop over tracks in the event
	222
	223	totalNevents++;
	224	cout<<" Accepted "<<naccepted<<" tracks"<<endl;
	225	}//end of loop over events in current directory
	226
	227	gAliceFile->Close();
	228	delete gAliceFile;
	229	gAliceFile = 0;
	230
	231	file->Close();
	232	delete file;
	233	file = 0;
	234	currentdir++;
	235	}while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
	236
	237
	238	delete iotrack;
	239	fIsRead = kTRUE;
	240	return 0;
	241
	242	}
	243	/********************************************************************/
	244
	245	TFile* AliHBTReaderITSv1::OpenTrackFile(Int_t ndir)
	246	{
	247	//opens files to be read for given directoru nomber in fDirs Array
	248	const TString& dirname = GetDirName(ndir);
	249	if (dirname == "")
	250	{
	251	Error("OpenGAliceFile","Can not get directory name");
	252	return 0x0;
	253	}
	254	TString filename = dirname + "/" + fITSTracksFileName;
	255
	256	TFile *file = TFile::Open(filename.Data());
	257	if (!file)
	258	{
	259	Error("Read","Can not open file %s",filename.Data());
	260	return 0x0;
	261	}
	262	if (!file->IsOpen())
	263	{
	264	Error("Read","Can not open file %s",filename.Data());
	265	return 0x0;
	266	}
	267
	268	return file;
	269	}
	270
	271
	272	/********************************************************************/
	273	TFile* AliHBTReaderITSv1::OpenGAliceFile(Int_t ndir)
	274	{
	275	const TString& dirname = GetDirName(ndir);
	276	if (dirname == "")
	277	{
	278	Error("OpenGAliceFile","Can not get directory name");
	279	return 0x0;
	280	}
	281
	282	TString filename = dirname + "/" + fGAliceFileName;
	283
	284	TFile* gAliceFile = TFile::Open(filename.Data());
	285	if ( gAliceFile== 0x0)
	286	{
	287	Error("OpenFiles","Can't open file named %s",filename.Data());
	288	return 0x0;
	289	}
	290	if (!gAliceFile->IsOpen())
	291	{
	292	Error("OpenFiles","Can't open file named %s",filename.Data());
	293	return 0x0;
	294	}
	295
	296	if (!(gAlice=(AliRun*)gAliceFile->Get("gAlice")))
	297	{
	298	Error("OpenFiles","gAlice have not been found on %s !\n",filename.Data());
	299	gAliceFile->Close();
	300	delete gAliceFile;
	301	return 0x0;
	302	}
	303
	304	return gAliceFile;
	305	}
	306
	307	/********************************************************************/
	308	/********************************************************************/
	309
	310