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 <AliESDtrack.h>
  11 #include <AliESD.h>
  12
  13 #include "AliHBTRun.h"
  14 #include "AliHBTEvent.h"
  15 #include "AliHBTParticle.h"
  16 #include "AliHBTParticleCut.h"
  17
  18 ClassImp(AliHBTReaderESD)
  19
  20 AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
  21  fParticles(new AliHBTRun()),
  22  fTracks(new AliHBTRun()),
  23  fESDFileName(esdfilename),
  24  fIsRead(kFALSE)
  25 {
  26   //cosntructor
  27   if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
  28     Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
  29 }
  30 /********************************************************************/
  31
  32 AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
  33  AliHBTReader(dirs),
  34  fParticles(new AliHBTRun()),
  35  fTracks(new AliHBTRun()),
  36  fESDFileName(esdfilename),
  37  fIsRead(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()
  46 {
  47  //desctructor
  48   delete fParticles;
  49   delete fTracks;
  50 }
  51 /********************************************************************/
  52
  53 AliHBTEvent* AliHBTReaderESD::GetParticleEvent(Int_t n)
  54  {
  55  //returns Nth event with simulated particles
  56    if (!fIsRead)
  57     if(Read(fParticles,fTracks))
  58      {
  59        Error("GetParticleEvent","Error in reading");
  60        return 0x0;
  61      }
  62    return fParticles->GetEvent(n);
  63  }
  64 /********************************************************************/
  65
  66 AliHBTEvent* AliHBTReaderESD::GetTrackEvent(Int_t n)
  67  {
  68  //returns Nth event with reconstructed tracks
  69    if (!fIsRead)
  70     if(Read(fParticles,fTracks))
  71      {
  72        Error("GetTrackEvent","Error in reading");
  73        return 0x0;
  74      }
  75    return fTracks->GetEvent(n);
  76  }
  77 /********************************************************************/
  78
  79 Int_t AliHBTReaderESD::GetNumberOfPartEvents()
  80  {
  81  //returns number of events of particles
  82    if (!fIsRead)
  83     if ( Read(fParticles,fTracks))
  84      {
  85        Error("GetNumberOfPartEvents","Error in reading");
  86        return 0;
  87      }
  88    return fParticles->GetNumberOfEvents();
  89  }
  90
  91 /********************************************************************/
  92 Int_t AliHBTReaderESD::GetNumberOfTrackEvents()
  93  {
  94  //returns number of events of tracks
  95   if (!fIsRead)
  96     if(Read(fParticles,fTracks))
  97      {
  98        Error("GetNumberOfTrackEvents","Error in reading");
  99        return 0;
 100      }
 101   return fTracks->GetNumberOfEvents();
 102  }
 103 /********************************************************************/
 104
 105
 106 Int_t AliHBTReaderESD::Read(AliHBTRun* particles, AliHBTRun *tracks)
 107 {
 108  //reads data and puts put to the particles and tracks objects
 109  //reurns 0 if everything is OK
 110  //
 111   Info("Read","");
 112   Int_t totalNevents = 0;
 113   Int_t currentdir = 0; //number of current directory name is fDirs array
 114
 115   Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
 116   Double_t w[kNSpecies];
 117   Double_t mom[3];//momentum
 118   Double_t pos[3];//position
 119    //****** Tentative particle type "concentrations"
 120   const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
 121
 122   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
 123   if (pdgdb == 0x0)
 124    {
 125      Error("Read","Can not get PDG Database Instance.");
 126      return 1;
 127    }
 128   if (!particles) //check if an object is instatiated
 129    {
 130      Error("Read"," particles object must instatiated before passing it to the reader");
 131      return 1;
 132    }
 133   if (!tracks)  //check if an object is instatiated
 134    {
 135      Error("Read"," tracks object must instatiated before passing it to the reader");
 136      return 1;
 137    }
 138   particles->Reset();//clear runs == delete all old events
 139   tracks->Reset();
 140
 141   Int_t Ndirs;
 142   if (fDirs) //if array with directories is supplied by user
 143    {
 144      Ndirs = fDirs->GetEntries(); //get the number if directories
 145    }
 146   else
 147    {
 148      Ndirs = 0; //if the array is not supplied read only from current directory
 149    }
 150
 151   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 152    {
 153      TFile* file = OpenFile(currentdir);
 154      if (file == 0x0)
 155        {
 156          Error("Read","Cannot get File for dir no. %d",currentdir);
 157          currentdir++;
 158          continue;
 159        }
 160
 161      for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events
 162       {
 163        TString esdname;
 164        esdname+=currentEvent;
 165        AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname));
 166        if (esd == 0x0)
 167          {
 168            if (AliHBTParticle::fgDebug > 2 )
 169             {
 170               Info("Read","Can not find ESD object named %s.",esdname.Data());
 171             }
 172            currentdir++;
 173            break;//we have to assume there is no more ESD objects in the file
 174          }
 175
 176        Info("Read","Reading Event %d",currentEvent);
 177
 178        Int_t ntr = esd->GetNumberOfTracks();
 179        Info("Read","Found %d tracks.",ntr);
 180        for (Int_t i = 0;i<ntr; i++)
 181         {
 182           AliESDtrack *esdtrack = esd->GetTrack(i);
 183           if (esdtrack == 0x0)
 184            {
 185              Error("Read","Can not get track %d", i);
 186              continue;
 187            }
 188           if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
 189            {
 190              if (AliHBTParticle::fgDebug > 2)
 191                Info("Read","Particle skipped: PID BIT is not set.");
 192              continue;
 193            }
 194
 195           esdtrack->GetESDpid(pidtable);
 196           esdtrack->GetPxPyPz(mom);
 197           esdtrack->GetXYZ(pos);
 198           Double_t rc=0.;
 199
 200            //Here we apply Bayes' formula
 201           for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
 202           if (rc==0.0)
 203            {
 204              if (AliHBTParticle::fgDebug > 2)
 205                Info("Read","Particle rejected since total bayessian PID probab. is zero.");
 206              continue;
 207            }
 208
 209           for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 210
 211           if (AliHBTParticle::fgDebug > 4)
 212            {
 213              Info("Read","###########################################################################");
 214              Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
 215              Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]);
 216              Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc);
 217              TString msg("Pid list got from track:");
 218              for (Int_t s = 0;s<kNSpecies;s++)
 219               {
 220                 msg+="\n    ";
 221                 msg+=s;
 222                 msg+="(";
 223                 msg+=GetSpeciesPdgCode((ESpecies)s);
 224                 msg+="): ";
 225                 msg+=w[s];
 226                 msg+=" (";
 227                 msg+=pidtable[s];
 228                 msg+=")";
 229               }
 230              Info("Read","%s",msg.Data());
 231            }//if (AliHBTParticle::fgDebug>4)
 232
 233           for (Int_t s = 0; s<kNSpecies; s++)
 234            {
 235              if (w[s] == 0.0) continue;
 236
 237              Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
 238              if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
 239
 240              Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
 241              Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
 242
 243              AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
 244                                                         mom[0], mom[1], mom[2], tEtot,
 245                                                         pos[0], pos[1], pos[2], 0.);
 246              //copy probabilitis of other species (if not zero)
 247              for (Int_t k = 0; k<kNSpecies; k++)
 248               {
 249                 if (k == s) continue;
 250                 if (w[k] == 0.0) continue;
 251                 track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
 252               }
 253              if(Pass(track))//check if meets all criteria of any of our cuts
 254                             //if it does not delete it and take next good track
 255                {
 256                  delete track;
 257                  continue;
 258                }
 259              tracks->AddParticle(totalNevents,track);
 260              if (AliHBTParticle::fgDebug > 4 )
 261               {
 262                 Info("Read","\n\nAdding Particle with incarnation %d",pdgcode);
 263                 track->Print();
 264               }
 265            }//for (Int_t s = 0; s<kNSpecies; s++)
 266
 267         }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
 268        totalNevents++;
 269       }//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file
 270
 271    }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 272
 273   fIsRead = kTRUE;
 274   return 0;
 275 }
 276 /**********************************************************/
 277
 278 TFile* AliHBTReaderESD::OpenFile(Int_t n)
 279 {
 280 //opens file with kine tree
 281
 282  const TString& dirname = GetDirName(n);
 283  if (dirname == "")
 284   {
 285    Error("OpenFiles","Can not get directory name");
 286    return 0x0;
 287   }
 288  TString filename = dirname +"/"+ fESDFileName;
 289  TFile *ret = TFile::Open(filename.Data());
 290
 291  if ( ret == 0x0)
 292   {
 293     Error("OpenFiles","Can't open file %s",filename.Data());
 294     return 0x0;
 295   }
 296  if (!ret->IsOpen())
 297   {
 298     Error("OpenFiles","Can't open file  %s",filename.Data());
 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 }