Bugfixes again for HLT.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderESD.cxx
CommitLineData
d7c6ab14 1// $Id$
2
3//____________________________________________________________________
4//////////////////////////////////////////////////////////////////////
5// //
6// class AliHBTReaderESD //
7// //
8// Reader for ALICE Event Summary Data (ESD). //
9// //
10// Piotr.Skowronski@cern.ch //
11// //
12//////////////////////////////////////////////////////////////////////
13
d9a296b1 14#include <Riostream.h>
d7c6ab14 15#include <TMath.h>
16#include <TPDGCode.h>
17#include <TString.h>
18#include <TObjString.h>
19#include <TTree.h>
20#include <TFile.h>
21#include <TKey.h>
22#include <TH1.h>
23#include <AliESD.h>
24#include <AliESDtrack.h>
04a02430 25#include <AliKalmanTrack.h>
d7c6ab14 26#include <AliJetEventParticles.h>
27#include "AliJetParticlesReaderESD.h"
28
29ClassImp(AliJetParticlesReaderESD)
30
301a24f1 31AliJetParticlesReaderESD::AliJetParticlesReaderESD(Bool_t constrained,
32 const Char_t* esdfilename) :
d7c6ab14 33 AliJetParticlesReader(),
301a24f1 34 fConstrained(constrained),
d7c6ab14 35 fESDFileName(esdfilename),
5d60c8f9 36 fESD(0),
d7c6ab14 37 fFile(0),
5d60c8f9 38 fTree(0),
d7c6ab14 39 fKeyIterator(0),
40 fPassFlag(AliESDtrack::kTPCrefit)
41{
42 //constructor
43}
44
45/********************************************************************/
46
47AliJetParticlesReaderESD::AliJetParticlesReaderESD(
301a24f1 48 Bool_t constrained,
d7c6ab14 49 TObjArray* dirs,
50 const Char_t* esdfilename) :
51 AliJetParticlesReader(dirs),
301a24f1 52 fConstrained(constrained),
d7c6ab14 53 fESDFileName(esdfilename),
5d60c8f9 54 fESD(0),
d7c6ab14 55 fFile(0),
5d60c8f9 56 fTree(0),
d7c6ab14 57 fKeyIterator(0),
58 fPassFlag(AliESDtrack::kTPCrefit)
59{
60 //constructor
61}
62
63/********************************************************************/
64
65AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
66{
67 //desctructor
5d60c8f9 68 if(fESD) delete fESD;
69 if(fTree) delete fTree;
d7c6ab14 70 if(fKeyIterator) delete fKeyIterator;
04a02430 71 if(fFile) delete fFile;
d7c6ab14 72}
73
d7c6ab14 74/**********************************************************/
75
76Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
77{
78 //Reads one ESD
79
80 if (esd == 0)
81 {
82 Error("ReadESD","ESD is NULL");
83 return kFALSE;
84 }
85
53f1300c 86 /*
87 TDatabasePDG* pdgdb = TDatabasePDG::Instance();
88 if (pdgdb == 0)
89 {
90 Error("ReadESD","Can not get PDG Database Instance.");
91 return kFALSE;
92 }
93 */
d7c6ab14 94
04a02430 95 Float_t mf = esd->GetMagneticField();
96 if (mf <= 0.0)
97 {
98 Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
99 return kFALSE;
100 }
101 AliKalmanTrack::SetMagneticField(mf/10.);
d7c6ab14 102
a86edc4e 103 Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
d7c6ab14 104 if((!fOwner) || (fEventParticles==0))
105 fEventParticles = new AliJetEventParticles();
106
04a02430 107 const Int_t kntr = esd->GetNumberOfTracks();
108 Info("ReadESD","Found %d tracks.",kntr);
109 fEventParticles->Reset(kntr);
110
111 TString headdesc="";
112 headdesc+="Run ";
113 headdesc+=esd->GetRunNumber();
114 headdesc+=": Ev ";
115 headdesc+=esd->GetEventNumber();
116 fEventParticles->SetHeader(headdesc);
117
d7c6ab14 118 Double_t vertexpos[3];//vertex position
119 const AliESDVertex* kvertex = esd->GetVertex();
120 if (kvertex == 0)
121 {
122 Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
123 vertexpos[0] = 0.0;
124 vertexpos[1] = 0.0;
125 vertexpos[2] = 0.0;
126 }
127 else
128 {
129 kvertex->GetXYZ(vertexpos);
130 }
131 fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
132
04a02430 133 //loop over tracks
d7c6ab14 134 for (Int_t i = 0;i<kntr; i++)
135 {
bad753b1 136
d7c6ab14 137 const AliESDtrack *kesdtrack = esd->GetTrack(i);
138 if (kesdtrack == 0)
139 {
04a02430 140 Error("ReadESD","Can't get track %d", i);
d7c6ab14 141 continue;
142 }
45cf33d2 143
144 ULong_t cmptest=(kesdtrack->GetStatus() & fPassFlag);
145 if (cmptest!=fPassFlag)
d7c6ab14 146 {
45cf33d2 147 Info("ReadNext","Particle %d skipped: %u.",i,kesdtrack->GetStatus());
9771a1a1 148 cout << i << " "; PrintESDtrack(kesdtrack); cout << endl;
d7c6ab14 149 continue;
150 }
151
152 Double_t mom[3]; //momentum
301a24f1 153 Double_t xyz[3]; //position
154 if (fConstrained) {
155 if (kesdtrack->GetConstrainedChi2() > 25) continue;
156 kesdtrack->GetConstrainedPxPyPz(mom);
157 kesdtrack->GetConstrainedXYZ(xyz);
158 } else {
9771a1a1 159 if(!kesdtrack->GetPxPyPzAt(0,mom)) continue;
53f1300c 160 kesdtrack->GetXYZAt(0, xyz);
301a24f1 161 }
d7c6ab14 162 const Float_t kmass=kesdtrack->GetMass();
163 const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
164 const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
165 const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
166 const Float_t kp=TMath::Sqrt(kp2);
167 const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30));
168 const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
301a24f1 169 //Double_t dx = xyz[0]-vertexpos[0];
170 //Double_t dy = xyz[1]-vertexpos[1];
171 //Float_t dca = TMath::Sqrt(dx*dx + dy*dy);
172 //Float_t dz = xyz[2]-vertexpos[2];
173 UInt_t index[6];
174 const Int_t kncl=kesdtrack->GetITSclusters(index)
175 +kesdtrack->GetTPCclusters(NULL)
176 +kesdtrack->GetTRDclusters(NULL);
d7c6ab14 177 if(IsAcceptedParticle(kpt,kphi,keta))
301a24f1 178 fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
d7c6ab14 179
180 } // loop over tracks
181
182 return kTRUE;
183}
184
185/**********************************************************/
186
187void AliJetParticlesReaderESD::Rewind()
188{
189 //rewinds reading
190
191 if(fFile) delete fFile;
192 if(fKeyIterator) delete fKeyIterator;
193 fKeyIterator = 0;
194 fFile = 0;
195 fCurrentDir = 0;
196 fNEventsRead = 0;
197}
198
199/**********************************************************/
200
d9a296b1 201Int_t AliJetParticlesReaderESD::ReadNext()
202{
203 //reads next event from fFile
204
205 do // is OK even if 0 dirs specified,
206 { // in that case we try to read from "./"
207 if (fFile == 0)
208 {
209 fFile = OpenFile(fCurrentDir);
210 if (fFile == 0)
211 {
212 Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
213 fCurrentDir++;
214 continue;
215 }
216
217 fCurrentEvent = 0;
218 //fFile->GetListOfKeys()->Print();
219
220 if(fTree) delete fTree;
221 fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
222 if(fTree)
223 fTree->SetBranchAddress("ESD",&fESD);
224 else
225 fKeyIterator = new TIter(fFile->GetListOfKeys());
226 }
227
228 if(fTree)
229 {
230 if(AliKalmanTrack::GetConvConst()<=0.)
231 AliKalmanTrack::SetMagneticField(0.5);
232 if(fCurrentEvent>=fTree->GetEntries())
233 {
234 fCurrentDir++;
235 delete fTree;
236 fTree = 0;
237 delete fFile;
238 fFile = 0;
239 continue;
240 }
241 fTree->GetEvent(fCurrentEvent);
242 }
243 else
244 { // "old" way via ESD objects stored in root file
245 TKey* key = (TKey*)fKeyIterator->Next();
246 if (key == 0)
247 {
248 fCurrentDir++;
249 delete fKeyIterator;
250 fKeyIterator = 0;
251 delete fFile; //we have to assume there are no more ESD objects in the fFile
252 fFile = 0;
253 continue;
254 }
255 TString esdname = "ESD";
256 esdname+=fCurrentEvent;
257 if(fESD) delete fESD;
258 fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
259 if (fESD == 0)
260 {
261 Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
262 fCurrentDir++;
263 delete fKeyIterator;
264 fKeyIterator = 0;
265 delete fFile;//we have to assume there is no more ESD objects in the fFile
266 fFile = 0;
267 continue;
268 }
269 }
270 ReadESD(fESD);
271 fCurrentEvent++;
272 fNEventsRead++;
273 return kTRUE;//success -> read one event
274 } while(fCurrentDir < GetNumberOfDirs());
275 //end of loop over directories specified in fDirs Obj Array
276 return kFALSE; //no more directories to read
277}
278
279/**********************************************************/
280
d7c6ab14 281TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
282{
283 //opens fFile with kine tree
284
285 const TString& dirname = GetDirName(n);
286 if (dirname == "")
287 {
04a02430 288 Error("OpenFiles","Can't get directory name");
d7c6ab14 289 return 0;
290 }
291 TString filename = dirname +"/"+ fESDFileName;
292 TFile *ret = TFile::Open(filename.Data());
293 if (ret == 0)
294 {
295 Error("OpenFiles","Can't open fFile %s",filename.Data());
296 return 0;
297 }
298 if (!ret->IsOpen())
299 {
300 Error("OpenFiles","Can't open fFile %s",filename.Data());
301 return 0;
302 }
303
304 return ret;
305}
45cf33d2 306
307/**********************************************************/
308
309void AliJetParticlesReaderESD::PrintESDtrack(const AliESDtrack *kesdtrack) const
310{
311 ULong_t status=kesdtrack->GetStatus();
312 cout << hex << status << dec << ": ";
313 if((status & AliESDtrack::kITSin) == AliESDtrack::kITSin) cout << "ITSin ";
314 if((status & AliESDtrack::kITSout) == AliESDtrack::kITSout) cout << "ITSout ";
315 if((status & AliESDtrack::kITSrefit) == AliESDtrack::kITSrefit) cout << "ITSrefit ";
316 if((status & AliESDtrack::kITSpid) == AliESDtrack::kITSpid) cout << "ITSpid ";
317
318 if((status & AliESDtrack::kTPCin) == AliESDtrack::kTPCin) cout << "TPCin ";
319 if((status & AliESDtrack::kTPCout) == AliESDtrack::kTPCout) cout << "TPCout ";
320 if((status & AliESDtrack::kTPCrefit) == AliESDtrack::kTPCrefit) cout << "TPCrefit ";
321 if((status & AliESDtrack::kTPCpid) == AliESDtrack::kTPCpid) cout << "TPCpid ";
322
323 if((status & AliESDtrack::kTRDin) == AliESDtrack::kTRDin) cout << "TRDin ";
324 if((status & AliESDtrack::kTRDout) == AliESDtrack::kTRDout) cout << "TRDout ";
325 if((status & AliESDtrack::kTRDrefit) == AliESDtrack::kTRDrefit) cout << "TRDrefit ";
326 if((status & AliESDtrack::kTRDpid) == AliESDtrack::kTRDpid) cout << "TRDpid ";
327 if((status & AliESDtrack::kTRDbackup) == AliESDtrack::kTRDbackup) cout << "TRDbackup ";
328 if((status & AliESDtrack::kTRDStop) == AliESDtrack::kTRDStop) cout << "TRDstop ";
329
330 if((status & AliESDtrack::kTOFin) == AliESDtrack::kTOFin) cout << "TOFin ";
331 if((status & AliESDtrack::kTOFout) == AliESDtrack::kTOFout) cout << "TOFout ";
332 if((status & AliESDtrack::kTOFrefit) == AliESDtrack::kTOFrefit) cout << "TOFrefit ";
333 if((status & AliESDtrack::kTOFpid) == AliESDtrack::kTOFpid) cout << "TOFpid ";
334
335 if((status & AliESDtrack::kPHOSpid) == AliESDtrack::kPHOSpid) cout << "PHOSpid ";
336 if((status & AliESDtrack::kRICHpid) == AliESDtrack::kRICHpid) cout << "RICHpid ";
337 if((status & AliESDtrack::kEMCALpid) == AliESDtrack::kEMCALpid) cout << "EMCALpid ";
338 if((status & AliESDtrack::kESDpid) == AliESDtrack::kESDpid) cout << "ESDpid ";
339 if((status & AliESDtrack::kTIME) == AliESDtrack::kTIME) cout << "TIME ";
340 cout << endl;
341 /*
342 cout << kesdtrack->GetConstrainedChi2()
343 << " " << kesdtrack->GetITSchi2()
344 << " " << kesdtrack->GetTPCchi2()
345 << " " << kesdtrack->GetTRDchi2()<< endl;
346 */
347}