Bugfixes again for HLT.
[u/mrichter/AliRoot.git] / JETAN / AliJetParticlesReaderESD.cxx
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
14 #include <Riostream.h>
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>
25 #include <AliKalmanTrack.h>
26 #include <AliJetEventParticles.h>
27 #include "AliJetParticlesReaderESD.h"
28
29 ClassImp(AliJetParticlesReaderESD)
30
31 AliJetParticlesReaderESD::AliJetParticlesReaderESD(Bool_t constrained,
32                                                    const Char_t* esdfilename) :
33   AliJetParticlesReader(),
34   fConstrained(constrained),
35   fESDFileName(esdfilename),
36   fESD(0),
37   fFile(0),
38   fTree(0),
39   fKeyIterator(0),
40   fPassFlag(AliESDtrack::kTPCrefit)
41 {
42   //constructor
43 }
44
45 /********************************************************************/
46   
47 AliJetParticlesReaderESD::AliJetParticlesReaderESD(
48                                       Bool_t constrained,
49                                       TObjArray* dirs,
50                                       const Char_t* esdfilename) :
51   AliJetParticlesReader(dirs),
52   fConstrained(constrained),
53   fESDFileName(esdfilename),
54   fESD(0),
55   fFile(0),
56   fTree(0),
57   fKeyIterator(0),
58   fPassFlag(AliESDtrack::kTPCrefit)
59 {
60   //constructor
61 }
62
63 /********************************************************************/
64
65 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
66 {
67   //desctructor
68   if(fESD) delete fESD;
69   if(fTree) delete fTree;
70   if(fKeyIterator) delete fKeyIterator;
71   if(fFile) delete fFile;
72 }
73
74 /**********************************************************/
75
76 Int_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
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   */
94    
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.);
102
103   Info("ReadESD","Reading Event %d",fCurrentDir*1000+fCurrentEvent);
104   if((!fOwner) || (fEventParticles==0)) 
105     fEventParticles = new AliJetEventParticles();
106
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
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
133   //loop over tracks
134   for (Int_t i = 0;i<kntr; i++)
135    {
136
137      const AliESDtrack *kesdtrack = esd->GetTrack(i);
138      if (kesdtrack == 0)
139       {
140         Error("ReadESD","Can't get track %d", i);
141         continue;
142       }
143      
144      ULong_t cmptest=(kesdtrack->GetStatus() & fPassFlag);
145      if (cmptest!=fPassFlag)
146       {
147         Info("ReadNext","Particle %d skipped: %u.",i,kesdtrack->GetStatus());
148         cout << i << " "; PrintESDtrack(kesdtrack); cout << endl;
149         continue;
150       }
151
152      Double_t mom[3];  //momentum
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 {
159        if(!kesdtrack->GetPxPyPzAt(0,mom)) continue;
160        kesdtrack->GetXYZAt(0, xyz);
161      }
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]);
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);
177      if(IsAcceptedParticle(kpt,kphi,keta))
178        fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
179
180    } // loop over tracks
181
182   return kTRUE;
183 }
184
185 /**********************************************************/
186
187 void 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
201 Int_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
281 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
282 {
283   //opens fFile with kine tree
284
285  const TString& dirname = GetDirName(n);
286  if (dirname == "")
287   {
288    Error("OpenFiles","Can't get directory name");
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 }
306
307 /**********************************************************/
308
309 void 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 }