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