Dummy copy constructor
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderESD.cxx
CommitLineData
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
18ClassImp(AliHBTReaderESD)
19
20AliHBTReaderESD::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
32AliHBTReaderESD::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
45AliHBTReaderESD::~AliHBTReaderESD()
46{
47 //desctructor
48 delete fParticles;
49 delete fTracks;
50}
51/********************************************************************/
52
53AliHBTEvent* 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
66AliHBTEvent* 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
79Int_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/********************************************************************/
92Int_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
106Int_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] + mass*mass);//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
278TFile* 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
306Int_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}