Inctrementing class versions
[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
bed069a4 10#include <AliRunLoader.h>
11#include <AliStack.h>
7be13e68 12#include <AliESDtrack.h>
13#include <AliESD.h>
14
15#include "AliHBTRun.h"
16#include "AliHBTEvent.h"
17#include "AliHBTParticle.h"
18#include "AliHBTParticleCut.h"
19
bed069a4 20AliHBTReaderESD r;
21
7be13e68 22ClassImp(AliHBTReaderESD)
23
bed069a4 24AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
7be13e68 25 fESDFileName(esdfilename),
bed069a4 26 fGAlFileName(galfilename),
27 fFile(0x0),
28 fRunLoader(0x0),
29 fReadParticles(kFALSE)
7be13e68 30{
31 //cosntructor
32 if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
33 Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
34}
35/********************************************************************/
36
bed069a4 37AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
7be13e68 38 AliHBTReader(dirs),
7be13e68 39 fESDFileName(esdfilename),
bed069a4 40 fGAlFileName(galfilename),
41 fFile(0x0),
42 fRunLoader(0x0),
43 fReadParticles(kFALSE)
7be13e68 44{
45 //cosntructor
46 if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
47 Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
48}
49/********************************************************************/
50
51AliHBTReaderESD::~AliHBTReaderESD()
52{
53 //desctructor
bed069a4 54 delete fFile;
55 delete fRunLoader;
7be13e68 56}
bed069a4 57/**********************************************************/
58Int_t AliHBTReaderESD::ReadNext()
7be13e68 59{
bed069a4 60//reads next event from fFile
61
62 //****** Tentative particle type "concentrations"
63 static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
7be13e68 64
65 Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
66 Double_t w[kNSpecies];
67 Double_t mom[3];//momentum
68 Double_t pos[3];//position
7be13e68 69
bed069a4 70 if (AliHBTParticle::fgDebug)
71 Info("ReadNext","Entered");
72
7be13e68 73 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
74 if (pdgdb == 0x0)
75 {
bed069a4 76 Error("Next","Can not get PDG Database Instance.");
7be13e68 77 return 1;
78 }
bed069a4 79
80 if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
81 if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
82
83 fParticlesEvent->Reset();
84 fTracksEvent->Reset();
85
7be13e68 86 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
87 {
bed069a4 88 if (fFile == 0x0)
7be13e68 89 {
bed069a4 90 fFile = OpenFile(fCurrentDir);//rl is opened here
91 if (fFile == 0x0)
7be13e68 92 {
bed069a4 93 Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
94 fCurrentDir++;
95 continue;
96 }
97 fCurrentEvent = 0;
98 }
99
100 //try to read
101 TString esdname;
102 esdname+=fCurrentEvent;
103 AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
104 if (esd == 0x0)
105 {
106 if (AliHBTParticle::fgDebug > 2 )
107 {
108 Info("ReadNext","Can not find ESD object named %s.",esdname.Data());
109 }
110 fCurrentDir++;
111 delete fFile;//we have to assume there is no more ESD objects in the fFile
112 delete fRunLoader;
113 fFile = 0x0;
114 continue;
115 }
116 AliStack* stack = 0x0;
117 if (fReadParticles && fRunLoader)
118 {
119 fRunLoader->GetEvent(fCurrentEvent);
120 stack = fRunLoader->Stack();
121 }
122 Info("ReadNext","Reading Event %d",fCurrentEvent);
7be13e68 123
bed069a4 124 Int_t ntr = esd->GetNumberOfTracks();
125 Info("ReadNext","Found %d tracks.",ntr);
126 for (Int_t i = 0;i<ntr; i++)
127 {
128 AliESDtrack *esdtrack = esd->GetTrack(i);
129 if (esdtrack == 0x0)
130 {
131 Error("Next","Can not get track %d", i);
132 continue;
133 }
134 if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
135 {
136 if (AliHBTParticle::fgDebug > 2)
137 Info("ReadNext","Particle skipped: PID BIT is not set.");
138 continue;
139 }
7be13e68 140
bed069a4 141 esdtrack->GetESDpid(pidtable);
142 esdtrack->GetPxPyPz(mom);
143 esdtrack->GetXYZ(pos);
144
145 //Particle from kinematics
146 AliHBTParticle* particle = 0;
147 Bool_t keeppart = kFALSE;
148 if ( fReadParticles && stack )
149 {
150 if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
151 TParticle *p = stack->Particle(esdtrack->GetLabel());
152 if (p==0x0)
153 {
154 Error("ReadNext","Can not find track with such label.");
155 continue;
156 }
157 particle = new AliHBTParticle(*p,i);
158 }
159
160 //Here we apply Bayes' formula
161 Double_t rc=0.;
162 for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
163 if (rc==0.0)
164 {
165 if (AliHBTParticle::fgDebug > 2)
166 Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
167 continue;
168 }
7be13e68 169
bed069a4 170 for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
7be13e68 171
bed069a4 172 if (AliHBTParticle::fgDebug > 4)
173 {
174 Info("ReadNext","###########################################################################");
175 Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
176 Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
177 TString msg("Pid list got from track:");
178 for (Int_t s = 0;s<kNSpecies;s++)
179 {
180 msg+="\n ";
181 msg+=s;
182 msg+="(";
183 msg+=GetSpeciesPdgCode((ESpecies)s);
184 msg+="): ";
185 msg+=w[s];
186 msg+=" (";
187 msg+=pidtable[s];
188 msg+=")";
189 }
190 Info("ReadNext","%s",msg.Data());
191 }//if (AliHBTParticle::fgDebug>4)
7be13e68 192
bed069a4 193 for (Int_t s = 0; s<kNSpecies; s++)
194 {
195 if (w[s] == 0.0) continue;
7be13e68 196
bed069a4 197 Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
198 if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
7be13e68 199
bed069a4 200 Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
201 Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
7be13e68 202
bed069a4 203 AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
204 mom[0], mom[1], mom[2], tEtot,
205 pos[0], pos[1], pos[2], 0.);
206 //copy probabilitis of other species (if not zero)
207 for (Int_t k = 0; k<kNSpecies; k++)
208 {
209 if (k == s) continue;
210 if (w[k] == 0.0) continue;
211 track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
212 }
213 if(Pass(track))//check if meets all criteria of any of our cuts
214 //if it does not delete it and take next good track
215 {
216 delete track;
217 continue;
218 }
219
220 fTracksEvent->AddParticle(track);
221 if (particle) fParticlesEvent->AddParticle(particle);
222 keeppart = kTRUE;
223
224 if (AliHBTParticle::fgDebug > 4 )
225 {
226 Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
227 track->Print();
228 }
229 }//for (Int_t s = 0; s<kNSpecies; s++)
230
231 if (keeppart == kFALSE) delete particle;//particle was not stored in event
232
233 }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
234
235 Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
236 fParticlesEvent->GetNumberOfParticles(), fTracksEvent->GetNumberOfParticles(),
237 fNEventsRead,fCurrentEvent,fCurrentDir);
7be13e68 238
bed069a4 239 fCurrentEvent++;
240 fNEventsRead++;
241 delete esd;
242 return 0;//success -> read one event
243 }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
244
245 return 1; //no more directories to read
246}
247/**********************************************************/
248void AliHBTReaderESD::Rewind()
249{
250 delete fFile;
251 fFile = 0;
252 delete fRunLoader;
253 fCurrentDir = 0;
254 fNEventsRead = 0;
255 fCurrentEvent++;
7be13e68 256
bed069a4 257}
7be13e68 258/**********************************************************/
bed069a4 259
7be13e68 260TFile* AliHBTReaderESD::OpenFile(Int_t n)
261{
bed069a4 262//opens fFile with kine tree
7be13e68 263
264 const TString& dirname = GetDirName(n);
265 if (dirname == "")
266 {
267 Error("OpenFiles","Can not get directory name");
268 return 0x0;
269 }
270 TString filename = dirname +"/"+ fESDFileName;
271 TFile *ret = TFile::Open(filename.Data());
272
273 if ( ret == 0x0)
274 {
bed069a4 275 Error("OpenFiles","Can't open fFile %s",filename.Data());
7be13e68 276 return 0x0;
277 }
278 if (!ret->IsOpen())
279 {
bed069a4 280 Error("OpenFiles","Can't open fFile %s",filename.Data());
7be13e68 281 return 0x0;
282 }
283
bed069a4 284 if (fReadParticles)
285 {
286 fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
287 if (fRunLoader == 0x0)
288 {
289 Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
290 delete ret;
291 return 0x0;
292 }
293 fRunLoader->LoadHeader();
294 if (fRunLoader->LoadKinematics())
295 {
296 Error("Next","Error occured while loading kinematics.");
297 delete fRunLoader;
298 delete ret;
299 return 0x0;
300 }
301 }
7be13e68 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}