3 //____________________________________________________________________
4 //////////////////////////////////////////////////////////////////////
6 // class AliHBTReaderESD //
8 // Reader for ALICE Event Summary Data (ESD). //
10 // Piotr.Skowronski@cern.ch //
12 //////////////////////////////////////////////////////////////////////
14 #include <Riostream.h>
18 #include <TObjString.h>
24 #include <AliESDtrack.h>
25 #include <AliKalmanTrack.h>
26 #include <AliJetEventParticles.h>
27 #include "AliJetParticlesReaderESD.h"
29 ClassImp(AliJetParticlesReaderESD)
31 AliJetParticlesReaderESD::AliJetParticlesReaderESD(Bool_t constrained,
32 const Char_t* esdfilename) :
33 AliJetParticlesReader(),
34 fConstrained(constrained),
35 fESDFileName(esdfilename),
40 fPassFlag(AliESDtrack::kTPCrefit)
45 /********************************************************************/
47 AliJetParticlesReaderESD::AliJetParticlesReaderESD(
50 const Char_t* esdfilename) :
51 AliJetParticlesReader(dirs),
52 fConstrained(constrained),
53 fESDFileName(esdfilename),
58 fPassFlag(AliESDtrack::kTPCrefit)
63 /********************************************************************/
65 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
69 if(fTree) delete fTree;
70 if(fKeyIterator) delete fKeyIterator;
71 if(fFile) delete fFile;
74 /**********************************************************/
76 Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
82 Error("ReadESD","ESD is NULL");
86 //TDatabasePDG* pdgdb = TDatabasePDG::Instance();
89 // Error("ReadESD","Can not get PDG Database Instance.");
93 Float_t mf = esd->GetMagneticField();
96 Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
99 AliKalmanTrack::SetMagneticField(mf/10.);
101 Info("ReadESD","Reading Event %d",fCurrentEvent);
102 if((!fOwner) || (fEventParticles==0))
103 fEventParticles = new AliJetEventParticles();
105 const Int_t kntr = esd->GetNumberOfTracks();
106 Info("ReadESD","Found %d tracks.",kntr);
107 fEventParticles->Reset(kntr);
111 headdesc+=esd->GetRunNumber();
113 headdesc+=esd->GetEventNumber();
114 fEventParticles->SetHeader(headdesc);
116 Double_t vertexpos[3];//vertex position
117 const AliESDVertex* kvertex = esd->GetVertex();
120 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
127 kvertex->GetXYZ(vertexpos);
129 fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
132 for (Int_t i = 0;i<kntr; i++)
135 const AliESDtrack *kesdtrack = esd->GetTrack(i);
138 Error("ReadESD","Can't get track %d", i);
142 if ((kesdtrack->GetStatus() & fPassFlag)!=fPassFlag)
144 Info("ReadNext","Particle skipped: %u.",kesdtrack->GetStatus());
148 Double_t mom[3]; //momentum
149 Double_t xyz[3]; //position
151 if (kesdtrack->GetConstrainedChi2() > 25) continue;
152 kesdtrack->GetConstrainedPxPyPz(mom);
153 kesdtrack->GetConstrainedXYZ(xyz);
155 kesdtrack->GetPxPyPz(mom);
156 kesdtrack->GetXYZ(xyz);
158 const Float_t kmass=kesdtrack->GetMass();
159 const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
160 const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
161 const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
162 const Float_t kp=TMath::Sqrt(kp2);
163 const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30));
164 const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
165 //Double_t dx = xyz[0]-vertexpos[0];
166 //Double_t dy = xyz[1]-vertexpos[1];
167 //Float_t dca = TMath::Sqrt(dx*dx + dy*dy);
168 //Float_t dz = xyz[2]-vertexpos[2];
170 const Int_t kncl=kesdtrack->GetITSclusters(index)
171 +kesdtrack->GetTPCclusters(NULL)
172 +kesdtrack->GetTRDclusters(NULL);
173 if(IsAcceptedParticle(kpt,kphi,keta))
174 fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
176 } // loop over tracks
181 /**********************************************************/
183 void AliJetParticlesReaderESD::Rewind()
187 if(fFile) delete fFile;
188 if(fKeyIterator) delete fKeyIterator;
195 /**********************************************************/
197 Int_t AliJetParticlesReaderESD::ReadNext()
199 //reads next event from fFile
201 do // is OK even if 0 dirs specified,
202 { // in that case we try to read from "./"
205 fFile = OpenFile(fCurrentDir);
208 Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
214 //fFile->GetListOfKeys()->Print();
216 if(fTree) delete fTree;
217 fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
219 fTree->SetBranchAddress("ESD",&fESD);
221 fKeyIterator = new TIter(fFile->GetListOfKeys());
226 if(AliKalmanTrack::GetConvConst()<=0.)
227 AliKalmanTrack::SetMagneticField(0.5);
228 if(fCurrentEvent>=fTree->GetEntries())
237 fTree->GetEvent(fCurrentEvent);
240 { // "old" way via ESD objects stored in root file
241 TKey* key = (TKey*)fKeyIterator->Next();
247 delete fFile; //we have to assume there are no more ESD objects in the fFile
251 TString esdname = "ESD";
252 esdname+=fCurrentEvent;
253 if(fESD) delete fESD;
254 fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
257 Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
261 delete fFile;//we have to assume there is no more ESD objects in the fFile
269 return kTRUE;//success -> read one event
270 } while(fCurrentDir < GetNumberOfDirs());
271 //end of loop over directories specified in fDirs Obj Array
272 return kFALSE; //no more directories to read
275 /**********************************************************/
277 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
279 //opens fFile with kine tree
281 const TString& dirname = GetDirName(n);
284 Error("OpenFiles","Can't get directory name");
287 TString filename = dirname +"/"+ fESDFileName;
288 TFile *ret = TFile::Open(filename.Data());
291 Error("OpenFiles","Can't open fFile %s",filename.Data());
296 Error("OpenFiles","Can't open fFile %s",filename.Data());