1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 //____________________________________________________________________
19 //////////////////////////////////////////////////////////////////////
21 // class AliReaderESD //
23 // reader for ALICE Event Summary Data (ESD). //
25 // Piotr.Skowronski@cern.ch //
27 //////////////////////////////////////////////////////////////////////
30 #include <TGliteXmlEventlist.h>
33 #include <TObjString.h>
35 #include <TParticle.h>
39 #include "AliAnalysis.h"
40 #include "AliAODRun.h"
43 #include "AliAODParticle.h"
44 #include "AliAODParticleCut.h"
45 #include "AliAODRun.h"
46 #include "AliAnalysis.h"
47 #include "AliClusterMap.h"
49 #include "AliESDtrack.h"
50 #include "AliKalmanTrack.h"
52 #include "AliReaderESD.h"
54 #include "AliRunLoader.h"
57 ClassImp(AliReaderESD)
59 AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
60 fESDFileName(esdfilename),
61 fGAlFileName(galfilename),
66 fCheckParticlePID(kFALSE),
67 fReadMostProbableOnly(kFALSE),
71 fITSTrackPoints(kFALSE),
72 fITSTrackPointsType(AliTrackPoints::kITS),
74 fReadCentralBarrel(kTRUE),
79 fTPCChi2PerClustMin(0.0),
80 fTPCChi2PerClustMax(10e5),
106 if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
107 Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
109 /********************************************************************/
111 AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
113 fESDFileName(esdfilename),
114 fGAlFileName(galfilename),
119 fCheckParticlePID(kFALSE),
120 fReadMostProbableOnly(kFALSE),
124 fITSTrackPoints(kFALSE),
125 fITSTrackPointsType(AliTrackPoints::kITS),
127 fReadCentralBarrel(kTRUE),
132 fTPCChi2PerClustMin(0.0),
133 fTPCChi2PerClustMax(10e5),
158 if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
159 Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
161 /********************************************************************/
163 AliReaderESD::~AliReaderESD()
171 /**********************************************************/
172 Int_t AliReaderESD::ReadNext()
174 //reads next event from fFile
175 //fRunLoader is for reading Kine
177 AliDebug(1,"Entered");
179 if (fEventSim == 0x0) fEventSim = new AliAOD();
180 if (fEventRec == 0x0) fEventRec = new AliAOD();
185 do //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
189 fFile = OpenFile(fCurrentDir);//rl is opened here
192 Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
198 TString esdname = "ESD";
199 esdname+=fCurrentEvent;
200 AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
203 AliDebug(3,Form("Can not find AliESD object named %s",esdname.Data()));
205 delete fFile;//we have to assume there is no more ESD objects in the fFile
216 return 0;//success -> read one event
217 }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
219 return 1; //no more directories to read
222 /**********************************************************/
223 Int_t AliReaderESD::ReadESD(AliESD* esd)
228 Error("ReadESD","ESD is NULL");
232 // seperate each method
233 if (fReadCentralBarrel) ReadESDCentral(esd);
235 if (fReadMuon) ReadESDMuon(esd);
237 if (fReadPHOS) ReadESDPHOS(esd);
242 /**********************************************************/
243 Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
245 //****** Tentative particle type "concentrations"
246 static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
248 Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
249 Double_t w[kNSpecies];
250 Double_t mom[3];//momentum
251 Double_t pos[3];//position
252 Double_t vertexpos[3];//vertex position
255 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
258 Error("ReadESD","Can not get PDG Database Instance.");
262 Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
264 if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
266 Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
272 Info("ReadESD","Magnetic Field is %f",mf);
273 //AliKalmanTrack::SetMagneticField(mf);
276 AliStack* stack = 0x0;
277 if (fReadSim && fRunLoader)
279 fRunLoader->GetEvent(fCurrentEvent);
280 stack = fRunLoader->Stack();
283 const AliESDVertex* vertex = esd->GetVertex();
286 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
293 vertex->GetXYZ(vertexpos);
296 AliDebug(1,Form("Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]));
298 Info("ReadESD","Reading Event %d",fCurrentEvent);
300 Int_t ntr = esd->GetNumberOfTracks();
301 Info("ReadESD","Found %d tracks.",ntr);
302 for (Int_t i = 0;i<ntr; i++)
304 AliESDtrack *esdtrack = esd->GetTrack(i);
307 Error("Next","Can not get track %d", i);
311 //if (esdtrack->HasVertexParameters() == kFALSE)
312 if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
314 AliDebug(3,Form("Particle skipped: Data at vertex not available."));
320 if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
322 AliDebug(3,"Particle skipped: Was not reconstructed in TPC.");
327 if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE)
329 AliDebug(3,"Particle skipped: PID BIT is not set.");
336 esdtrack->GetConstrainedExternalParameters(extx,extp);
339 AliDebug(3,"Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
342 esdtrack->GetESDpid(pidtable);
343 esdtrack->GetConstrainedPxPyPz(mom);
344 esdtrack->GetConstrainedXYZ(pos);
345 pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
346 pos[1] -= vertexpos[1];
347 pos[2] -= vertexpos[2];
349 Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
351 //Particle from kinematics
352 AliAODParticle* particle = 0;
353 Bool_t keeppart = kFALSE;
354 if ( fReadSim && stack )
356 if (esdtrack->GetLabel() < 0) continue;//this is fake - we are not able to match any track
357 TParticle *p = stack->Particle(esdtrack->GetLabel());
360 Error("ReadNext","Can not find track with such label.");
364 if (fCheckParticlePID)
366 if(Rejected(p->GetPdgCode()))
368 AliDebug(6,Form("Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()));
369 continue; //check if we are intersted with particles of this type
372 // if(p->GetPdgCode()<0) charge = -1;
373 particle = new AliAODParticle(*p,i);
377 if(CheckTrack(esdtrack)) continue;
379 //Here we apply Bayes' formula
381 for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
384 AliDebug(3,"Particle rejected since total bayessian PID probab. is zero.");
388 for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
390 if (AliDebugLevel() > 4)
392 Info("ReadNext","###########################################################################");
393 Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
394 Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
395 TString msg("Pid list got from track:");
396 for (Int_t s = 0;s<kNSpecies;s++)
401 msg+=charge*GetSpeciesPdgCode((ESpecies)s);
408 Info("ReadNext","%s",msg.Data());
409 }//if (AliDebugLevel()>4)
411 AliTrackPoints* tpts = 0x0;
412 if (fNTrackPoints > 0)
414 tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
415 // tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
418 AliTrackPoints* itstpts = 0x0;
421 itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
422 // itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
425 AliClusterMap* cmap = 0x0;
428 cmap = new AliClusterMap(esdtrack);
431 //If this flag fReadMostProbableOnly is false the
432 //loop over species (see "LOOPPIDS") is over all possible PIDs
433 //in other case the most probablle PID is searched
434 //and the loop is limited to that PID only
436 Int_t firstspecie = 0;
437 Int_t lastspecie = kNSpecies;
439 if (fReadMostProbableOnly)
441 //find the most probable PID
443 Float_t maxprob = w[0];
444 for (Int_t s=1; s<AliPID::kSPECIES; s++)
453 lastspecie = spec + 1;
456 for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS
458 Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
462 AliDebug(6,Form("Probability of being PID %d is zero. Continuing.",pdgcode));
466 if(Rejected(pdgcode))
468 AliDebug(6,Form("PID (%d) did not pass the cut.",pdgcode));
469 continue; //check if we are intersted with particles of this type
472 Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
473 Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
475 AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i,
476 mom[0], mom[1], mom[2], tEtot,
477 pos[0], pos[1], pos[2], 0.);
478 //copy probabilitis of other species (if not zero)
479 for (Int_t k = 0; k<kNSpecies; k++)
481 if (k == s) continue;
482 if (w[k] == 0.0) continue;
483 track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
486 if(Rejected(track))//check if meets all criteria of any of our cuts
487 //if it does not delete it and take next good track
489 AliDebug(5,"Track did not pass the cut");
494 //Single Particle cuts on cluster map and track points rather do not have sense
497 track->SetTPCTrackPoints(tpts);
502 track->SetITSTrackPoints(itstpts);
507 track->SetClusterMap(cmap);
510 fEventRec->AddParticle(track);
511 if (particle) fEventSim->AddParticle(particle);
514 if (AliDebugLevel() > 4 )
516 Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
518 if (particle) particle->Print();
519 Info("ReadNext","\n----------------------------------------------\n");
521 }//for (Int_t s = 0; s<kNSpecies; s++)
523 if (keeppart == kFALSE)
525 delete particle;//particle was not stored in event
532 if ( fReadSim && stack )
534 if (particle->P() < 0.00001)
536 Info("ReadNext","###################################");
537 Info("ReadNext","###################################");
538 Info("ReadNext","Track Label %d",esdtrack->GetLabel());
539 TParticle *p = stack->Particle(esdtrack->GetLabel());
548 }//for (Int_t i = 0;i<ntr; i++) -- loop over tracks
550 Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
551 fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
552 fNEventsRead,fCurrentEvent,fCurrentDir);
553 fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
555 /******************************************************/
556 /****** Setting glevet properties *************/
557 /******************************************************/
559 if (fEventRec->GetNumberOfParticles() > 0)
561 fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
566 /**********************************************************/
567 Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
570 Double_t vertexpos[3];//vertex position, assuming no secondary decay
572 const AliESDVertex* vertex = esd->GetVertex();
575 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
580 vertex->GetXYZ(vertexpos);
583 Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
585 AliDebug(1,Form("Reading Event %d \nFound %d tracks.",fCurrentEvent,nTracks));
588 Float_t Chi2Cut = 100.;
589 Float_t PtCutMin = 1.;
590 Float_t PtCutMax = 10000.;
591 Float_t muonMass = 0.105658389;
593 Double_t thetaX, thetaY, pYZ;
594 Double_t pxRec1, pyRec1, pzRec1, E1;
602 for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
604 AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
606 thetaX = muonTrack->GetThetaX();
607 thetaY = muonTrack->GetThetaY();
609 pYZ = 1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
610 pzRec1 = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
611 pxRec1 = pzRec1 * TMath::Tan(thetaX);
612 pyRec1 = pzRec1 * TMath::Tan(thetaY);
613 charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
614 E1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
615 fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, E1);
617 ntrackhits = muonTrack->GetNHit();
618 fitfmin = muonTrack->GetChi2();
620 // transverse momentum
621 Float_t pt1 = fV1.Pt();
624 Float_t ch1 = fitfmin / (2.0 * ntrackhits - 5);
626 if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
627 AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack,
628 pxRec1, pyRec1,pzRec1, E1,
629 vertexpos[0], vertexpos[1], vertexpos[2], 0.);
630 fEventRec->AddParticle(track);
634 fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
638 /**********************************************************/
640 void AliReaderESD::Rewind()
643 // delete fKeyIterator;
646 // fKeyIterator = 0x0;
651 if (fEventList) fEventList->Reset();
652 if (fTrackCounter) fTrackCounter->Reset();
654 /**********************************************************/
656 TFile* AliReaderESD::OpenFile(Int_t n)
658 //opens fFile with tree
667 while (fCurrentDir < n)
675 const TString& dirname = GetDirName(n);
678 Error("OpenFiles","Can not get directory name");
685 filename = fEventList->GetURL(fESDFileName);
689 filename = dirname +"/"+ fESDFileName;
691 Info("OpenFile","%s ==> %s",fESDFileName.Data(),filename.Data());
692 TFile *ret = TFile::Open(filename.Data());
696 Error("OpenFiles","Can't open fFile %s",filename.Data());
701 Error("OpenFiles","Can't open fFile %s",filename.Data());
710 gafilename = fEventList->GetURL(fGAlFileName);
714 gafilename = dirname +"/"+ fGAlFileName;
716 Info("OpenFile","%s ==> %s",fGAlFileName.Data(),gafilename.Data());
718 fRunLoader = AliRunLoader::Open(gafilename);
720 if (fRunLoader == 0x0)
722 Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
727 fRunLoader->LoadHeader();
731 TString kinefilename = fEventList->GetURL("Kinematics.root");
732 fRunLoader->SetKineFileName(kinefilename);
733 Info("OpenFile","%s ==> %s","Kinematics.root",kinefilename.Data());
736 if (fRunLoader->LoadKinematics())
738 Error("Next","Error occured while loading kinematics.");
747 /**********************************************************/
749 Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
751 //returns pdg code from the PID index
752 //ask jura about charge
771 ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
776 /********************************************************************/
777 Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const
779 //Performs check of the track
781 if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE;
783 if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE;
785 if (t->GetTPCclusters(0x0) > 0)
787 Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0));
788 if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE;
792 t->GetConstrainedExternalCovariance(cc);
794 if ( (cc[0] < fC00Min) || (cc[0] > fC00Max) ) return kTRUE;
795 if ( (cc[2] < fC11Min) || (cc[2] > fC11Max) ) return kTRUE;
796 if ( (cc[5] < fC22Min) || (cc[5] > fC22Max) ) return kTRUE;
797 if ( (cc[9] < fC33Min) || (cc[9] > fC33Max) ) return kTRUE;
798 if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
801 t->GetInnerExternalCovariance(cc);
803 if ( (cc[0] < fTPCC00Min) || (cc[0] > fTPCC00Max) ) return kTRUE;
804 if ( (cc[2] < fTPCC11Min) || (cc[2] > fTPCC11Max) ) return kTRUE;
805 if ( (cc[5] < fTPCC22Min) || (cc[5] > fTPCC22Max) ) return kTRUE;
806 if ( (cc[9] < fTPCC33Min) || (cc[9] > fTPCC33Max) ) return kTRUE;
807 if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE;
812 /********************************************************************/
814 void AliReaderESD::SetChi2Range(Float_t min, Float_t max)
816 //sets range of Chi2 per Cluster
820 /********************************************************************/
822 void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
824 //sets range of Number Of Clusters that tracks have to have
828 /********************************************************************/
830 void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
832 //sets range of Chi2 per Cluster
833 fTPCChi2PerClustMin = min;
834 fTPCChi2PerClustMax = max;
836 /********************************************************************/
838 void AliReaderESD::SetC00Range(Float_t min, Float_t max)
840 //Sets range of C00 parameter of covariance matrix of the track
841 //it defines uncertainty of the momentum
845 /********************************************************************/
847 void AliReaderESD::SetC11Range(Float_t min, Float_t max)
849 //Sets range of C11 parameter of covariance matrix of the track
850 //it defines uncertainty of the momentum
854 /********************************************************************/
856 void AliReaderESD::SetC22Range(Float_t min, Float_t max)
858 //Sets range of C22 parameter of covariance matrix of the track
859 //it defines uncertainty of the momentum
863 /********************************************************************/
865 void AliReaderESD::SetC33Range(Float_t min, Float_t max)
867 //Sets range of C33 parameter of covariance matrix of the track
868 //it defines uncertainty of the momentum
872 /********************************************************************/
874 void AliReaderESD::SetC44Range(Float_t min, Float_t max)
876 //Sets range of C44 parameter of covariance matrix of the track
877 //it defines uncertainty of the momentum