03f18847f9c161374c8af15d113fbc416ee755e8
[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          {
146            if (AliHBTParticle::fgDebug > 2) 
147              Info("ReadNext","Particle skipped: Data at vertex not available.");
148            continue;
149          }
150          
151         if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
152          {
153            if (AliHBTParticle::fgDebug > 2) 
154              Info("ReadNext","Particle skipped: PID BIT is not set.");
155            continue;
156          }
157        
158         esdtrack->GetESDpid(pidtable);
159         esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]);
160         esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
161         
162         Double_t extx;
163         Double_t extp[5];
164         esdtrack->GetExternalParameters(extx,extp);
165         
166         Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
167         
168         //Particle from kinematics
169         AliHBTParticle* particle = 0;
170         Bool_t keeppart = kFALSE;
171         if ( fReadParticles && stack  )
172          {
173            if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
174            TParticle *p = stack->Particle(esdtrack->GetLabel());
175            if (p==0x0) 
176             {
177               Error("ReadNext","Can not find track with such label.");
178               continue;
179             }
180 //           if(p->GetPdgCode()<0) charge = -1;
181            particle = new AliHBTParticle(*p,i);
182            
183          }
184          
185         //Here we apply Bayes' formula
186         Double_t rc=0.;
187         for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
188         if (rc==0.0) 
189          {
190            if (AliHBTParticle::fgDebug > 2) 
191              Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
192            continue;
193          }
194
195         for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
196
197         if (AliHBTParticle::fgDebug > 4)
198          { 
199            Info("ReadNext","###########################################################################");
200            Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
201            Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
202            TString msg("Pid list got from track:");
203            for (Int_t s = 0;s<kNSpecies;s++)
204             {
205               msg+="\n    ";
206               msg+=s;
207               msg+="(";
208               msg+=charge*GetSpeciesPdgCode((ESpecies)s);
209               msg+="): ";
210               msg+=w[s];
211               msg+=" (";
212               msg+=pidtable[s];
213               msg+=")";
214             }
215            Info("ReadNext","%s",msg.Data());
216          }//if (AliHBTParticle::fgDebug>4)
217
218         for (Int_t s = 0; s<kNSpecies; s++)
219          {
220            Float_t pp = w[s];
221            if (pp == 0.0) continue;
222
223            Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
224            if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 
225
226            Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
227            Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
228
229            AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
230                                                       mom[0], mom[1], mom[2], tEtot,
231                                                       pos[0], pos[1], pos[2], 0.);
232            //copy probabilitis of other species (if not zero)
233            for (Int_t k = 0; k<kNSpecies; k++)
234             {
235               if (k == s) continue;
236               if (w[k] == 0.0) continue;
237               track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
238             }
239            if(Pass(track))//check if meets all criteria of any of our cuts
240                           //if it does not delete it and take next good track
241             { 
242               delete track;
243               continue;
244             }
245            
246            fTracksEvent->AddParticle(track);
247            if (particle) fParticlesEvent->AddParticle(particle);
248            keeppart = kTRUE;
249                
250            if (AliHBTParticle::fgDebug > 4 )
251             {
252               Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
253               track->Print();
254             }
255          }//for (Int_t s = 0; s<kNSpecies; s++)
256          
257         if (keeppart == kFALSE) delete particle;//particle was not stored in event
258         
259       }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
260      
261      Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
262             fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
263             fNEventsRead,fCurrentEvent,fCurrentDir);
264       
265      fCurrentEvent++;
266      fNEventsRead++;
267      delete esd;
268      return 0;//success -> read one event
269    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
270    
271   return 1; //no more directories to read
272 }
273 /**********************************************************/
274 void AliHBTReaderESD::Rewind()
275 {
276   //rewinds reading 
277   delete fFile;
278   fFile = 0;
279   delete fRunLoader;
280   fCurrentDir = 0;
281   fNEventsRead = 0;
282   fCurrentEvent++;
283 }
284 /**********************************************************/
285
286 TFile* AliHBTReaderESD::OpenFile(Int_t n)
287 {
288 //opens fFile with kine tree
289
290  const TString& dirname = GetDirName(n);
291  if (dirname == "")
292   {
293    Error("OpenFiles","Can not get directory name");
294    return 0x0;
295   }
296  TString filename = dirname +"/"+ fESDFileName;
297  TFile *ret = TFile::Open(filename.Data()); 
298
299  if ( ret == 0x0)
300   {
301     Error("OpenFiles","Can't open fFile %s",filename.Data());
302     return 0x0;
303   }
304  if (!ret->IsOpen())
305   {
306     Error("OpenFiles","Can't open fFile  %s",filename.Data());
307     return 0x0;
308   }
309  
310  if (fReadParticles)
311   {
312    fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
313    if (fRunLoader == 0x0)
314     {
315       Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
316       delete ret;
317       return 0x0;
318     }
319    fRunLoader->LoadHeader();
320    if (fRunLoader->LoadKinematics())
321     {
322       Error("Next","Error occured while loading kinematics.");
323       delete fRunLoader;
324       delete ret;
325       return 0x0;
326     }
327   }
328  return ret;
329 }
330 /**********************************************************/
331
332 Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
333 {
334   //returns pdg code from the PID index
335   //ask jura about charge
336   switch (spec)
337    {
338      case kESDElectron:
339        return kPositron;
340        break;
341      case kESDMuon:
342        return kMuonPlus;
343        break;
344      case kESDPion:
345        return kPiPlus;
346        break;
347      case kESDKaon:
348        return kKPlus;
349        break;
350      case kESDProton:
351        return kProton;
352        break;
353      default:
354        ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
355        break;
356    }
357   return 0;
358 }