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

#include "AliHBTReaderInternal.h"

#include <iostream.h>
//#include <fstream.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>

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

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


AliHBTReaderInternal nnn;

ClassImp(AliHBTReaderInternal)
/********************************************************************/

AliHBTReaderInternal::AliHBTReaderInternal()
{
  fParticles = 0x0; 
  fTracks = 0x0;
  fIsRead = kFALSE;
}
/********************************************************************/

AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):fFileName(filename)
{ 
 fParticles = new AliHBTRun();
 fTracks    = new AliHBTRun();
 fIsRead = kFALSE;
}
/********************************************************************/
AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
  AliHBTReader(dirs),fFileName(filename)
{ 
 fParticles = new AliHBTRun();
 fTracks    = new AliHBTRun();
 fIsRead = kFALSE;
}
AliHBTReaderInternal::~AliHBTReaderInternal()
 {
 //desctructor
   delete fParticles;
   delete fTracks;
 }
/********************************************************************/
AliHBTEvent* AliHBTReaderInternal::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* AliHBTReaderInternal::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 AliHBTReaderInternal::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 AliHBTReaderInternal::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 AliHBTReaderInternal::Read(AliHBTRun* particles, AliHBTRun *tracks)
 {
 //reads data and puts put to the particles and tracks objects
 //reurns 0 if everything is OK
 //
  cout<<"AliHBTReaderInternal::Read()"<<endl;
  Int_t i; //iterator and some temprary values
  Int_t Nevents = 0;
  TFile *aFile;//file with tracks
  AliHBTParticle* p = 0x0;
  
  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();
 
  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
   }

  TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
  TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);

  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
   {
    
    if( (i=OpenFile(aFile, currentdir)) )
     {
       Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
       return i;
     }
   /***************************/
   /***************************/
   /***************************/
    
     TTree* tree = (TTree*)aFile->Get("data");
     if (tree == 0x0)
      {
       Error("Read","Can not get the tree");
       return 1;
      }
     
     TBranch *trackbranch=tree->GetBranch("tracks");//get the branch with tracks
     if (trackbranch == 0x0) ////check if we got the branch
       {//if not return with error
         Warning("Read","Can't find a branch with tracks !\n"); 
       }
     else
      {
        trackbranch->SetAddress(&tbuffer);
      }

     TBranch *partbranch=tree->GetBranch("particles");//get the branch with particles
     if (partbranch == 0x0) ////check if we got the branch
       {//if not return with error
         Warning("Read","Can't find a branch with particles !\n"); 
       }
     else
      {
        partbranch->SetAddress(&pbuffer);
      }
     
     Nevents = (Int_t)tree->GetEntries();
     cout<<"________________________________________________________\n";
     cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
     
     for (Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
       {
         tree->GetEvent(currentEvent);
         if (partbranch)
          {
            for(i = 0; i < pbuffer->GetEntries(); i++)
             {
               p = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
               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)) continue; //check if we are intersted with particles of this type 
                                                   //if not take next partilce
               AliHBTParticle* part = new AliHBTParticle(*p);
               particles->AddParticle(currentEvent,part);//put track and particle on the run
             }
            cout<<"Read: "<<particles->GetNumberOfParticlesInEvent(currentEvent)<<" particles  ";
          }
         else cout<<"Read: 0 particles  ";
         
         if (trackbranch)
          {
            for(i = 0; i < tbuffer->GetEntries(); i++)
             {
               p = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
               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)) continue; //check if we are intersted with particles of this type 
                                                   //if not take next partilce
               AliHBTParticle* part = new AliHBTParticle(*p);
               tracks->AddParticle(currentEvent,part);//put track and particle on the run
             }
            cout<<tracks->GetNumberOfParticlesInEvent(currentEvent)<<" tracks"<<endl;
          }
         else cout<<" 0 tracks"<<endl;

         
       }
    

   /***************************/
   /***************************/
   /***************************/
   currentdir++;
   aFile->Close();
   aFile = 0x0;
   }while(currentdir < Ndirs);

  fIsRead = kTRUE;
  return 0;
}

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

