Added ESD reader for HLT tracks.
[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 <TMath.h>
15 #include <TPDGCode.h>
16 #include <TString.h>
17 #include <TObjString.h>
18 #include <TTree.h>
19 #include <TFile.h>
20 #include <TKey.h>
21 #include <TH1.h>
22 #include <AliESD.h>
23 #include <AliESDtrack.h>
24 #include <AliJetEventParticles.h>
25 #include "AliJetParticlesReaderESD.h"
26
27 ClassImp(AliJetParticlesReaderESD)
28
29 AliJetParticlesReaderESD::AliJetParticlesReaderESD(const Char_t* esdfilename) :
30   AliJetParticlesReader(),
31   fESDFileName(esdfilename),
32   fESD(0),
33   fFile(0),
34   fTree(0),
35   fKeyIterator(0),
36   fPassFlag(AliESDtrack::kTPCrefit)
37 {
38   //constructor
39 }
40
41 /********************************************************************/
42   
43 AliJetParticlesReaderESD::AliJetParticlesReaderESD(
44                                       TObjArray* dirs,
45                                       const Char_t* esdfilename) :
46   AliJetParticlesReader(dirs),
47   fESDFileName(esdfilename),
48   fESD(0),
49   fFile(0),
50   fTree(0),
51   fKeyIterator(0),
52   fPassFlag(AliESDtrack::kTPCrefit)
53 {
54   //constructor
55 }
56
57 /********************************************************************/
58
59 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
60 {
61   //desctructor
62   if(fFile) delete fFile;
63   if(fESD) delete fESD;
64   if(fTree) delete fTree;
65   if(fKeyIterator) delete fKeyIterator;
66 }
67
68
69 Int_t AliJetParticlesReaderESD::ReadNext()
70 {
71   //reads next event from fFile
72
73   do   // is OK even if 0 dirs specified, 
74     {  // in that case we try to read from "./"
75       if (fFile == 0)
76         {
77           fFile = OpenFile(fCurrentDir);
78           if (fFile == 0)
79             {
80               Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
81               fCurrentDir++;
82               continue;
83             }
84      
85           fCurrentEvent = 0;
86           //fFile->GetListOfKeys()->Print();
87           
88           if(fTree) delete fTree;
89           fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
90           if(fTree)
91             fTree->SetBranchAddress("ESD",&fESD);
92           else
93             fKeyIterator = new TIter(fFile->GetListOfKeys());  
94         } 
95
96       if(fTree)
97         {
98           if(fCurrentEvent>=fTree->GetEntries())
99             {
100               fCurrentDir++;
101               delete fTree;
102               fTree = 0;
103               delete fFile;
104               fFile = 0;
105               continue;
106             }
107           fTree->GetEvent(fCurrentEvent);
108         } 
109       else 
110         { // "old" way via ESD objects stored in root file
111           TKey* key = (TKey*)fKeyIterator->Next();
112           if (key == 0)
113             {
114               fCurrentDir++;
115               delete fKeyIterator;
116               fKeyIterator = 0;
117               delete fFile; //we have to assume there are no more ESD objects in the fFile
118               fFile = 0;
119               continue;
120             }
121           TString esdname = "ESD";
122           esdname+=fCurrentEvent;
123           if(fESD) delete fESD;
124           fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
125           if (fESD == 0)
126             {
127               Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
128               fCurrentDir++;
129               delete fKeyIterator;
130               fKeyIterator = 0;
131               delete fFile;//we have to assume there is no more ESD objects in the fFile
132               fFile = 0;
133               continue;
134             }
135         }
136       ReadESD(fESD);
137       fCurrentEvent++;
138       fNEventsRead++;
139       return kTRUE;//success -> read one event
140     }  while(fCurrentDir < GetNumberOfDirs());
141       //end of loop over directories specified in fDirs Obj Array  
142   return kFALSE; //no more directories to read
143 }
144
145 /**********************************************************/
146
147 Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
148 {
149   //Reads one ESD
150
151   if (esd == 0)
152    {
153      Error("ReadESD","ESD is NULL");
154      return kFALSE;
155    }
156
157   //TDatabasePDG* pdgdb = TDatabasePDG::Instance();
158   //if (pdgdb == 0)
159   //{
160   //   Error("ReadESD","Can not get PDG Database Instance.");
161   //   return kFALSE;
162   // }
163    
164   //Float_t mf = esd->GetMagneticField(); 
165   //if ( (mf == 0.0) && (fNTrackPoints > 0) )
166   //{
167   //   Error("ReadESD","Magnetic Field is 0 and Track Points demanded. Skipping to next event.");
168   //   return kFALSE;
169   //}
170
171   Info("ReadESD","Reading Event %d",fCurrentEvent);
172   if((!fOwner) || (fEventParticles==0)) 
173     fEventParticles = new AliJetEventParticles();
174
175   Double_t vertexpos[3];//vertex position
176   const AliESDVertex* kvertex = esd->GetVertex();
177   if (kvertex == 0)
178    {
179      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
180      vertexpos[0] = 0.0;
181      vertexpos[1] = 0.0;
182      vertexpos[2] = 0.0;
183    }
184   else
185    {
186      kvertex->GetXYZ(vertexpos);
187    }
188   fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
189
190   const Int_t kntr = esd->GetNumberOfTracks();
191   Info("ReadESD","Found %d tracks.",kntr);
192   for (Int_t i = 0;i<kntr; i++)
193    {
194      const AliESDtrack *kesdtrack = esd->GetTrack(i);
195      if (kesdtrack == 0)
196       {
197         Error("ReadESD","Can not get track %d", i);
198         continue;
199       }
200
201      if ((kesdtrack->GetStatus() & fPassFlag) == kFALSE)
202       {
203         Info("ReadNext","Particle skipped: %ud.",kesdtrack->GetStatus());
204         continue;
205       }
206
207      Double_t mom[3];  //momentum
208      kesdtrack->GetPxPyPz(mom);
209      //kesdtrack->GetConstrainedPxPyPz(mom);
210      //Double_t pos[3];//position
211      //kesdtrack->GetXYZ(pos);
212      //kesdtrack->GetConstrainedXYZ(pos);
213      const Float_t kmass=kesdtrack->GetMass();
214      const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
215      const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
216      const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
217      const Float_t kp=TMath::Sqrt(kp2);
218      const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30)); 
219      const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
220      if(IsAcceptedParticle(kpt,kphi,keta))
221        fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kpt,kphi,keta);
222
223    } // loop over tracks
224
225   return kTRUE;
226 }
227
228 /**********************************************************/
229
230 void AliJetParticlesReaderESD::Rewind()
231 {
232   //rewinds reading 
233
234   if(fFile) delete fFile;
235   if(fKeyIterator) delete fKeyIterator;
236   fKeyIterator = 0;
237   fFile = 0;
238   fCurrentDir = 0;
239   fNEventsRead = 0;
240 }
241
242 /**********************************************************/
243
244 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
245 {
246   //opens fFile with kine tree
247
248  const TString& dirname = GetDirName(n);
249  if (dirname == "")
250   {
251    Error("OpenFiles","Can not get directory name");
252    return 0;
253   }
254  TString filename = dirname +"/"+ fESDFileName;
255  TFile *ret = TFile::Open(filename.Data()); 
256  if (ret == 0)
257   {
258     Error("OpenFiles","Can't open fFile %s",filename.Data());
259     return 0;
260   }
261  if (!ret->IsOpen())
262   {
263     Error("OpenFiles","Can't open fFile  %s",filename.Data());
264     return 0;
265   }
266    
267  return ret;
268 }