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

#include "AliHBTReaderESD.h"

#include <TPDGCode.h>
#include <TString.h>
#include <TObjString.h>
#include <TTree.h>
#include <TFile.h>
#include <TParticle.h>

#include <AliESDtrack.h>
#include <AliESD.h>

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

ClassImp(AliHBTReaderESD)

AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
 fParticles(new AliHBTRun()),
 fTracks(new AliHBTRun()),
 fESDFileName(esdfilename),
 fIsRead(kFALSE)
{
  //cosntructor
  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
    Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/
  
AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
 AliHBTReader(dirs), 
 fParticles(new AliHBTRun()),
 fTracks(new AliHBTRun()),
 fESDFileName(esdfilename),
 fIsRead(kFALSE)
{
  //cosntructor
  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
    Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
}
/********************************************************************/

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

AliHBTEvent* AliHBTReaderESD::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* AliHBTReaderESD::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 AliHBTReaderESD::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 AliHBTReaderESD::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 AliHBTReaderESD::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 totalNevents = 0;
  Int_t currentdir = 0; //number of current directory name is fDirs array
  
  Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
  Double_t w[kNSpecies];
  Double_t mom[3];//momentum
  Double_t pos[3];//position
   //****** Tentative particle type "concentrations"
  const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
  
  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
  if (pdgdb == 0x0)
   {
     Error("Read","Can not get PDG Database Instance.");
     return 1;
   }
  if (!particles) //check if an object is instatiated
   {
     Error("Read"," particles object must instatiated before passing it to the reader");
     return 1;
   }
  if (!tracks)  //check if an object is instatiated
   {
     Error("Read"," tracks object must instatiated before passing it to the reader");
     return 1;
   }
  particles->Reset();//clear runs == delete all old events
  tracks->Reset();

  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 "./"
   {
     TFile* file = OpenFile(currentdir);
     if (file == 0x0)
       {
         Error("Read","Cannot get File for dir no. %d",currentdir);
         currentdir++;
         continue;
       } 
       
     for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events
      {
       TString esdname;
       esdname+=currentEvent;
       AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname));
       if (esd == 0x0)
         {
           if (AliHBTParticle::fgDebug > 2 )
            {
              Info("Read","Can not find ESD object named %s.",esdname.Data());
            }
           currentdir++;
           break;//we have to assume there is no more ESD objects in the file
         } 
       
       Info("Read","Reading Event %d",currentEvent);
  
       Int_t ntr = esd->GetNumberOfTracks();
       Info("Read","Found %d tracks.",ntr);
       for (Int_t i = 0;i<ntr; i++)
        {
          AliESDtrack *esdtrack = esd->GetTrack(i);
          if (esdtrack == 0x0)
           {
             Error("Read","Can not get track %d", i);
             continue;
           }
          if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
           {
             if (AliHBTParticle::fgDebug > 2) 
               Info("Read","Particle skipped: PID BIT is not set.");
             continue;
           }

          esdtrack->GetESDpid(pidtable);
          esdtrack->GetPxPyPz(mom);
          esdtrack->GetXYZ(pos);
          Double_t rc=0.;

           //Here we apply Bayes' formula
          for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
          if (rc==0.0) 
           {
             if (AliHBTParticle::fgDebug > 2) 
               Info("Read","Particle rejected since total bayessian PID probab. is zero.");
             continue;
           }

          for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;

          if (AliHBTParticle::fgDebug > 4)
           { 
             Info("Read","###########################################################################");
             Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
             Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]);
             Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc);
             TString msg("Pid list got from track:");
             for (Int_t s = 0;s<kNSpecies;s++)
              {
                msg+="\n    ";
                msg+=s;
                msg+="(";
                msg+=GetSpeciesPdgCode((ESpecies)s);
                msg+="): ";
                msg+=w[s];
                msg+=" (";
                msg+=pidtable[s];
                msg+=")";
              }
             Info("Read","%s",msg.Data());
           }//if (AliHBTParticle::fgDebug>4)
           
          for (Int_t s = 0; s<kNSpecies; s++)
           {
             if (w[s] == 0.0) continue;

             Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
             if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 

             Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
             Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track

             AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
                                                        mom[0], mom[1], mom[2], tEtot,
                                                        pos[0], pos[1], pos[2], 0.);
             //copy probabilitis of other species (if not zero)
             for (Int_t k = 0; k<kNSpecies; k++)
              {
                if (k == s) continue;
                if (w[k] == 0.0) continue;
                track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
              }
             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;
                 continue;
               }
             tracks->AddParticle(totalNevents,track);
             if (AliHBTParticle::fgDebug > 4 )
              {
                Info("Read","\n\nAdding Particle with incarnation %d",pdgcode);
                track->Print();
              }
           }//for (Int_t s = 0; s<kNSpecies; s++)

        }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
       totalNevents++;
      }//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file
      
   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
  
  fIsRead = kTRUE;
  return 0;
}      
/**********************************************************/
   
