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

//____________________________________________________________________
//////////////////////////////////////////////////////////////////////
//                                                                  //
//  class AliHBTReaderITSv1                                         //
//                                                                  //
//  Reader for ITSv1 tracks. Not maintained since v1 is not         //
//  supposed to be used                                             //
//                                                                  //
//////////////////////////////////////////////////////////////////////

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