Int_t AliHBTReaderInternal::OpenFile(TFile*& aFile,Int_t event)
{

   const TString& dirname = GetDirName(event); 
   if (dirname == "")
    {
      Error("OpenFile","Can not get directory name");
      return 4;
    }
   
   TString filename = dirname +"/"+ fFileName;
   aFile = TFile::Open(filename.Data());
   if ( aFile  == 0x0 ) 
     {
       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
       return 1;
     }
   if (!aFile->IsOpen())
     {
       Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
       return 1;
     }
   return 0; 
}


Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile)
 {
  //reads tracks from runs and writes them to file
    Int_t i,j;
  
  
  TFile *histoOutput = TFile::Open(outfile,"recreate");
  
  if (!histoOutput->IsOpen())
   {
     cout<<"File is not opened"<<endl;
     return 1;
   }
    

  TTree *tracktree = new TTree("data","Tree with tracks");

  TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
  TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
  tbuffer->SetOwner();
  pbuffer->SetOwner();

  TClonesArray &particles = *pbuffer;
  TClonesArray &tracks = *tbuffer;
    
  TString name("Tracks");
 
  
  Int_t NT = reader->GetNumberOfTrackEvents();
  Int_t NP = reader->GetNumberOfPartEvents();
  
  Bool_t trck = (NT > 0) ? kTRUE : kFALSE;
  Bool_t part = (NP > 0) ? kTRUE : kFALSE;

  TBranch *trackbranch = 0x0, *partbranch = 0x0;
  
  
  if (trck) trackbranch = tracktree->Branch("tracks","TClonesArray",&tbuffer);
  if (part) partbranch = tracktree->Branch("particles","TClonesArray",&pbuffer);


  if ( (trck) && (part) && (NP != NT))
   {
     cerr<<"Warning number of track and particle events is different"<<endl;
   }
  
  Int_t N;
  if (NT >= NP ) N = NT; else N = NP;

  for ( i =0;i< N; i++)
    {
      if (trck && (i<=NT))
       {
         AliHBTEvent* trackev = reader->GetTrackEvent(i);
         for ( j = 0; j< trackev->GetNumberOfParticles();j++)
          {
            cout<<j<<"\r";
	    new (tracks[j]) AliHBTParticle(*(trackev->GetParticle(j)));
          }

       }
      cout<<endl;
      
      if (part && (i<=NP))
       {
        AliHBTEvent* partev = reader->GetParticleEvent(i);
        for ( j = 0; j< partev->GetNumberOfParticles();j++)
         {
           cout<<j<<"\r";
           new (particles[j]) AliHBTParticle(*(partev->GetParticle(j)));
         }
       
       }

      histoOutput->cd();
      tracktree->Fill();
      tbuffer->Delete();
      pbuffer->Delete();

     }

  
  histoOutput->cd();
  tracktree->Write(0,TObject::kOverwrite);
  histoOutput->Close();

  return 0;
 }
Commit	Line	Data
91923a38	1	#include "AliHBTReaderInternal.h"
	2
	3	#include <iostream.h>
	4	//#include <fstream.h>
	5	#include <TTree.h>
	6	#include <TFile.h>
	7	#include <TParticle.h>
	8
	9	#include <AliRun.h>
	10	#include <AliMagF.h>
	11
	12	#include "AliHBTRun.h"
	13	#include "AliHBTEvent.h"
	14	#include "AliHBTParticle.h"
	15	#include "AliHBTParticleCut.h"
	16
	17
	18	AliHBTReaderInternal nnn;
	19
	20	ClassImp(AliHBTReaderInternal)
	21	/********************************************************************/
	22
	23	AliHBTReaderInternal::AliHBTReaderInternal()
	24	{
	25	fParticles = 0x0;
	26	fTracks = 0x0;
	27	fIsRead = kFALSE;
	28	}
	29	/********************************************************************/
	30
	31	AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):fFileName(filename)
	32	{
	33	fParticles = new AliHBTRun();
	34	fTracks = new AliHBTRun();
	35	fIsRead = kFALSE;
	36	}
	37	/********************************************************************/
	38	AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
	39	AliHBTReader(dirs),fFileName(filename)
	40	{
	41	fParticles = new AliHBTRun();
	42	fTracks = new AliHBTRun();
	43	fIsRead = kFALSE;
	44	}
	45	AliHBTReaderInternal::~AliHBTReaderInternal()
	46	{
	47	//desctructor
	48	delete fParticles;
	49	delete fTracks;
	50	}
	51	/********************************************************************/
	52	AliHBTEvent* AliHBTReaderInternal::GetParticleEvent(Int_t n)
	53	{
	54	//returns Nth event with simulated particles
	55	if (!fIsRead)
	56	if(Read(fParticles,fTracks))
	57	{
	58	Error("GetParticleEvent","Error in reading");
	59	return 0x0;
	60	}
	61	return fParticles->GetEvent(n);
	62	}
	63	/********************************************************************/
	64	AliHBTEvent* AliHBTReaderInternal::GetTrackEvent(Int_t n)
