Tracking in non-uniform nmagnetic field (Yu.Belikov)
[u/mrichter/AliRoot.git] / ANALYSIS / AliReaderESD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //____________________________________________________________________
19 //////////////////////////////////////////////////////////////////////
20 //                                                                  //
21 // class AliReaderESD                                               //
22 //                                                                  //
23 // reader for ALICE Event Summary Data (ESD).                       //
24 //                                                                  //
25 // Piotr.Skowronski@cern.ch                                         //
26 //                                                                  //
27 //////////////////////////////////////////////////////////////////////
28
29 #include <TFile.h>
30 #include <TGliteXmlEventlist.h>
31 #include <TH1.h>
32 #include <TKey.h>
33 #include <TObjString.h>
34 #include <TPDGCode.h>
35 #include <TParticle.h>
36 #include <TString.h>
37 #include <TTree.h>
38
39 #include "AliAnalysis.h"
40 #include "AliAODRun.h"
41
42 #include "AliAOD.h"
43 #include "AliAODParticle.h"
44 #include "AliAODParticleCut.h"
45 #include "AliAODRun.h"
46 #include "AliAnalysis.h"
47 #include "AliClusterMap.h"
48 #include "AliESD.h"
49 #include "AliESDtrack.h"
50 #include "AliKalmanTrack.h"
51 #include "AliLog.h"
52 #include "AliReaderESD.h"
53 #include "AliRun.h"
54 #include "AliRunLoader.h"
55 #include "AliStack.h"
56
57 ClassImp(AliReaderESD)
58
59 AliReaderESD::AliReaderESD(const Char_t* esdfilename, const Char_t* galfilename):
60  fESDFileName(esdfilename),
61  fGAlFileName(galfilename),
62  fFile(0x0),
63  fRunLoader(0x0),
64  fKeyIterator(0x0),
65  fReadSim(kFALSE),
66  fCheckParticlePID(kFALSE),
67  fReadMostProbableOnly(kFALSE),
68  fNTrackPoints(0),
69  fdR(0.0),
70  fClusterMap(kFALSE),
71  fITSTrackPoints(kFALSE),
72  fITSTrackPointsType(AliTrackPoints::kITS),
73  fMustTPC(kFALSE),
74  fReadCentralBarrel(kTRUE),
75  fReadMuon(kFALSE),
76  fReadPHOS(kFALSE),
77  fNTPCClustMin(0),
78  fNTPCClustMax(1500),
79  fTPCChi2PerClustMin(0.0),
80  fTPCChi2PerClustMax(10e5),
81  fChi2Min(0.0),
82  fChi2Max(10e5),
83  fC00Min(0.0),
84  fC00Max(10e5),
85  fC11Min(0.0),
86  fC11Max(10e5),
87  fC22Min(0.0),
88  fC22Max(10e5),
89  fC33Min(0.0),
90  fC33Max(10e5),
91  fC44Min(0.0),
92  fC44Max(10e5),
93  fTPCC00Min(0.0),
94  fTPCC00Max(10e5),
95  fTPCC11Min(0.0),
96  fTPCC11Max(10e5),
97  fTPCC22Min(0.0),
98  fTPCC22Max(10e5),
99  fTPCC33Min(0.0),
100  fTPCC33Max(10e5),
101  fTPCC44Min(0.0),
102  fTPCC44Max(10e5)
103
104 {
105   //cosntructor
106   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
107     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
108 }
109 /********************************************************************/
110   
111 AliReaderESD::AliReaderESD(TObjArray* dirs,const Char_t* esdfilename, const Char_t* galfilename):
112  AliReader(dirs), 
113  fESDFileName(esdfilename),
114  fGAlFileName(galfilename),
115  fFile(0x0),
116  fRunLoader(0x0),
117  fKeyIterator(0x0),
118  fReadSim(kFALSE),
119  fCheckParticlePID(kFALSE),
120  fReadMostProbableOnly(kFALSE),
121  fNTrackPoints(0),
122  fdR(0.0),
123  fClusterMap(kFALSE),
124  fITSTrackPoints(kFALSE),
125  fITSTrackPointsType(AliTrackPoints::kITS),
126  fMustTPC(kFALSE),
127  fReadCentralBarrel(kTRUE),
128  fReadMuon(kFALSE),
129  fReadPHOS(kFALSE),
130  fNTPCClustMin(0),
131  fNTPCClustMax(150),
132  fTPCChi2PerClustMin(0.0),
133  fTPCChi2PerClustMax(10e5),
134  fChi2Min(0.0),
135  fChi2Max(10e5),
136  fC00Min(0.0),
137  fC00Max(10e5),
138  fC11Min(0.0),
139  fC11Max(10e5),
140  fC22Min(0.0),
141  fC22Max(10e5),
142  fC33Min(0.0),
143  fC33Max(10e5),
144  fC44Min(0.0),
145  fC44Max(10e5),
146  fTPCC00Min(0.0),
147  fTPCC00Max(10e5),
148  fTPCC11Min(0.0),
149  fTPCC11Max(10e5),
150  fTPCC22Min(0.0),
151  fTPCC22Max(10e5),
152  fTPCC33Min(0.0),
153  fTPCC33Max(10e5),
154  fTPCC44Min(0.0),
155  fTPCC44Max(10e5)
156 {
157   //cosntructor
158   if ( ((Int_t)kNSpecies) != ((Int_t)AliPID::kSPECIES))
159     Fatal("AliReaderESD","ESD defintions probobly changed. Ask Youra.");
160 }
161 /********************************************************************/
162
163 AliReaderESD::~AliReaderESD()
164 {
165  //desctructor
166     delete fRunLoader;
167     delete fKeyIterator;
168     delete fFile;
169 }
170
171 /**********************************************************/
172 Int_t AliReaderESD::ReadNext()
173 {
174 //reads next event from fFile
175   //fRunLoader is for reading Kine
176   
177   AliDebug(1,"Entered");
178     
179   if (fEventSim == 0x0)  fEventSim = new AliAOD();
180   if (fEventRec == 0x0)  fEventRec = new AliAOD();
181   
182   fEventSim->Reset();
183   fEventRec->Reset();
184         
185   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
186     {
187       if (fFile == 0x0)
188        {
189          fFile = OpenFile(fCurrentDir);//rl is opened here
190          if (fFile == 0x0)
191            {
192              Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
193              fCurrentDir++;
194              continue;
195            }
196          fCurrentEvent = 0;
197        }
198      TString esdname = "ESD";
199      esdname+=fCurrentEvent;
200      AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
201      if (esd == 0x0)
202       {
203         AliDebug(3,Form("Can not find AliESD object named %s",esdname.Data()));
204         fCurrentDir++;
205         delete fFile;//we have to assume there is no more ESD objects in the fFile
206         fFile = 0x0;
207         delete fRunLoader;
208         fRunLoader = 0x0;
209         continue;
210       }
211       ReadESD(esd);
212       
213       fCurrentEvent++;
214       fNEventsRead++;
215       delete esd;
216       return 0;//success -> read one event
217     }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array  
218    
219   return 1; //no more directories to read
220 }
221
222 /**********************************************************/
223 Int_t AliReaderESD::ReadESD(AliESD* esd)
224 {
225 //Reads esd data
226   if (esd == 0x0)
227    {
228      Error("ReadESD","ESD is NULL");
229      return 1;
230    }
231   
232   // seperate each method
233   if (fReadCentralBarrel) ReadESDCentral(esd);
234
235   if (fReadMuon) ReadESDMuon(esd);
236
237   if (fReadPHOS) ReadESDPHOS(esd);
238
239   return 1;
240 }
241
242 /**********************************************************/
243 Int_t AliReaderESD::ReadESDCentral(AliESD* esd)
244 {
245   //****** Tentative particle type "concentrations"
246   static const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
247   
248   Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
249   Double_t w[kNSpecies];
250   Double_t mom[3];//momentum
251   Double_t pos[3];//position
252   Double_t vertexpos[3];//vertex position
253   //Reads one ESD
254
255   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
256   if (pdgdb == 0x0)
257    {
258      Error("ReadESD","Can not get PDG Database Instance.");
259      return 1;
260    }
261    
262   Float_t mf = esd->GetMagneticField()/10.; //AliESD::GetMagnField returns mf in kG
263
264   if ( (mf == 0.0) && ((fNTrackPoints > 0) || fITSTrackPoints) )
265    {
266       Error("ReadESD","Magnetic Field is 0 and Track Points Demended. Skipping to next event.");
267
268    }
269
270   if (fITSTrackPoints)
271    {
272      Info("ReadESD","Magnetic Field is %f",mf);
273      //AliKalmanTrack::SetMagneticField(mf);
274    }
275  
276   AliStack* stack = 0x0;
277   if (fReadSim && fRunLoader)
278    {
279      fRunLoader->GetEvent(fCurrentEvent);
280      stack = fRunLoader->Stack();
281    }
282
283   const AliESDVertex* vertex = esd->GetVertex();
284   if (vertex == 0x0)
285    {
286      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
287      vertexpos[0] = 0.0;
288      vertexpos[1] = 0.0;
289      vertexpos[2] = 0.0;
290    }
291   else
292    {
293      vertex->GetXYZ(vertexpos);
294    }
295    
296   AliDebug(1,Form("Primary Vertex is (%f,%f,%f)",vertexpos[0],vertexpos[1],vertexpos[2]));
297    
298   Info("ReadESD","Reading Event %d",fCurrentEvent);
299
300   Int_t ntr = esd->GetNumberOfTracks();
301   Info("ReadESD","Found %d tracks.",ntr);
302   for (Int_t i = 0;i<ntr; i++)
303    {
304      AliESDtrack *esdtrack = esd->GetTrack(i);
305      if (esdtrack == 0x0)
306       {
307         Error("Next","Can not get track %d", i);
308         continue;
309       }
310
311      //if (esdtrack->HasVertexParameters() == kFALSE) 
312      if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
313       {
314         AliDebug(3,Form("Particle skipped: Data at vertex not available."));
315         continue;
316       }
317
318      if (fMustTPC)
319       {
320        if ((esdtrack->GetStatus() & AliESDtrack::kTPCin) == kFALSE)
321         {
322           AliDebug(3,"Particle skipped: Was not reconstructed in TPC.");
323           continue;
324         }
325       }     
326
327      if ((esdtrack->GetStatus() & AliESDtrack::kESDpid) == kFALSE) 
328       {
329         AliDebug(3,"Particle skipped: PID BIT is not set.");
330         continue;
331       }
332
333
334      Double_t extx;
335      Double_t extp[5];
336      esdtrack->GetConstrainedExternalParameters(extx,extp);
337      if (extp[4] == 0.0)
338       {
339         AliDebug(3,"Track has 0 contrianed curvature -> Probobly parameters never updated. Skipping.");
340         continue;
341       } 
342      esdtrack->GetESDpid(pidtable);
343      esdtrack->GetConstrainedPxPyPz(mom);
344      esdtrack->GetConstrainedXYZ(pos);
345      pos[0] -= vertexpos[0];//we are interested only in relative position to Primary vertex at this point
346      pos[1] -= vertexpos[1];
347      pos[2] -= vertexpos[2];
348
349      Int_t charge = (extp[4] > 0)?1:-1;//if curvature=charg/Pt is positive charge is positive
350
351      //Particle from kinematics
352      AliAODParticle* particle = 0;
353      Bool_t keeppart = kFALSE;
354      if ( fReadSim && stack  )
355       {
356         if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track 
357         TParticle *p = stack->Particle(esdtrack->GetLabel());
358         if (p==0x0) 
359          {
360            Error("ReadNext","Can not find track with such label.");
361            continue;
362          }
363          
364         if (fCheckParticlePID)
365          {
366            if(Rejected(p->GetPdgCode())) 
367             {
368               AliDebug(6,Form("Simulated Particle PID (%d) did not pass the cut.",p->GetPdgCode()));
369               continue; //check if we are intersted with particles of this type 
370             }
371          }  
372 //           if(p->GetPdgCode()<0) charge = -1;
373         particle = new AliAODParticle(*p,i);
374
375       }
376       
377      if(CheckTrack(esdtrack)) continue;
378       
379      //Here we apply Bayes' formula
380      Double_t rc=0.;
381      for (Int_t s=0; s<AliPID::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
382      if (rc==0.0) 
383       {
384         AliDebug(3,"Particle rejected since total bayessian PID probab. is zero.");
385         continue;
386       }
387
388      for (Int_t s=0; s<AliPID::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
389
390      if (AliDebugLevel() > 4)
391       { 
392         Info("ReadNext","###########################################################################");
393         Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
394         Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
395         TString msg("Pid list got from track:");
396         for (Int_t s = 0;s<kNSpecies;s++)
397          {
398            msg+="\n    ";
399            msg+=s;
400            msg+="(";
401            msg+=charge*GetSpeciesPdgCode((ESpecies)s);
402            msg+="): ";
403            msg+=w[s];
404            msg+=" (";
405            msg+=pidtable[s];
406            msg+=")";
407          }
408         Info("ReadNext","%s",msg.Data());
409       }//if (AliDebugLevel()>4)
410
411       AliTrackPoints* tpts = 0x0;
412       if (fNTrackPoints > 0) 
413        {
414          tpts = new AliTrackPoints(fNTrackPoints,esdtrack,mf*10.0,fdR);
415 //         tpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
416        }
417
418       AliTrackPoints* itstpts = 0x0;
419       if (fITSTrackPoints) 
420        {
421          itstpts = new AliTrackPoints(fITSTrackPointsType,esdtrack,mf*10.0);
422 //         itstpts->Move(-vertexpos[0],-vertexpos[1],-vertexpos[2]);
423        }
424
425       AliClusterMap* cmap = 0x0; 
426       if (  fClusterMap ) 
427        {
428          cmap = new AliClusterMap(esdtrack);
429        }
430      
431      //If this flag fReadMostProbableOnly is false the 
432      //loop over species (see "LOOPPIDS") is over all possible PIDs
433      //in other case the most probablle PID is searched
434      //and the loop is limited to that PID only
435
436      Int_t firstspecie = 0;
437      Int_t lastspecie = kNSpecies;
438      
439      if (fReadMostProbableOnly)
440       { 
441         //find the most probable PID
442         Int_t spec = 0;
443         Float_t maxprob = w[0];
444         for (Int_t s=1; s<AliPID::kSPECIES; s++) 
445          {
446            if (w[s]>maxprob)
447             {
448               maxprob = w[s];
449               spec = s;
450             }
451          }
452         firstspecie = spec;
453         lastspecie = spec + 1;
454       }
455      
456      for (Int_t s = firstspecie; s<lastspecie; s++)//LOOPPIDS
457       {
458         Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
459         Float_t pp = w[s];
460         if (pp == 0.0) 
461          {
462            AliDebug(6,Form("Probability of being PID %d is zero. Continuing.",pdgcode));
463            continue;
464          }
465
466         if(Rejected(pdgcode)) 
467          {
468            AliDebug(6,Form("PID (%d) did not pass the cut.",pdgcode));
469            continue; //check if we are intersted with particles of this type 
470          }
471
472         Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
473         Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
474
475         AliAODParticle* track = new AliAODParticle(pdgcode, w[s],i, 
476                                                    mom[0], mom[1], mom[2], tEtot,
477                                                    pos[0], pos[1], pos[2], 0.);
478         //copy probabilitis of other species (if not zero)
479         for (Int_t k = 0; k<kNSpecies; k++)
480          {
481            if (k == s) continue;
482            if (w[k] == 0.0) continue;
483            track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
484          }
485
486         if(Rejected(track))//check if meets all criteria of any of our cuts
487                        //if it does not delete it and take next good track
488          { 
489            AliDebug(5,"Track did not pass the cut");
490            delete track;
491            continue;
492          }
493
494         //Single Particle cuts on cluster map and track points rather do not have sense
495         if (tpts)
496          {
497            track->SetTPCTrackPoints(tpts);
498          }
499          
500         if (itstpts)
501          {
502            track->SetITSTrackPoints(itstpts); 
503          }
504
505         if (cmap) 
506          { 
507            track->SetClusterMap(cmap);
508          }
509
510         fEventRec->AddParticle(track);
511         if (particle) fEventSim->AddParticle(particle);
512         keeppart = kTRUE;
513
514         if (AliDebugLevel() > 4 )
515          {
516            Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
517            track->Print();
518            if (particle) particle->Print();
519            Info("ReadNext","\n----------------------------------------------\n");
520          }
521       }//for (Int_t s = 0; s<kNSpecies; s++)
522
523      if (keeppart == kFALSE) 
524       {
525         delete particle;//particle was not stored in event
526         delete tpts;
527         delete itstpts;
528         delete cmap;
529       }
530      else
531       {
532         if ( fReadSim && stack  )
533          {
534            if (particle->P() < 0.00001)
535             {
536               Info("ReadNext","###################################");
537               Info("ReadNext","###################################");
538               Info("ReadNext","Track Label %d",esdtrack->GetLabel());
539               TParticle *p = stack->Particle(esdtrack->GetLabel());
540               Info("ReadNext","");
541               p->Print();
542               Info("ReadNext","");
543               particle->Print();
544             }
545          }   
546       } 
547
548    }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
549
550   Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
551          fEventRec->GetNumberOfParticles(), fEventSim->GetNumberOfParticles(),
552          fNEventsRead,fCurrentEvent,fCurrentDir);
553   fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
554   
555   /******************************************************/
556   /******     Setting glevet properties     *************/
557   /******************************************************/
558     
559   if (fEventRec->GetNumberOfParticles() > 0)
560    {
561      fEventRec->SetPrimaryVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
562    }
563   return 0;
564 }
565
566 /**********************************************************/
567 Int_t AliReaderESD::ReadESDMuon(AliESD* esd)
568 {
569
570   Double_t vertexpos[3];//vertex position, assuming no secondary decay
571
572   const AliESDVertex* vertex = esd->GetVertex();
573
574   if (vertex == 0x0) {
575     Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
576     vertexpos[0] = 0.0;
577     vertexpos[1] = 0.0;
578     vertexpos[2] = 0.0;
579   } else {
580     vertex->GetXYZ(vertexpos);
581   }
582
583  Int_t nTracks = (Int_t)esd->GetNumberOfMuonTracks() ;
584
585  AliDebug(1,Form("Reading Event %d \nFound %d tracks.",fCurrentEvent,nTracks));
586
587  // settings
588   Float_t Chi2Cut = 100.;
589   Float_t PtCutMin = 1.;
590   Float_t PtCutMax = 10000.;
591   Float_t muonMass = 0.105658389;
592   Int_t pdgcode = -13;
593   Double_t thetaX, thetaY, pYZ;
594   Double_t pxRec1, pyRec1, pzRec1, E1;
595   Int_t charge;
596
597   Int_t ntrackhits;
598   Double_t fitfmin;
599
600   TLorentzVector fV1;
601   fEventRec->Reset();
602   for (Int_t iTrack = 0; iTrack <  nTracks;  iTrack++) {
603
604       AliESDMuonTrack* muonTrack = esd->GetMuonTrack(iTrack);
605
606       thetaX = muonTrack->GetThetaX();
607       thetaY = muonTrack->GetThetaY();
608
609       pYZ     =  1./TMath::Abs(muonTrack->GetInverseBendingMomentum());
610       pzRec1  = - pYZ / TMath::Sqrt(1.0 + TMath::Tan(thetaY)*TMath::Tan(thetaX));
611       pxRec1  = pzRec1 * TMath::Tan(thetaX);
612       pyRec1  = pzRec1 * TMath::Tan(thetaY);
613       charge = Int_t(TMath::Sign(1.,muonTrack->GetInverseBendingMomentum()));
614       E1 = TMath::Sqrt(muonMass * muonMass + pxRec1 * pxRec1 + pyRec1 * pyRec1 + pzRec1 * pzRec1);
615       fV1.SetPxPyPzE(pxRec1, pyRec1, pzRec1, E1);
616
617       ntrackhits = muonTrack->GetNHit();
618       fitfmin    = muonTrack->GetChi2();
619
620       // transverse momentum
621       Float_t pt1 = fV1.Pt();
622
623       // chi2 per d.o.f.
624       Float_t ch1 =  fitfmin / (2.0 * ntrackhits - 5);
625
626       if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
627         AliAODParticle* track = new AliAODParticle(pdgcode*charge,1,iTrack, 
628                                                    pxRec1, pyRec1,pzRec1, E1,
629                                                    vertexpos[0], vertexpos[1], vertexpos[2], 0.);
630         fEventRec->AddParticle(track);
631       }
632  
633   }
634   fTrackCounter->Fill(fEventRec->GetNumberOfParticles());
635   return 0;
636 }
637
638 /**********************************************************/
639
640 void AliReaderESD::Rewind()
641 {
642   //rewinds reading 
643   //  delete fKeyIterator;
644   delete fFile;
645   fFile = 0x0;
646   // fKeyIterator = 0x0;
647   delete fRunLoader;
648   fRunLoader = 0x0;
649   fCurrentDir = 0;
650   fNEventsRead = 0;
651   if (fEventList) fEventList->Reset(); 
652   if (fTrackCounter) fTrackCounter->Reset();
653 }
654 /**********************************************************/
655
656 TFile* AliReaderESD::OpenFile(Int_t n)
657 {
658 //opens fFile with  tree
659  if (fEventList)
660   { 
661     if (fCurrentDir > n)
662      {
663        fEventList->Reset();
664        fCurrentDir = 0;
665      }
666     
667     while (fCurrentDir < n)
668      {
669        fEventList->Next();
670        fCurrentDir++;
671      }
672     fEventList->Next();
673   }
674  
675  const TString& dirname = GetDirName(n);
676  if (dirname == "")
677   {
678    Error("OpenFiles","Can not get directory name");
679    return 0x0;
680   }
681  TString filename;
682  
683  if (fEventList)
684   {
685     filename = fEventList->GetURL(fESDFileName);
686   }
687  else
688   {
689     filename = dirname +"/"+ fESDFileName;
690   }
691  Info("OpenFile","%s ==> %s",fESDFileName.Data(),filename.Data()); 
692  TFile *ret = TFile::Open(filename.Data()); 
693
694  if ( ret == 0x0)
695   {
696     Error("OpenFiles","Can't open fFile %s",filename.Data());
697     return 0x0;
698   }
699  if (!ret->IsOpen())
700   {
701     Error("OpenFiles","Can't open fFile  %s",filename.Data());
702     return 0x0;
703   }
704  
705  if (fReadSim )
706   {
707     TString gafilename;
708     if (fEventList)
709      {
710        gafilename = fEventList->GetURL(fGAlFileName);
711      }
712     else
713      {
714        gafilename = dirname +"/"+ fGAlFileName;
715      }
716    Info("OpenFile","%s ==> %s",fGAlFileName.Data(),gafilename.Data()); 
717
718    fRunLoader = AliRunLoader::Open(gafilename);
719    
720    if (fRunLoader == 0x0)
721     {
722       Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
723       delete ret;
724       return 0x0;
725     }
726     
727    fRunLoader->LoadHeader();
728    
729    if (fEventList)
730     {
731       TString kinefilename = fEventList->GetURL("Kinematics.root");
732       fRunLoader->SetKineFileName(kinefilename);
733       Info("OpenFile","%s ==> %s","Kinematics.root",kinefilename.Data()); 
734     } 
735    
736    if (fRunLoader->LoadKinematics())
737     {
738       Error("Next","Error occured while loading kinematics.");
739       delete fRunLoader;
740       delete ret;
741       return 0x0;
742     }
743   }
744    
745  return ret;
746 }
747 /**********************************************************/
748
749 Int_t AliReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
750 {
751   //returns pdg code from the PID index
752   //ask jura about charge
753   switch (spec)
754    {
755      case kESDElectron:
756        return kPositron;
757        break;
758      case kESDMuon:
759        return kMuonPlus;
760        break;
761      case kESDPion:
762        return kPiPlus;
763        break;
764      case kESDKaon:
765        return kKPlus;
766        break;
767      case kESDProton:
768        return kProton;
769        break;
770      default:
771        ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
772        break;
773    }
774   return 0;
775 }
776 /********************************************************************/
777 Bool_t AliReaderESD::CheckTrack(AliESDtrack* t) const
778 {
779   //Performs check of the track
780   
781   if ( (t->GetConstrainedChi2() < fChi2Min) || (t->GetConstrainedChi2() > fChi2Max) ) return kTRUE;
782   
783   if ( (t->GetTPCclusters(0x0) < fNTPCClustMin) || (t->GetTPCclusters(0x0) > fNTPCClustMax) ) return kTRUE;
784
785   if (t->GetTPCclusters(0x0) > 0)
786    {
787      Float_t chisqpercl = t->GetTPCchi2()/((Double_t)t->GetTPCclusters(0x0));
788      if ( (chisqpercl < fTPCChi2PerClustMin) || (chisqpercl > fTPCChi2PerClustMax) ) return kTRUE;
789    }
790
791   Double_t cc[15];
792   t->GetConstrainedExternalCovariance(cc);
793
794   if ( (cc[0]  < fC00Min) || (cc[0]  > fC00Max) ) return kTRUE;
795   if ( (cc[2]  < fC11Min) || (cc[2]  > fC11Max) ) return kTRUE;
796   if ( (cc[5]  < fC22Min) || (cc[5]  > fC22Max) ) return kTRUE;
797   if ( (cc[9]  < fC33Min) || (cc[9]  > fC33Max) ) return kTRUE;
798   if ( (cc[14] < fC44Min) || (cc[14] > fC44Max) ) return kTRUE;
799
800
801   t->GetInnerExternalCovariance(cc);
802
803   if ( (cc[0]  < fTPCC00Min) || (cc[0]  > fTPCC00Max) ) return kTRUE;
804   if ( (cc[2]  < fTPCC11Min) || (cc[2]  > fTPCC11Max) ) return kTRUE;
805   if ( (cc[5]  < fTPCC22Min) || (cc[5]  > fTPCC22Max) ) return kTRUE;
806   if ( (cc[9]  < fTPCC33Min) || (cc[9]  > fTPCC33Max) ) return kTRUE;
807   if ( (cc[14] < fTPCC44Min) || (cc[14] > fTPCC44Max) ) return kTRUE;
808
809   return kFALSE;
810
811 }
812 /********************************************************************/
813
814 void AliReaderESD::SetChi2Range(Float_t min, Float_t max)
815 {
816   //sets range of Chi2 per Cluster
817   fChi2Min = min;
818   fChi2Max = max;
819 }
820 /********************************************************************/
821
822 void AliReaderESD::SetTPCNClustersRange(Int_t min,Int_t max)
823 {
824  //sets range of Number Of Clusters that tracks have to have
825  fNTPCClustMin = min;
826  fNTPCClustMax = max;
827 }
828 /********************************************************************/
829
830 void AliReaderESD::SetTPCChi2PerCluserRange(Float_t min, Float_t max)
831 {
832   //sets range of Chi2 per Cluster
833   fTPCChi2PerClustMin = min;
834   fTPCChi2PerClustMax = max;
835 }
836 /********************************************************************/
837
838 void AliReaderESD::SetC00Range(Float_t min, Float_t max)
839 {
840  //Sets range of C00 parameter of covariance matrix of the track
841  //it defines uncertainty of the momentum
842   fC00Min = min;
843   fC00Max = max;
844 }
845 /********************************************************************/
846
847 void AliReaderESD::SetC11Range(Float_t min, Float_t max)
848 {
849  //Sets range of C11 parameter of covariance matrix of the track
850  //it defines uncertainty of the momentum
851   fC11Min = min;
852   fC11Max = max;
853 }
854 /********************************************************************/
855
856 void AliReaderESD::SetC22Range(Float_t min, Float_t max)
857 {
858  //Sets range of C22 parameter of covariance matrix of the track
859  //it defines uncertainty of the momentum
860   fC22Min = min;
861   fC22Max = max;
862 }
863 /********************************************************************/
864
865 void AliReaderESD::SetC33Range(Float_t min, Float_t max)
866 {
867  //Sets range of C33 parameter of covariance matrix of the track
868  //it defines uncertainty of the momentum
869   fC33Min = min;
870   fC33Max = max;
871 }
872 /********************************************************************/
873
874 void AliReaderESD::SetC44Range(Float_t min, Float_t max)
875 {
876  //Sets range of C44 parameter of covariance matrix of the track
877  //it defines uncertainty of the momentum
878   fC44Min = min;
879   fC44Max = max;
880 }