Info guarded by debug flag
[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>
88378f71 18#include <TKey.h>
7be13e68 19#include <TParticle.h>
20
88378f71 21#include <AliRun.h>
bed069a4 22#include <AliRunLoader.h>
23#include <AliStack.h>
7be13e68 24#include <AliESDtrack.h>
25#include <AliESD.h>
26
27#include "AliHBTRun.h"
28#include "AliHBTEvent.h"
29#include "AliHBTParticle.h"
30#include "AliHBTParticleCut.h"
88378f71 31#include "AliHBTTrackPoints.h"
32#include "AliHBTClusterMap.h"
7be13e68 33
34ClassImp(AliHBTReaderESD)
35
bed069a4 36AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
7be13e68 37 fESDFileName(esdfilename),
bed069a4 38 fGAlFileName(galfilename),
39 fFile(0x0),
40 fRunLoader(0x0),
88378f71 41 fKeyIterator(0x0),
42 fReadParticles(kFALSE),
43 fNTrackPoints(0),
44 fdR(0.0),
45 fClusterMap(kFALSE),
46 fNTPCClustMin(0),
47 fNTPCClustMax(150),
48 fTPCChi2PerClustMin(0.0),
49 fTPCChi2PerClustMax(10e5),
50 fC00Min(0.0),
51 fC00Max(10e5),
52 fC11Min(0.0),
53 fC11Max(10e5),
54 fC22Min(0.0),
55 fC22Max(10e5),
56 fC33Min(0.0),
57 fC33Max(10e5),
58 fC44Min(0.0),
59 fC44Max(10e5),
60 fTPCC00Min(0.0),
61 fTPCC00Max(10e5),
62 fTPCC11Min(0.0),
63 fTPCC11Max(10e5),
64 fTPCC22Min(0.0),
65 fTPCC22Max(10e5),
66 fTPCC33Min(0.0),
67 fTPCC33Max(10e5),
68 fTPCC44Min(0.0),
69 fTPCC44Max(10e5)
70
7be13e68 71{
72 //cosntructor
73 if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
74 Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
75}
76/********************************************************************/
77
bed069a4 78AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
7be13e68 79 AliHBTReader(dirs),
7be13e68 80 fESDFileName(esdfilename),
bed069a4 81 fGAlFileName(galfilename),
82 fFile(0x0),
83 fRunLoader(0x0),
88378f71 84 fKeyIterator(0x0),
85 fReadParticles(kFALSE),
86 fNTrackPoints(0),
87 fdR(0.0),
88 fClusterMap(kFALSE),
89 fNTPCClustMin(0),
90 fNTPCClustMax(150),
91 fTPCChi2PerClustMin(0.0),
92 fTPCChi2PerClustMax(10e5),
93 fC00Min(0.0),
94 fC00Max(10e5),
95 fC11Min(0.0),
96 fC11Max(10e5),
97 fC22Min(0.0),
98 fC22Max(10e5),
99 fC33Min(0.0),
100 fC33Max(10e5),
101 fC44Min(0.0),
102 fC44Max(10e5),
103 fTPCC00Min(0.0),
104 fTPCC00Max(10e5),
105 fTPCC11Min(0.0),
106 fTPCC11Max(10e5),
107 fTPCC22Min(0.0),
108 fTPCC22Max(10e5),
109 fTPCC33Min(0.0),
110 fTPCC33Max(10e5),
111 fTPCC44Min(0.0),
112 fTPCC44Max(10e5)
7be13e68 113{
114 //cosntructor
115 if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
116 Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
117}
118/********************************************************************/
119
120AliHBTReaderESD::~AliHBTReaderESD()
121{
122 //desctructor
bed069a4 123 delete fRunLoader;
88378f71 124 delete fKeyIterator;
7e35966a 125 delete fFile;
7be13e68 126}
bed069a4 127/**********************************************************/
128Int_t AliHBTReaderESD::ReadNext()
7be13e68 129{
bed069a4 130//reads next event from fFile
88378f71 131 //fRunLoader is for reading Kine
bed069a4 132 //****** Tentative particle type "concentrations"
133 static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
7be13e68 134
135 Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
136 Double_t w[kNSpecies];
137 Double_t mom[3];//momentum
138 Double_t pos[3];//position
7be13e68 139
83a9f895 140 if (AliHBTParticle::GetDebug())
bed069a4 141 Info("ReadNext","Entered");
142
7be13e68 143 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
144 if (pdgdb == 0x0)
145 {
bed069a4 146 Error("Next","Can not get PDG Database Instance.");
7be13e68 147 return 1;
148 }
bed069a4 149
150 if (fParticlesEvent == 0x0) fParticlesEvent = new AliHBTEvent();
151 if (fTracksEvent == 0x0) fTracksEvent = new AliHBTEvent();
152
153 fParticlesEvent->Reset();
154 fTracksEvent->Reset();
155
7be13e68 156 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
157 {
bed069a4 158 if (fFile == 0x0)
7be13e68 159 {
bed069a4 160 fFile = OpenFile(fCurrentDir);//rl is opened here
161 if (fFile == 0x0)
7be13e68 162 {
bed069a4 163 Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
164 fCurrentDir++;
165 continue;
166 }
167 fCurrentEvent = 0;
88378f71 168 fKeyIterator = new TIter(fFile->GetListOfKeys());
169 fFile->Dump();
170 fFile->GetListOfKeys()->Print();
bed069a4 171 }
88378f71 172 TKey* key = (TKey*)fKeyIterator->Next();
173 if (key == 0x0)
174 {
175 if (AliHBTParticle::GetDebug() > 2 )
176 {
177 Info("ReadNext","No more keys.");
178 }
179 fCurrentDir++;
180 delete fKeyIterator;
181 fKeyIterator = 0x0;
182 delete fFile;//we have to assume there is no more ESD objects in the fFile
183 fFile = 0x0;
184 delete fRunLoader;
185 fRunLoader = 0x0;
186 continue;
187 }
bed069a4 188 //try to read
88378f71 189
190
191// TObject* esdobj = key->ReadObj();
192// if (esdobj == 0x0)
193// {
194// if (AliHBTParticle::GetDebug() > 2 )
195// {
196// Info("ReadNext","Key read NULL. Key Name is %s",key->GetName());
197// key->Dump();
198// }
199// continue;
200// }
201// esdobj->Dump();
202// AliESD* esd = dynamic_cast<AliESD*>(esdobj);
203
204 TString esdname = "ESD";
bed069a4 205 esdname+=fCurrentEvent;
206 AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
207 if (esd == 0x0)
208 {
88378f71 209// if (AliHBTParticle::GetDebug() > 2 )
210// {
211// Info("ReadNext","This key is not an AliESD object %s",key->GetName());
212// }
83a9f895 213 if (AliHBTParticle::GetDebug() > 2 )
bed069a4 214 {
88378f71 215 Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
bed069a4 216 }
217 fCurrentDir++;
88378f71 218 delete fKeyIterator;
219 fKeyIterator = 0x0;
bed069a4 220 delete fFile;//we have to assume there is no more ESD objects in the fFile
bed069a4 221 fFile = 0x0;
88378f71 222 delete fRunLoader;
7e35966a 223 fRunLoader = 0x0;
bed069a4 224 continue;
225 }
88378f71 226
227
228 Float_t mf = esd->GetMagneticField();
229
230 if ( (mf == 0.0) && (fNTrackPoints > 0) )
231 {
232 Error("ReadNext","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
233 fCurrentEvent++;
234 continue;
235 }
236
bed069a4 237 AliStack* stack = 0x0;
238 if (fReadParticles && fRunLoader)
239 {
240 fRunLoader->GetEvent(fCurrentEvent);
241 stack = fRunLoader->Stack();
242 }
243 Info("ReadNext","Reading Event %d",fCurrentEvent);
7be13e68 244
bed069a4 245 Int_t ntr = esd->GetNumberOfTracks();
246 Info("ReadNext","Found %d tracks.",ntr);
247 for (Int_t i = 0;i<ntr; i++)
248 {
249 AliESDtrack *esdtrack = esd->GetTrack(i);
250 if (esdtrack == 0x0)
251 {
252 Error("Next","Can not get track %d", i);
253 continue;
254 }
7e35966a 255
aafd28f6 256 //if (esdtrack->HasVertexParameters() == kFALSE)
257 if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
7e35966a 258 {
83a9f895 259 if (AliHBTParticle::GetDebug() > 2)
7e35966a 260 Info("ReadNext","Particle skipped: Data at vertex not available.");
261 continue;
262 }
263
bed069a4 264 if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
265 {
83a9f895 266 if (AliHBTParticle::GetDebug() > 2)
bed069a4 267 Info("ReadNext","Particle skipped: PID BIT is not set.");
268 continue;
269 }
7e35966a 270
bed069a4 271 esdtrack->GetESDpid(pidtable);
aafd28f6 272 //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]);
273 esdtrack->GetPxPyPz(mom);
274 //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
275 esdtrack->GetXYZ(pos);
bed069a4 276
3c69895b 277 Double_t extx;
278 Double_t extp[5];
279 esdtrack->GetExternalParameters(extx,extp);
280
281 Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
282
bed069a4 283 //Particle from kinematics
284 AliHBTParticle* particle = 0;
285 Bool_t keeppart = kFALSE;
286 if ( fReadParticles && stack )
287 {
288 if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
289 TParticle *p = stack->Particle(esdtrack->GetLabel());
290 if (p==0x0)
291 {
292 Error("ReadNext","Can not find track with such label.");
293 continue;
294 }
3c69895b 295// if(p->GetPdgCode()<0) charge = -1;
bed069a4 296 particle = new AliHBTParticle(*p,i);
3c69895b 297
bed069a4 298 }
299
300 //Here we apply Bayes' formula
301 Double_t rc=0.;
302 for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
303 if (rc==0.0)
304 {
83a9f895 305 if (AliHBTParticle::GetDebug() > 2)
bed069a4 306 Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
307 continue;
308 }
7be13e68 309
bed069a4 310 for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
7be13e68 311
83a9f895 312 if (AliHBTParticle::GetDebug() > 4)
bed069a4 313 {
314 Info("ReadNext","###########################################################################");
315 Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
316 Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
317 TString msg("Pid list got from track:");
318 for (Int_t s = 0;s<kNSpecies;s++)
319 {
320 msg+="\n ";
321 msg+=s;
322 msg+="(";
3c69895b 323 msg+=charge*GetSpeciesPdgCode((ESpecies)s);
bed069a4 324 msg+="): ";
325 msg+=w[s];
326 msg+=" (";
327 msg+=pidtable[s];
328 msg+=")";
329 }
330 Info("ReadNext","%s",msg.Data());
83a9f895 331 }//if (AliHBTParticle::GetDebug()>4)
88378f71 332
333 AliHBTTrackPoints* tpts = 0x0;
334 if (fNTrackPoints > 0)
335 {
336 tpts = new AliHBTTrackPoints(fNTrackPoints,esdtrack,mf,fdR);
337 }
338
339 AliHBTClusterMap* cmap = 0x0;
340 if ( fClusterMap )
341 {
342 cmap = new AliHBTClusterMap(esdtrack);
343 }
7be13e68 344
bed069a4 345 for (Int_t s = 0; s<kNSpecies; s++)
346 {
dd82cadc 347 Float_t pp = w[s];
348 if (pp == 0.0) continue;
7be13e68 349
3c69895b 350 Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
bed069a4 351 if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
7be13e68 352
bed069a4 353 Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
354 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 355
bed069a4 356 AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
357 mom[0], mom[1], mom[2], tEtot,
358 pos[0], pos[1], pos[2], 0.);
359 //copy probabilitis of other species (if not zero)
360 for (Int_t k = 0; k<kNSpecies; k++)
361 {
362 if (k == s) continue;
363 if (w[k] == 0.0) continue;
3c69895b 364 track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
bed069a4 365 }
88378f71 366
bed069a4 367 if(Pass(track))//check if meets all criteria of any of our cuts
368 //if it does not delete it and take next good track
369 {
370 delete track;
371 continue;
372 }
88378f71 373
374 //Single Particle cuts on cluster map and track points rather do not have sense
375 if (tpts)
376 {
377 track->SetTrackPoints(tpts);
378 }
379
380 if (cmap)
381 {
382 track->SetClusterMap(cmap);
383 }
bed069a4 384
385 fTracksEvent->AddParticle(track);
386 if (particle) fParticlesEvent->AddParticle(particle);
387 keeppart = kTRUE;
388
83a9f895 389 if (AliHBTParticle::GetDebug() > 4 )
bed069a4 390 {
391 Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
392 track->Print();
393 }
394 }//for (Int_t s = 0; s<kNSpecies; s++)
395
88378f71 396 if (keeppart == kFALSE)
397 {
398 delete particle;//particle was not stored in event
399 delete tpts;
400 delete cmap;
401 }
bed069a4 402
403 }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
404
405 Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
dd82cadc 406 fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
bed069a4 407 fNEventsRead,fCurrentEvent,fCurrentDir);
7be13e68 408
bed069a4 409 fCurrentEvent++;
410 fNEventsRead++;
411 delete esd;
412 return 0;//success -> read one event
413 }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
414
415 return 1; //no more directories to read
416}
417/**********************************************************/
418void AliHBTReaderESD::Rewind()
419{
7e35966a 420 //rewinds reading
88378f71 421 delete fKeyIterator;
bed069a4 422 delete fFile;
88378f71 423 fFile = 0x0;
424 fKeyIterator = 0x0;
bed069a4 425 delete fRunLoader;
88378f71 426 fRunLoader = 0x0;
bed069a4 427 fCurrentDir = 0;
428 fNEventsRead = 0;
429 fCurrentEvent++;
bed069a4 430}
7be13e68 431/**********************************************************/
bed069a4 432
7be13e68 433TFile* AliHBTReaderESD::OpenFile(Int_t n)
434{
bed069a4 435//opens fFile with kine tree
7be13e68 436
437 const TString& dirname = GetDirName(n);
438 if (dirname == "")
439 {
440 Error("OpenFiles","Can not get directory name");
441 return 0x0;
442 }
443 TString filename = dirname +"/"+ fESDFileName;
444 TFile *ret = TFile::Open(filename.Data());
445
446 if ( ret == 0x0)
447 {
bed069a4 448 Error("OpenFiles","Can't open fFile %s",filename.Data());
7be13e68 449 return 0x0;
450 }
451 if (!ret->IsOpen())
452 {
bed069a4 453 Error("OpenFiles","Can't open fFile %s",filename.Data());
7be13e68 454 return 0x0;
455 }
456
88378f71 457 if (fReadParticles )
bed069a4 458 {
459 fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
460 if (fRunLoader == 0x0)
461 {
462 Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
463 delete ret;
464 return 0x0;
465 }
88378f71 466
bed069a4 467 fRunLoader->LoadHeader();
468 if (fRunLoader->LoadKinematics())
469 {
470 Error("Next","Error occured while loading kinematics.");
471 delete fRunLoader;
472 delete ret;
473 return 0x0;
474 }
475 }
88378f71 476
7be13e68 477 return ret;
478}
479/**********************************************************/
480
481Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
482{
483 //returns pdg code from the PID index
484 //ask jura about charge
485 switch (spec)
486 {
487 case kESDElectron:
488 return kPositron;
489 break;
490 case kESDMuon:
491 return kMuonPlus;
492 break;
493 case kESDPion:
494 return kPiPlus;
495 break;
496 case kESDKaon:
497 return kKPlus;
498 break;
499 case kESDProton:
500 return kProton;
501 break;
502 default:
503 ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
504 break;
505 }
506 return 0;
507}
88378f71 508/********************************************************************/
509
510void AliHBTReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
511{
512 //sets range of Number Of Clusters that tracks have to have
513 fNTPCClustMin = min;
514 fNTPCClustMax = max;
515}
516/********************************************************************/
517
518void AliHBTReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
519{
520 //sets range of Chi2 per Cluster
521 fTPCChi2PerClustMin = min;
522 fTPCChi2PerClustMax = max;
523}
524/********************************************************************/
525
526void AliHBTReaderESD::SetC00Range(Float_t min, Float_t max)
527{
528 //Sets range of C00 parameter of covariance matrix of the track
529 //it defines uncertainty of the momentum
530 fC00Min = min;
531 fC00Max = max;
532}
533/********************************************************************/
534
535void AliHBTReaderESD::SetC11Range(Float_t min, Float_t max)
536{
537 //Sets range of C11 parameter of covariance matrix of the track
538 //it defines uncertainty of the momentum
539 fC11Min = min;
540 fC11Max = max;
541}
542/********************************************************************/
543
544void AliHBTReaderESD::SetC22Range(Float_t min, Float_t max)
545{
546 //Sets range of C22 parameter of covariance matrix of the track
547 //it defines uncertainty of the momentum
548 fC22Min = min;
549 fC22Max = max;
550}
551/********************************************************************/
552
553void AliHBTReaderESD::SetC33Range(Float_t min, Float_t max)
554{
555 //Sets range of C33 parameter of covariance matrix of the track
556 //it defines uncertainty of the momentum
557 fC33Min = min;
558 fC33Max = max;
559}
560/********************************************************************/
561
562void AliHBTReaderESD::SetC44Range(Float_t min, Float_t max)
563{
564 //Sets range of C44 parameter of covariance matrix of the track
565 //it defines uncertainty of the momentum
566 fC44Min = min;
567 fC44Max = max;
568}