Removing some temporary solutions for the track parameters at the primary vertex...
[u/mrichter/AliRoot.git] / HBTAN / AliHBTReaderESD.cxx
1 #include "AliHBTReaderESD.h"
2 //____________________________________________________________________
3 //////////////////////////////////////////////////////////////////////
4 //                                                                  //
5 // class AliHBTReaderESD                                            //
6 //                                                                  //
7 // reader for ALICE Event Summary Data (ESD).                       //
8 //                                                                  //
9 // Piotr.Skowronski@cern.ch                                         //
10 //                                                                  //
11 //////////////////////////////////////////////////////////////////////
12
13 #include <TPDGCode.h>
14 #include <TString.h>
15 #include <TObjString.h>
16 #include <TTree.h>
17 #include <TFile.h>
18 #include <TParticle.h>
19
20 #include <AliRunLoader.h>
21 #include <AliStack.h>
22 #include <AliESDtrack.h>
23 #include <AliESD.h>
24
25 #include "AliHBTRun.h"
26 #include "AliHBTEvent.h"
27 #include "AliHBTParticle.h"
28 #include "AliHBTParticleCut.h"
29
30 ClassImp(AliHBTReaderESD)
31
32 AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
33  fESDFileName(esdfilename),
34  fGAlFileName(galfilename),
35  fFile(0x0),
36  fRunLoader(0x0),
37  fReadParticles(kFALSE)
38 {
39   //cosntructor
40   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
41     Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
42 }
43 /********************************************************************/
44   
45 AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
46  AliHBTReader(dirs), 
47  fESDFileName(esdfilename),
48  fGAlFileName(galfilename),
49  fFile(0x0),
50  fRunLoader(0x0),
51  fReadParticles(kFALSE)
52 {
53   //cosntructor
54   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
55     Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
56 }
57 /********************************************************************/
58
59 AliHBTReaderESD::~AliHBTReaderESD()
60 {
61  //desctructor
62   delete fRunLoader;
63   delete fFile;
64 }
65 /**********************************************************/
66 Int_t AliHBTReaderESD::ReadNext()
67 {
68 //reads next event from fFile
69   
70   //****** Tentative particle type "concentrations"
71   static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
72   
73   Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
74   Double_t w[kNSpecies];
75   Double_t mom[3];//momentum
76   Double_t pos[3];//position
77   
78   if (AliHBTParticle::fgDebug)
79     Info("ReadNext","Entered");
80     
81   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
82   if (pdgdb == 0x0)
83    {
84      Error("Next","Can not get PDG Database Instance.");
85      return 1;
86    }
87   
88   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
89   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
90   
91   fParticlesEvent->Reset();
92   fTracksEvent->Reset();
93         
94   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
95    {
96      if (fFile == 0x0)
97       {
98        fFile = OpenFile(fCurrentDir);//rl is opened here
99        if (fFile == 0x0)
100          {
101            Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
102            fCurrentDir++;
103            continue;
104          }
105        fCurrentEvent = 0;
106       } 
107       
108      //try to read
109      TString esdname;
110      esdname+=fCurrentEvent;
111      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
112      if (esd == 0x0)
113       {
114         if (AliHBTParticle::fgDebug > 2 )
115           {
116             Info("ReadNext","Can not find ESD object named %s.",esdname.Data());
117           }
118         fCurrentDir++;
119         delete fFile;//we have to assume there is no more ESD objects in the fFile
120         delete fRunLoader;
121         fFile = 0x0;
122         fRunLoader = 0x0;
123         continue;
124       }
125      AliStack* stack = 0x0;
126      if (fReadParticles && fRunLoader)
127       {
128         fRunLoader->GetEvent(fCurrentEvent);
129         stack = fRunLoader->Stack();
130       }
131      Info("ReadNext","Reading Event %d",fCurrentEvent);
132   
133      Int_t ntr = esd->GetNumberOfTracks();
134      Info("ReadNext","Found %d tracks.",ntr);
135      for (Int_t i = 0;i<ntr; i++)
136       {
137         AliESDtrack *esdtrack = esd->GetTrack(i);
138         if (esdtrack == 0x0)
139          {
140            Error("Next","Can not get track %d", i);
141            continue;
142          }
143         
144         //if (esdtrack->HasVertexParameters() == kFALSE) 
145         if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
146          {
147            if (AliHBTParticle::fgDebug > 2) 
148              Info("ReadNext","Particle skipped: Data at vertex not available.");
149            continue;
150          }
151          
152         if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
153          {
154            if (AliHBTParticle::fgDebug > 2) 
155              Info("ReadNext","Particle skipped: PID BIT is not set.");
156            continue;
157          }
158        
159         esdtrack->GetESDpid(pidtable);
160         //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]); 
161         esdtrack->GetPxPyPz(mom);
162         //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
163         esdtrack->GetXYZ(pos);
164         
165         Double_t extx;
166         Double_t extp[5];
167         esdtrack->GetExternalParameters(extx,extp);
168         
169         Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
170         
171         //Particle from kinematics
172         AliHBTParticle* particle = 0;
173         Bool_t keeppart = kFALSE;
174         if ( fReadParticles && stack  )
175          {
176            if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
177            TParticle *p = stack->Particle(esdtrack->GetLabel());
178            if (p==0x0) 
179             {
180               Error("ReadNext","Can not find track with such label.");
181               continue;
182             }
183 //           if(p->GetPdgCode()<0) charge = -1;
184            particle = new AliHBTParticle(*p,i);
185            
186          }
187          
188         //Here we apply Bayes' formula
189         Double_t rc=0.;
190         for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
191         if (rc==0.0) 
192          {
193            if (AliHBTParticle::fgDebug > 2) 
194              Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
195            continue;
196          }
197
198         for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
199
200         if (AliHBTParticle::fgDebug > 4)
201          { 
202            Info("ReadNext","###########################################################################");
203            Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
204            Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
205            TString msg("Pid list got from track:");
206            for (Int_t s = 0;s<kNSpecies;s++)
207             {
208               msg+="\n    ";
209               msg+=s;
210               msg+="(";
211               msg+=charge*GetSpeciesPdgCode((ESpecies)s);
212               msg+="): ";
213               msg+=w[s];
214               msg+=" (";
215               msg+=pidtable[s];
216               msg+=")";
217             }
218            Info("ReadNext","%s",msg.Data());
219          }//if (AliHBTParticle::fgDebug>4)
220
221         for (Int_t s = 0; s<kNSpecies; s++)
222          {
223            Float_t pp = w[s];
224            if (pp == 0.0) continue;
225
226            Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
227            if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 
228
229            Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
230            Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
231
232            AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
233                                                       mom[0], mom[1], mom[2], tEtot,
234                                                       pos[0], pos[1], pos[2], 0.);
235            //copy probabilitis of other species (if not zero)
236            for (Int_t k = 0; k<kNSpecies; k++)
237             {
238               if (k == s) continue;
239               if (w[k] == 0.0) continue;
240               track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
241             }
242            if(Pass(track))//check if meets all criteria of any of our cuts
243                           //if it does not delete it and take next good track
244             { 
245               delete track;
246               continue;
247             }
248            
249            fTracksEvent->AddParticle(track);
250            if (particle) fParticlesEvent->AddParticle(particle);
251            keeppart = kTRUE;
252                
253            if (AliHBTParticle::fgDebug > 4 )
254             {
255               Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
256               track->Print();
257             }
258          }//for (Int_t s = 0; s<kNSpecies; s++)
259          
260         if (keeppart == kFALSE) delete particle;//particle was not stored in event
261         
262       }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
263      
264      Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
265             fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
266             fNEventsRead,fCurrentEvent,fCurrentDir);
267       
268      fCurrentEvent++;
269      fNEventsRead++;
270      delete esd;
271      return 0;//success -> read one event
272    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
273    
274   return 1; //no more directories to read
275 }
276 /**********************************************************/
277 void AliHBTReaderESD::Rewind()
278 {
279   //rewinds reading 
280   delete fFile;
281   fFile = 0;
282   delete fRunLoader;
283   fCurrentDir = 0;
284   fNEventsRead = 0;
285   fCurrentEvent++;
286 }
287 /**********************************************************/
288
289 TFile* AliHBTReaderESD::OpenFile(Int_t n)
290 {
291 //opens fFile with kine tree
292
293  const TString& dirname = GetDirName(n);
294  if (dirname == "")
295   {
296    Error("OpenFiles","Can not get directory name");
297    return 0x0;
298   }
299  TString filename = dirname +"/"+ fESDFileName;
300  TFile *ret = TFile::Open(filename.Data()); 
301
302  if ( ret == 0x0)
303   {
304     Error("OpenFiles","Can't open fFile %s",filename.Data());
305     return 0x0;
306   }
307  if (!ret->IsOpen())
308   {
309     Error("OpenFiles","Can't open fFile  %s",filename.Data());
310     return 0x0;
311   }
312  
313  if (fReadParticles)
314   {
315    fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
316    if (fRunLoader == 0x0)
317     {
318       Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
319       delete ret;
320       return 0x0;
321     }
322    fRunLoader->LoadHeader();
323    if (fRunLoader->LoadKinematics())
324     {
325       Error("Next","Error occured while loading kinematics.");
326       delete fRunLoader;
327       delete ret;
328       return 0x0;
329     }
330   }
331  return ret;
332 }
333 /**********************************************************/
334
335 Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
336 {
337   //returns pdg code from the PID index
338   //ask jura about charge
339   switch (spec)
340    {
341      case kESDElectron:
342        return kPositron;
343        break;
344      case kESDMuon:
345        return kMuonPlus;
346        break;
347      case kESDPion:
348        return kPiPlus;
349        break;
350      case kESDKaon:
351        return kKPlus;
352        break;
353      case kESDProton:
354        return kProton;
355        break;
356      default:
357        ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
358        break;
359    }
360   return 0;
361 }