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