TFile* AliHBTReaderESD::OpenFile(Int_t n)
{
//opens file with kine tree

 const TString& dirname = GetDirName(n);
 if (dirname == "")
  {
   Error("OpenFiles","Can not get directory name");
   return 0x0;
  }
 TString filename = dirname +"/"+ fESDFileName;
 TFile *ret = TFile::Open(filename.Data()); 

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

Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
{
  //returns pdg code from the PID index
  //ask jura about charge
  switch (spec)
   {
     case kESDElectron:
       return kPositron;
       break;
     case kESDMuon:
       return kMuonPlus;
       break;
     case kESDPion:
       return kPiPlus;
       break;
     case kESDKaon:
       return kKPlus;
       break;
     case kESDProton:
       return kProton;
       break;
     default:
       ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
       break;
   }
  return 0;
}
Commit	Line	Data
7be13e68	1	#include "AliHBTReaderESD.h"
	2
	3	#include <TPDGCode.h>
	4	#include <TString.h>
	5	#include <TObjString.h>
	6	#include <TTree.h>
	7	#include <TFile.h>
	8	#include <TParticle.h>
	9
	10	#include <AliESDtrack.h>
	11	#include <AliESD.h>
	12
	13	#include "AliHBTRun.h"
	14	#include "AliHBTEvent.h"
	15	#include "AliHBTParticle.h"
	16	#include "AliHBTParticleCut.h"
	17
	18	ClassImp(AliHBTReaderESD)
	19
	20	AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
	21	fParticles(new AliHBTRun()),
	22	fTracks(new AliHBTRun()),
	23	fESDFileName(esdfilename),
	24	fIsRead(kFALSE)
	25	{
	26	//cosntructor
	27	if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
	28	Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
	29	}
	30	/********************************************************************/
	31
	32	AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
	33	AliHBTReader(dirs),
	34	fParticles(new AliHBTRun()),
	35	fTracks(new AliHBTRun()),
	36	fESDFileName(esdfilename),
	37	fIsRead(kFALSE)
	38	{
	39	//cosntructor
	40	if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
	41	Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
	42	}
	43	/********************************************************************/
	44
	45	AliHBTReaderESD::~AliHBTReaderESD()
	46	{
	47	//desctructor
	48	delete fParticles;
	49	delete fTracks;
	50	}
	51	/********************************************************************/
	52
	53	AliHBTEvent* AliHBTReaderESD::GetParticleEvent(Int_t n)
	54	{
	55	//returns Nth event with simulated particles
	56	if (!fIsRead)
	57	if(Read(fParticles,fTracks))
	58	{
	59	Error("GetParticleEvent","Error in reading");
	60	return 0x0;
	61	}
	62	return fParticles->GetEvent(n);
	63	}
	64	/********************************************************************/
65
66	AliHBTEvent* AliHBTReaderESD::GetTrackEvent(Int_t n)
67	{
68	//returns Nth event with reconstructed tracks
69	if (!fIsRead)
70	if(Read(fParticles,fTracks))
71	{
72	Error("GetTrackEvent","Error in reading");
73	return 0x0;
74	}
75	return fTracks->GetEvent(n);
76	}
77	/********************************************************************/
78
79	Int_t AliHBTReaderESD::GetNumberOfPartEvents()
80	{
81	//returns number of events of particles
82	if (!fIsRead)
83	if ( Read(fParticles,fTracks))
84	{
85	Error("GetNumberOfPartEvents","Error in reading");
86	return 0;
87	}
88	return fParticles->GetNumberOfEvents();
89	}
90
91	/********************************************************************/
92	Int_t AliHBTReaderESD::GetNumberOfTrackEvents()
93	{
94	//returns number of events of tracks
95	if (!fIsRead)
96	if(Read(fParticles,fTracks))
97	{
98	Error("GetNumberOfTrackEvents","Error in reading");
99	return 0;
100	}
101	return fTracks->GetNumberOfEvents();
102	}
103	/********************************************************************/
104
105
106	Int_t AliHBTReaderESD::Read(AliHBTRun* particles, AliHBTRun *tracks)
107	{
108	//reads data and puts put to the particles and tracks objects
109	//reurns 0 if everything is OK
110	//
111	Info("Read","");
112	Int_t totalNevents = 0;
113	Int_t currentdir = 0; //number of current directory name is fDirs array
114
115	Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
116	Double_t w[kNSpecies];
117	Double_t mom[3];//momentum
118	Double_t pos[3];//position
119	//****** Tentative particle type "concentrations"
120	const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
121
122	TDatabasePDG* pdgdb = TDatabasePDG::Instance();
123	if (pdgdb == 0x0)
124	{
125	Error("Read","Can not get PDG Database Instance.");
126	return 1;
127	}
128	if (!particles) //check if an object is instatiated
129	{
130	Error("Read"," particles object must instatiated before passing it to the reader");
131	return 1;
132	}
133	if (!tracks) //check if an object is instatiated
134	{
135	Error("Read"," tracks object must instatiated before passing it to the reader");
136	return 1;
137	}
138	particles->Reset();//clear runs == delete all old events
139	tracks->Reset();
140
141	Int_t Ndirs;
142	if (fDirs) //if array with directories is supplied by user
143	{
144	Ndirs = fDirs->GetEntries(); //get the number if directories
145	}
146	else
147	{
148	Ndirs = 0; //if the array is not supplied read only from current directory
149	}
150
151	do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
152	{
153	TFile* file = OpenFile(currentdir);
154	if (file == 0x0)
155	{
156	Error("Read","Cannot get File for dir no. %d",currentdir);
157	currentdir++;
158	continue;
159	}
160
161	for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events
162	{
163	TString esdname;
164	esdname+=currentEvent;
165	AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname));
166	if (esd == 0x0)
167	{
168	if (AliHBTParticle::fgDebug > 2 )
169	{
170	Info("Read","Can not find ESD object named %s.",esdname.Data());
171	}
172	currentdir++;
173	break;//we have to assume there is no more ESD objects in the file
174	}
175
176	Info("Read","Reading Event %d",currentEvent);
177
178	Int_t ntr = esd->GetNumberOfTracks();
179	Info("Read","Found %d tracks.",ntr);
180	for (Int_t i = 0;i<ntr; i++)
181	{
182	AliESDtrack *esdtrack = esd->GetTrack(i);
183	if (esdtrack == 0x0)
184	{
185	Error("Read","Can not get track %d", i);
186	continue;
187	}
188	if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
189	{
190	if (AliHBTParticle::fgDebug > 2)
191	Info("Read","Particle skipped: PID BIT is not set.");
192	continue;
193	}
194
195	esdtrack->GetESDpid(pidtable);
196	esdtrack->GetPxPyPz(mom);
197	esdtrack->GetXYZ(pos);
198	Double_t rc=0.;
199
200	//Here we apply Bayes' formula
201	for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
202	if (rc==0.0)
203	{
204	if (AliHBTParticle::fgDebug > 2)
205	Info("Read","Particle rejected since total bayessian PID probab. is zero.");
206	continue;
207	}
208
209	for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
210
211	if (AliHBTParticle::fgDebug > 4)
212	{
213	Info("Read","###########################################################################");
214	Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
215	Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]);
216	Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc);
217	TString msg("Pid list got from track:");
218	for (Int_t s = 0;s<kNSpecies;s++)
219	{
220	msg+="\n ";
221	msg+=s;
222	msg+="(";
223	msg+=GetSpeciesPdgCode((ESpecies)s);
224	msg+="): ";
225	msg+=w[s];
226	msg+=" (";
227	msg+=pidtable[s];
228	msg+=")";
229	}
230	Info("Read","%s",msg.Data());
231	}//if (AliHBTParticle::fgDebug>4)
232
233	for (Int_t s = 0; s<kNSpecies; s++)
234	{
235	if (w[s] == 0.0) continue;
236
237	Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
238	if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
239
240	Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
241	Double_t tEtot = TMath::Sqrt( mom[0]mom[0] + mom[1]mom[1] + mom[2]mom[2] + massmass);//total energy of the track
242
243	AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
244	mom[0], mom[1], mom[2], tEtot,
245	pos[0], pos[1], pos[2], 0.);
246	//copy probabilitis of other species (if not zero)
247	for (Int_t k = 0; k<kNSpecies; k++)
248	{
249	if (k == s) continue;
250	if (w[k] == 0.0) continue;
251	track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
252	}
253	if(Pass(track))//check if meets all criteria of any of our cuts
254	//if it does not delete it and take next good track
255	{
256	delete track;
257	continue;
258	}
259	tracks->AddParticle(totalNevents,track);
260	if (AliHBTParticle::fgDebug > 4 )
261	{
262	Info("Read","\n\nAdding Particle with incarnation %d",pdgcode);
263	track->Print();
264	}
265	}//for (Int_t s = 0; s<kNSpecies; s++)
266
267	}//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
268	totalNevents++;
269	}//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file
270
271	}while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
272
273	fIsRead = kTRUE;
274	return 0;
275	}
276	/**********************************************************/
277
278	TFile* AliHBTReaderESD::OpenFile(Int_t n)
279	{
280	//opens file with kine tree
281
282	const TString& dirname = GetDirName(n);
283	if (dirname == "")
284	{
285	Error("OpenFiles","Can not get directory name");
286	return 0x0;
287	}
288	TString filename = dirname +"/"+ fESDFileName;
289	TFile *ret = TFile::Open(filename.Data());
290
291	if ( ret == 0x0)
292	{
293	Error("OpenFiles","Can't open file %s",filename.Data());
294	return 0x0;
295	}
296	if (!ret->IsOpen())
297	{
298	Error("OpenFiles","Can't open file %s",filename.Data());
299	return 0x0;
300	}
301
302	return ret;
303	}
304	/**********************************************************/
305
306	Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
307	{
308	//returns pdg code from the PID index
309	//ask jura about charge
310	switch (spec)
311	{
312	case kESDElectron:
313	return kPositron;
314	break;
315	case kESDMuon:
316	return kMuonPlus;
317	break;
318	case kESDPion:
319	return kPiPlus;
320	break;
321	case kESDKaon:
322	return kKPlus;
323	break;
324	case kESDProton:
325	return kProton;
326	break;
327	default:
328	::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
329	break;
330	}
331	return 0;
332	}