65	{
66	//returns Nth event with reconstructed tracks
67	if (!fIsRead)
68	if(Read(fParticles,fTracks))
69	{
70	Error("GetTrackEvent","Error in reading");
71	return 0x0;
72	}
73	return fTracks->GetEvent(n);
74	}
75	/********************************************************************/
76
77	Int_t AliHBTReaderInternal::GetNumberOfPartEvents()
78	{
79	//returns number of events of particles
80	if (!fIsRead)
81	if ( Read(fParticles,fTracks))
82	{
83	Error("GetNumberOfPartEvents","Error in reading");
84	return 0;
85	}
86	return fParticles->GetNumberOfEvents();
87	}
88
89	/********************************************************************/
90	Int_t AliHBTReaderInternal::GetNumberOfTrackEvents()
91	{
92	//returns number of events of tracks
93	if (!fIsRead)
94	if(Read(fParticles,fTracks))
95	{
96	Error("GetNumberOfTrackEvents","Error in reading");
97	return 0;
98	}
99	return fTracks->GetNumberOfEvents();
100	}
101	/********************************************************************/
102
103
104	Int_t AliHBTReaderInternal::Read(AliHBTRun* particles, AliHBTRun *tracks)
105	{
106	//reads data and puts put to the particles and tracks objects
107	//reurns 0 if everything is OK
108	//
109	cout<<"AliHBTReaderInternal::Read()"<<endl;
110	Int_t i; //iterator and some temprary values
111	Int_t Nevents = 0;
112	TFile *aFile;//file with tracks
113	AliHBTParticle* p = 0x0;
114
115	if (!particles) //check if an object is instatiated
116	{
117	Error("Read"," particles object must instatiated before passing it to the reader");
118	}
119	if (!tracks) //check if an object is instatiated
120	{
121	Error("Read"," tracks object must instatiated before passing it to the reader");
122	}
123	particles->Reset();//clear runs == delete all old events
124	tracks->Reset();
125
126	Int_t currentdir = 0;
127
128	Int_t Ndirs;
129	if (fDirs) //if array with directories is supplied by user
130	{
131	Ndirs = fDirs->GetEntries(); //get the number if directories
132	}
133	else
134	{
135	Ndirs = 0; //if the array is not supplied read only from current directory
136	}
137
138	TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
139	TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
140
141	do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
142	{
143
144	if( (i=OpenFile(aFile, currentdir)) )
145	{
146	Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
147	return i;
148	}
149	/***************************/
150	/***************************/
151	/***************************/
152
153	TTree* tree = (TTree*)aFile->Get("data");
154	if (tree == 0x0)
155	{
156	Error("Read","Can not get the tree");
157	return 1;
158	}
159
160	TBranch *trackbranch=tree->GetBranch("tracks");//get the branch with tracks
161	if (trackbranch == 0x0) ////check if we got the branch
162	{//if not return with error
163	Warning("Read","Can't find a branch with tracks !\n");
164	}
165	else
166	{
167	trackbranch->SetAddress(&tbuffer);
168	}
169
170	TBranch *partbranch=tree->GetBranch("particles");//get the branch with particles
171	if (partbranch == 0x0) ////check if we got the branch
172	{//if not return with error
173	Warning("Read","Can't find a branch with particles !\n");
174	}
175	else
176	{
177	partbranch->SetAddress(&pbuffer);
178	}
179
180	Nevents = (Int_t)tree->GetEntries();
181	cout<<"________________________________________________________\n";
182	cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
183
184	for (Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
185	{
186	tree->GetEvent(currentEvent);
187	if (partbranch)
188	{
189	for(i = 0; i < pbuffer->GetEntries(); i++)
190	{
191	p = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
192	if(p == 0x0) continue; //if returned pointer is NULL
193	if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
194	if(Pass(p)) continue; //check if we are intersted with particles of this type
195	//if not take next partilce
196	AliHBTParticle* part = new AliHBTParticle(*p);
197	particles->AddParticle(currentEvent,part);//put track and particle on the run
198	}
199	cout<<"Read: "<<particles->GetNumberOfParticlesInEvent(currentEvent)<<" particles ";
200	}
201	else cout<<"Read: 0 particles ";
202
203	if (trackbranch)
204	{
205	for(i = 0; i < tbuffer->GetEntries(); i++)
206	{
207	p = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
208	if(p == 0x0) continue; //if returned pointer is NULL
209	if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
210	if(Pass(p)) continue; //check if we are intersted with particles of this type
211	//if not take next partilce
212	AliHBTParticle* part = new AliHBTParticle(*p);
213	tracks->AddParticle(currentEvent,part);//put track and particle on the run
214	}
215	cout<<tracks->GetNumberOfParticlesInEvent(currentEvent)<<" tracks"<<endl;
216	}
217	else cout<<" 0 tracks"<<endl;
218
219
220
221	}
222
223
224	/***************************/
225	/***************************/
226	/***************************/
227	currentdir++;
228	aFile->Close();
229	aFile = 0x0;
230	}while(currentdir < Ndirs);
231
232	fIsRead = kTRUE;
233	return 0;
234	}
235
236	/********************************************************************/
237
238	Int_t AliHBTReaderInternal::OpenFile(TFile*& aFile,Int_t event)
239	{
240
241	const TString& dirname = GetDirName(event);
242	if (dirname == "")
243	{
244	Error("OpenFile","Can not get directory name");
245	return 4;
246	}
247
248	TString filename = dirname +"/"+ fFileName;
249	aFile = TFile::Open(filename.Data());
250	if ( aFile == 0x0 )
251	{
252	Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
253	return 1;
254	}
255	if (!aFile->IsOpen())
256	{
257	Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
258	return 1;
259	}
260	return 0;
261	}
262
263
264
265
266	Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile)
267	{
268	//reads tracks from runs and writes them to file
269	Int_t i,j;
270
271
272	TFile *histoOutput = TFile::Open(outfile,"recreate");
273
274	if (!histoOutput->IsOpen())
275	{
276	cout<<"File is not opened"<<endl;
277	return 1;
278	}
279
280
281	TTree *tracktree = new TTree("data","Tree with tracks");
282
283	TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
284	TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
285	tbuffer->SetOwner();
286	pbuffer->SetOwner();
287
288	TClonesArray &particles = *pbuffer;
289	TClonesArray &tracks = *tbuffer;
290
291	TString name("Tracks");
292
293
294	Int_t NT = reader->GetNumberOfTrackEvents();
295	Int_t NP = reader->GetNumberOfPartEvents();
296
297	Bool_t trck = (NT > 0) ? kTRUE : kFALSE;
298	Bool_t part = (NP > 0) ? kTRUE : kFALSE;
299
300	TBranch trackbranch = 0x0, partbranch = 0x0;
301
302
303	if (trck) trackbranch = tracktree->Branch("tracks","TClonesArray",&tbuffer);
304	if (part) partbranch = tracktree->Branch("particles","TClonesArray",&pbuffer);
305
306
307
308	if ( (trck) && (part) && (NP != NT))
309	{
310	cerr<<"Warning number of track and particle events is different"<<endl;
311	}
312
313	Int_t N;
314	if (NT >= NP ) N = NT; else N = NP;
315
316	for ( i =0;i< N; i++)
317	{
318	if (trck && (i<=NT))
319	{
320	AliHBTEvent* trackev = reader->GetTrackEvent(i);
321	for ( j = 0; j< trackev->GetNumberOfParticles();j++)
322	{
323	cout<<j<<"\r";
324	new (tracks[j]) AliHBTParticle(*(trackev->GetParticle(j)));
325	}
326
327	}
328	cout<<endl;
329
330	if (part && (i<=NP))
331	{
91923a38	332	AliHBTEvent* partev = reader->GetParticleEvent(i);
91923a38	333	for ( j = 0; j< partev->GetNumberOfParticles();j++)
	334	{
	335	cout<<j<<"\r";
	336	new (particles[j]) AliHBTParticle(*(partev->GetParticle(j)));
	337	}
91923a38	338
	339	}
	340
91923a38	341	histoOutput->cd();
91923a38	342	tracktree->Fill();
91923a38	343	tbuffer->Delete();
	344	pbuffer->Delete();
	345
	346	}
	347
	348
	349
	350	histoOutput->cd();
91923a38	351	tracktree->Write(0,TObject::kOverwrite);
91923a38	352	histoOutput->Close();
	353
	354	return 0;
	355	}