30394dd38669409c96df4443e73f0614386789f8
[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   //TDatabasePDG* pdgdb = TDatabasePDG::Instance();
87   //if (pdgdb == 0)
88   //{
89   //   Error("ReadESD","Can not get PDG Database Instance.");
90   //   return kFALSE;
91   // }
92    
93   Float_t mf = esd->GetMagneticField(); 
94   if (mf <= 0.0)  
95   {
96      Error("ReadESD","Magnetic Field is 0. Skipping to next event.");
97      return kFALSE;
98   }
99   AliKalmanTrack::SetMagneticField(mf/10.);
100
101   Info("ReadESD","Reading Event %d",fCurrentEvent);
102   if((!fOwner) || (fEventParticles==0)) 
103     fEventParticles = new AliJetEventParticles();
104
105   const Int_t kntr = esd->GetNumberOfTracks();
106   Info("ReadESD","Found %d tracks.",kntr);
107   fEventParticles->Reset(kntr);
108
109   TString headdesc="";
110   headdesc+="Run ";
111   headdesc+=esd->GetRunNumber();
112   headdesc+=": Ev ";
113   headdesc+=esd->GetEventNumber();
114   fEventParticles->SetHeader(headdesc);
115
116   Double_t vertexpos[3];//vertex position
117   const AliESDVertex* kvertex = esd->GetVertex();
118   if (kvertex == 0)
119    {
120      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
121      vertexpos[0] = 0.0;
122      vertexpos[1] = 0.0;
123      vertexpos[2] = 0.0;
124    }
125   else
126    {
127      kvertex->GetXYZ(vertexpos);
128    }
129   fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
130
131   //loop over tracks
132   for (Int_t i = 0;i<kntr; i++)
133    {
134      const AliESDtrack *kesdtrack = esd->GetTrack(i);
135      if (kesdtrack == 0)
136       {
137         Error("ReadESD","Can't get track %d", i);
138         continue;
139       }
140
141      if ((kesdtrack->GetStatus() & fPassFlag)) // != fPassFlag)
142       {
143         Info("ReadNext","Particle skipped: %ud.",kesdtrack->GetStatus());
144         continue;
145       }
146
147      Double_t mom[3];  //momentum
148      Double_t xyz[3];  //position
149      if (fConstrained) {
150        if (kesdtrack->GetConstrainedChi2() > 25) continue;
151        kesdtrack->GetConstrainedPxPyPz(mom);
152        kesdtrack->GetConstrainedXYZ(xyz);
153      } else {
154        kesdtrack->GetPxPyPz(mom);
155        kesdtrack->GetXYZ(xyz);
156      }
157      const Float_t kmass=kesdtrack->GetMass();
158      const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
159      const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
160      const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
161      const Float_t kp=TMath::Sqrt(kp2);
162      const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30)); 
163      const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
164      //Double_t dx = xyz[0]-vertexpos[0];
165      //Double_t dy = xyz[1]-vertexpos[1];
166      //Float_t dca = TMath::Sqrt(dx*dx + dy*dy);
167      //Float_t dz = xyz[2]-vertexpos[2];
168      UInt_t index[6];
169      const Int_t kncl=kesdtrack->GetITSclusters(index)
170                       +kesdtrack->GetTPCclusters(NULL)
171                       +kesdtrack->GetTRDclusters(NULL);
172      if(IsAcceptedParticle(kpt,kphi,keta))
173        fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kncl,kpt,kphi,keta);
174
175    } // loop over tracks
176
177   return kTRUE;
178 }
179
180 /**********************************************************/
181
182 void AliJetParticlesReaderESD::Rewind()
183 {
184   //rewinds reading 
185
186   if(fFile) delete fFile;
187   if(fKeyIterator) delete fKeyIterator;
188   fKeyIterator = 0;
189   fFile = 0;
190   fCurrentDir = 0;
191   fNEventsRead = 0;
192 }
193
194 /**********************************************************/
195
196 Int_t AliJetParticlesReaderESD::ReadNext()
197 {
198   //reads next event from fFile
199
200   do   // is OK even if 0 dirs specified, 
201     {  // in that case we try to read from "./"
202       if (fFile == 0)
203         {
204           fFile = OpenFile(fCurrentDir);
205           if (fFile == 0)
206             {
207               Error("ReadNext","Can't get fFile for dir no. %d",fCurrentDir);
208               fCurrentDir++;
209               continue;
210             }
211      
212           fCurrentEvent = 0;
213           //fFile->GetListOfKeys()->Print();
214           
215           if(fTree) delete fTree;
216           fTree = dynamic_cast<TTree*>(fFile->Get("esdTree"));
217           if(fTree)
218             fTree->SetBranchAddress("ESD",&fESD);
219           else
220             fKeyIterator = new TIter(fFile->GetListOfKeys());  
221         } 
222
223       if(fTree)
224         {
225           if(AliKalmanTrack::GetConvConst()<=0.)
226             AliKalmanTrack::SetMagneticField(0.5);
227           if(fCurrentEvent>=fTree->GetEntries())
228             {
229               fCurrentDir++;
230               delete fTree;
231               fTree = 0;
232               delete fFile;
233               fFile = 0;
234               continue;
235             }
236           fTree->GetEvent(fCurrentEvent);
237         } 
238       else 
239         { // "old" way via ESD objects stored in root file
240           TKey* key = (TKey*)fKeyIterator->Next();
241           if (key == 0)
242             {
243               fCurrentDir++;
244               delete fKeyIterator;
245               fKeyIterator = 0;
246               delete fFile; //we have to assume there are no more ESD objects in the fFile
247               fFile = 0;
248               continue;
249             }
250           TString esdname = "ESD";
251           esdname+=fCurrentEvent;
252           if(fESD) delete fESD;
253           fESD = dynamic_cast<AliESD*>(fFile->Get(esdname));
254           if (fESD == 0)
255             {
256               Info("ReadNext","Can't find AliESD object named %s",esdname.Data());
257               fCurrentDir++;
258               delete fKeyIterator;
259               fKeyIterator = 0;
260               delete fFile;//we have to assume there is no more ESD objects in the fFile
261               fFile = 0;
262               continue;
263             }
264         }
265       ReadESD(fESD);
266       fCurrentEvent++;
267       fNEventsRead++;
268       return kTRUE;//success -> read one event
269     }  while(fCurrentDir < GetNumberOfDirs());
270       //end of loop over directories specified in fDirs Obj Array  
271   return kFALSE; //no more directories to read
272 }
273
274 /**********************************************************/
275
276 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
277 {
278   //opens fFile with kine tree
279
280  const TString& dirname = GetDirName(n);
281  if (dirname == "")
282   {
283    Error("OpenFiles","Can't get directory name");
284    return 0;
285   }
286  TString filename = dirname +"/"+ fESDFileName;
287  TFile *ret = TFile::Open(filename.Data()); 
288  if (ret == 0)
289   {
290     Error("OpenFiles","Can't open fFile %s",filename.Data());
291     return 0;
292   }
293  if (!ret->IsOpen())
294   {
295     Error("OpenFiles","Can't open fFile  %s",filename.Data());
296     return 0;
297   }
298    
299  return ret;
300 }