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