HBTAN/AliHBTReaderESD.cxx

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