HBTAN/AliHBTReaderESD.cxx

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