1 #include "AliReaderESD.h"
2 //____________________________________________________________________
3 //////////////////////////////////////////////////////////////////////
5 // class AliReaderESD //
7 // reader for ALICE Event Summary Data (ESD). //
9 // Piotr.Skowronski@cern.ch //
11 //////////////////////////////////////////////////////////////////////
15 #include <TObjString.h>
19 #include <TParticle.h>
22 #include <TGliteXmlEventlist.h>
25 #include <AliRunLoader.h>
27 #include <AliESDtrack.h>
28 #include <AliKalmanTrack.h>
31 #include "AliAnalysis.h"
32 #include "AliAODRun.h"
34 #include "AliAODParticle.h"
35 #include "AliAODParticleCut.h"
36 #include "AliClusterMap.h"
38 ClassImp(AliReaderESD)
40 AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
41 fESDFileName(esdfilename),
42 fGAlFileName(galfilename),
47 fCheckParticlePID(kFALSE),
48 fReadMostProbableOnly(kFALSE),
52 fITSTrackPoints(kFALSE),
53 fITSTrackPointsType(AliTrackPoints::kITS),
55 fReadCentralBarrel(kTRUE),
60 fTPCChi2PerClustMin(0.0),
61 fTPCChi2PerClustMax(10e5),
87 if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
88 Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
90 /********************************************************************/
92 AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
94 fESDFileName(esdfilename),
95 fGAlFileName(galfilename),
100 fCheckParticlePID(kFALSE),
101 fReadMostProbableOnly(kFALSE),
105 fITSTrackPoints(kFALSE),
106 fITSTrackPointsType(AliTrackPoints::kITS),
108 fReadCentralBarrel(kTRUE),
113 fTPCChi2PerClustMin(0.0),
114 fTPCChi2PerClustMax(10e5),
139 if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
140 Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
142 /********************************************************************/
144 AliReaderESD::~AliReaderESD()
152 /**********************************************************/
153 Int_t AliReaderESD::ReadNext()
155 //reads next event from fFile
156 //fRunLoader is for reading Kine
158 if (AliVAODParticle::GetDebug())
159 Info("ReadNext","Entered");
161 if (fEventSim == 0x0) fEventSim = new AliAOD();
162 if (fEventRec == 0x0) fEventRec = new AliAOD();
167 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
171 fFile = OpenFile(fCurrentDir);//rl is opened here
174 Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
180 TString esdname = "ESD";
181 esdname+=fCurrentEvent;
182 AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
185 if (AliVAODParticle::GetDebug() > 2 )
187 Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
190 delete fFile;//we have to assume there is no more ESD objects in the fFile
201 return 0;//success -> read one event
202 }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
204 return 1; //no more directories to read
207 /**********************************************************/
208 Int_t AliReaderESD::ReadESD(AliESD* esd)
213 Error("ReadESD","ESD is NULL");
217 // seperate each method
218 if (fReadCentralBarrel) ReadESDCentral(esd);
220 if (fReadMuon) ReadESDMuon(esd);
222 if (fReadPHOS) ReadESDPHOS(esd);
227 /**********************************************************/
228 Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
230 //****** Tentative particle type "concentrations"
231 static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
233 Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
234 Double_t w[kNSpecies];
235 Double_t mom[3];//momentum
236 Double_t pos[3];//position
237 Double_t vertexpos[3];//vertex position
240 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
243 Error("ReadESD","Can not get PDG Database Instance.");
247 Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
249 if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
251 Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
257 Info("ReadESD","Magnetic Field is %f",mf);
258 AliKalmanTrack::SetMagneticField(mf);
261 AliStack* stack = 0x0;
262 if (fReadSim && fRunLoader)
264 fRunLoader->GetEvent(fCurrentEvent);
265 stack = fRunLoader->Stack();
268 const AliESDVertex* vertex = esd->GetVertex();
271 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
278 vertex->GetXYZ(vertexpos);
281 if (AliVAODParticle::GetDebug() > 0)
283 Info("ReadESD","Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]);
286 Info("ReadESD","Reading Event %d",fCurrentEvent);
288 Int_t ntr = esd->GetNumberOfTracks();
289 Info("ReadESD","Found %d tracks.",ntr);
290 for (Int_t i = 0;i<ntr; i++)
292 AliESDtrack *esdtrack = esd->GetTrack(i);
295 Error("Next","Can not get track %d", i);
299 //if (esdtrack->HasVertexParameters() == kFALSE)
300 if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
302 if (AliVAODParticle::GetDebug() > 2)
303 Info("ReadNext","Particle skipped: Data at vertex not available.");
309 if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
311 if (AliVAODParticle::GetDebug() > 2)
312 Info("ReadNext","Particle skipped: Was not reconstructed in TPC.");
317 if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
319 if (AliVAODParticle::GetDebug() > 2)
320 Info("ReadNext","Particle skipped: PID BIT is not set.");
327 esdtrack->GetConstrainedExternalParameters(extx,extp);
330 if (AliVAODParticle::GetDebug() > 2)
331 Info("ReadNext","Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
334 esdtrack->GetESDpid(pidtable);
335 esdtrack->GetConstrainedPxPyPz(mom);
336 esdtrack->GetConstrainedXYZ(pos);
337 pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
338 pos[1] -= vertexpos[1];
339 pos[2] -= vertexpos[2];
341 Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
343 //Particle from kinematics
344 AliAODParticle* particle = 0;
345 Bool_t keeppart = kFALSE;
346 if ( fReadSim && stack )
348 if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
349 TParticle *p = stack->Particle(esdtrack->GetLabel());
352 Error("ReadNext","Can not find track with such label.");
356 if (fCheckParticlePID)
358 if(Rejected(p->GetPdgCode()))
360 if ( AliVAODParticle::GetDebug() > 5 )
361 Info("ReadNext","Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode());
362 continue; //check if we are intersted with particles of this type
365 // if(p->GetPdgCode()<0) charge = -1;
366 particle = new AliAODParticle(*p,i);
370 if(CheckTrack(esdtrack)) continue;
372 //Here we apply Bayes' formula
374 for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
377 if (AliVAODParticle::GetDebug() > 2)
378 Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
382 for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
384 if (AliVAODParticle::GetDebug() > 4)
386 Info("ReadNext","###########################################################################");
387 Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
388 Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
389 TString msg("Pid list got from track:");
390 for (Int_t s = 0;s<kNSpecies;s++)
395 msg+=charge*GetSpeciesPdgCode((ESpecies)s);
402 Info("ReadNext","%s",msg.Data());
403 }//if (AliVAODParticle::GetDebug()>4)
405 AliTrackPoints* tpts = 0x0;
406 if (fNTrackPoints > 0)
408 tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
409 // tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
412 AliTrackPoints* itstpts = 0x0;
415 itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
416 // itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
419 AliClusterMap* cmap = 0x0;
422 cmap = new AliClusterMap(esdtrack);
425 //If this flag fReadMostProbableOnly is false the
426 //loop over species (see "LOOPPIDS") is over all possible PIDs
427 //in other case the most probablle PID is searched
428 //and the loop is limited to that PID only
430 Int_t firstspecie = 0;
431 Int_t lastspecie = kNSpecies;
433 if (fReadMostProbableOnly)
435 //find the most probable PID
437 Float_t maxprob = w[0];
438 for (Int_t s=1; s<AliPID::kSPECIES; s++)
447 lastspecie = spec + 1;
450 for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS
452 Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
456 if ( AliVAODParticle::GetDebug() > 5 )
457 Info("ReadNext","Probability of being PID %d is zero. Continuing.",pdgcode);
461 if(Rejected(pdgcode))
463 if ( AliVAODParticle::GetDebug() > 5 )
464 Info("ReadNext","PID (%d) did not pass the cut.",pdgcode);
465 continue; //check if we are intersted with particles of this type
468 Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
469 Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
471 AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i,
472 mom[0], mom[1], mom[2], tEtot,
473 pos[0], pos[1], pos[2], 0.);
474 //copy probabilitis of other species (if not zero)
475 for (Int_t k = 0; k<kNSpecies; k++)
477 if (k == s) continue;
478 if (w[k] == 0.0) continue;
479 track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
482 if(Rejected(track))//check if meets all criteria of any of our cuts
483 //if it does not delete it and take next good track
485 if ( AliVAODParticle::GetDebug() > 4 )
486 Info("ReadNext","Track did not pass the cut");
491 //Single Particle cuts on cluster map and track points rather do not have sense
494 track->SetTPCTrackPoints(tpts);
499 track->SetITSTrackPoints(itstpts);
504 track->SetClusterMap(cmap);
507 fEventRec->AddParticle(track);
508 if (particle) fEventSim->AddParticle(particle);
511 if (AliVAODParticle::GetDebug() > 4 )
513 Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
515 if (particle) particle->Print();
516 Info("ReadNext","\n----------------------------------------------\n");
518 }//for (Int_t s = 0; s<kNSpecies; s++)
520 if (keeppart == kFALSE)
522 delete particle;//particle was not stored in event
529 if ( fReadSim && stack )
531 if (particle->P() < 0.00001)
533 Info("ReadNext","###################################");
534 Info("ReadNext","###################################");
535 Info("ReadNext","Track Label %d",esdtrack->GetLabel());
536 TParticle *p = stack->Particle(esdtrack->GetLabel());
545 }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
547 Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
548 fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
549 fNEventsRead,fCurrentEvent,fCurrentDir);
550 fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
552 /******************************************************/
553 /****** Setting glevet properties *************/
554 /******************************************************/
556 if (fEventRec->GetNumberOfParticles() > 0)
558 fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
563 /**********************************************************/
564 Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
567 Double_t vertexpos[3];//vertex position, assuming no secondary decay
569 const AliESDVertex* vertex = esd->GetVertex();
572 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
577 vertex->GetXYZ(vertexpos);
580 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
582 if (AliVAODParticle::GetDebug() > 0) {
583 Info("ReadESD","Reading Event %d",fCurrentEvent);
584 Info("ReadESD","Found %d tracks.",nTracks);
587 Float_t Chi2Cut = 100.;
588 Float_t PtCutMin = 1.;
589 Float_t PtCutMax = 10000.;
590 Float_t muonMass = 0.105658389;
592 Double_t thetaX, thetaY, pYZ;
593 Double_t pxRec1, pyRec1, pzRec1, E1;
601 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
603 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
605 thetaX = muonTrack->GetThetaX();
606 thetaY = muonTrack->GetThetaY();
608 pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
609 pzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
610 pxRec1 = pzRec1 * TMath::Tan(thetaX);
611 pyRec1 = pzRec1 * TMath::Tan(thetaY);
612 charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
613 E1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
614 fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, E1);
616 ntrackhits = muonTrack->GetNHit();
617 fitfmin = muonTrack->GetChi2();
619 // transverse momentum
620 Float_t pt1 = fV1.Pt();
623 Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5);
625 if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
626 AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack,
627 pxRec1, pyRec1,pzRec1, E1,
628 vertexpos[0], vertexpos[1], vertexpos[2], 0.);
629 fEventRec->AddParticle(track);
633 fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
637 /**********************************************************/
639 void AliReaderESD::Rewind()
642 // delete fKeyIterator;
645 // fKeyIterator = 0x0;
650 if (fEventList) fEventList->Reset();
651 if (fTrackCounter) fTrackCounter->Reset();
653 /**********************************************************/
655 TFile* AliReaderESD::OpenFile(Int_t n)
657 //opens fFile with tree
666 while (fCurrentDir < n)
674 const TString& dirname = GetDirName(n);
677 Error("OpenFiles","Can not get directory name");
684 filename = fEventList->GetURL(fESDFileName);
688 filename = dirname +"/"+ fESDFileName;
690 Info("OpenFile","%s ==> %s",fESDFileName.Data(),filename.Data());
691 TFile *ret = TFile::Open(filename.Data());
695 Error("OpenFiles","Can't open fFile %s",filename.Data());
700 Error("OpenFiles","Can't open fFile %s",filename.Data());
709 gafilename = fEventList->GetURL(fGAlFileName);
713 gafilename = dirname +"/"+ fGAlFileName;
715 Info("OpenFile","%s ==> %s",fGAlFileName.Data(),gafilename.Data());
717 fRunLoader = AliRunLoader::Open(gafilename);
719 if (fRunLoader == 0x0)
721 Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
726 fRunLoader->LoadHeader();
730 TString kinefilename = fEventList->GetURL("Kinematics.root");
731 fRunLoader->SetKineFileName(kinefilename);
732 Info("OpenFile","%s ==> %s","Kinematics.root",kinefilename.Data());
735 if (fRunLoader->LoadKinematics())
737 Error("Next","Error occured while loading kinematics.");
746 /**********************************************************/
748 Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
750 //returns pdg code from the PID index
751 //ask jura about charge
770 ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
775 /********************************************************************/
776 Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const
778 //Performs check of the track
780 if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE;
782 if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE;
784 if (t->GetTPCclusters(0x0) > 0)
786 Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0));
787 if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE;
791 t->GetConstrainedExternalCovariance(cc);
793 if ( (cc[0] < fC00Min) || (cc[0] > fC00Max) ) return kTRUE;
794 if ( (cc[2] < fC11Min) || (cc[2] > fC11Max) ) return kTRUE;
795 if ( (cc[5] < fC22Min) || (cc[5] > fC22Max) ) return kTRUE;
796 if ( (cc[9] < fC33Min) || (cc[9] > fC33Max) ) return kTRUE;
797 if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
800 t->GetInnerExternalCovariance(cc);
802 if ( (cc[0] < fTPCC00Min) || (cc[0] > fTPCC00Max) ) return kTRUE;
803 if ( (cc[2] < fTPCC11Min) || (cc[2] > fTPCC11Max) ) return kTRUE;
804 if ( (cc[5] < fTPCC22Min) || (cc[5] > fTPCC22Max) ) return kTRUE;
805 if ( (cc[9] < fTPCC33Min) || (cc[9] > fTPCC33Max) ) return kTRUE;
806 if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE;
811 /********************************************************************/
813 void AliReaderESD::SetChi2Range(Float_t min, Float_t max)
815 //sets range of Chi2 per Cluster
819 /********************************************************************/
821 void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
823 //sets range of Number Of Clusters that tracks have to have
827 /********************************************************************/
829 void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
831 //sets range of Chi2 per Cluster
832 fTPCChi2PerClustMin = min;
833 fTPCChi2PerClustMax = max;
835 /********************************************************************/
837 void AliReaderESD::SetC00Range(Float_t min, Float_t max)
839 //Sets range of C00 parameter of covariance matrix of the track
840 //it defines uncertainty of the momentum
844 /********************************************************************/
846 void AliReaderESD::SetC11Range(Float_t min, Float_t max)
848 //Sets range of C11 parameter of covariance matrix of the track
849 //it defines uncertainty of the momentum
853 /********************************************************************/
855 void AliReaderESD::SetC22Range(Float_t min, Float_t max)
857 //Sets range of C22 parameter of covariance matrix of the track
858 //it defines uncertainty of the momentum
862 /********************************************************************/
864 void AliReaderESD::SetC33Range(Float_t min, Float_t max)
866 //Sets range of C33 parameter of covariance matrix of the track
867 //it defines uncertainty of the momentum
871 /********************************************************************/
873 void AliReaderESD::SetC44Range(Float_t min, Float_t max)
875 //Sets range of C44 parameter of covariance matrix of the track
876 //it defines uncertainty of the momentum