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