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