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