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