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::GetDebug())
  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::GetDebug() > 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         if ((esdtrack->GetStatus() & AliESDtrack::kITSrefit) == kFALSE)
 146          {
 147            if (AliHBTParticle::GetDebug() > 2)
 148              Info("ReadNext","Particle skipped: Data at vertex not available.");
 149            continue;
 150          }
 151
 152         if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE)
 153          {
 154            if (AliHBTParticle::GetDebug() > 2)
 155              Info("ReadNext","Particle skipped: PID BIT is not set.");
 156            continue;
 157          }
 158
 159         esdtrack->GetESDpid(pidtable);
 160         //esdtrack->GetVertexPxPyPz(mom[0],mom[1],mom[2]);
 161         esdtrack->GetPxPyPz(mom);
 162         //esdtrack->GetVertexXYZ(pos[0],pos[1],pos[2]);
 163         esdtrack->GetXYZ(pos);
 164
 165         Double_t extx;
 166         Double_t extp[5];
 167         esdtrack->GetExternalParameters(extx,extp);
 168
 169         Int_t charge = (extp[4] > 0)?-1:1;//if curvature=-charg/Pt is positive charge is negative
 170
 171         //Particle from kinematics
 172         AliHBTParticle* particle = 0;
 173         Bool_t keeppart = kFALSE;
 174         if ( fReadParticles && stack  )
 175          {
 176            if (esdtrack->GetLabel() < 0) continue;//this is fake -  we are not able to match any track
 177            TParticle *p = stack->Particle(esdtrack->GetLabel());
 178            if (p==0x0)
 179             {
 180               Error("ReadNext","Can not find track with such label.");
 181               continue;
 182             }
 183 //           if(p->GetPdgCode()<0) charge = -1;
 184            particle = new AliHBTParticle(*p,i);
 185
 186          }
 187
 188         //Here we apply Bayes' formula
 189         Double_t rc=0.;
 190         for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
 191         if (rc==0.0)
 192          {
 193            if (AliHBTParticle::GetDebug() > 2)
 194              Info("ReadNext","Particle rejected since total bayessian PID probab. is zero.");
 195            continue;
 196          }
 197
 198         for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
 199
 200         if (AliHBTParticle::GetDebug() > 4)
 201          {
 202            Info("ReadNext","###########################################################################");
 203            Info("ReadNext","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
 204            Info("ReadNext","Position: %f %f %f",pos[0],pos[1],pos[2]);
 205            TString msg("Pid list got from track:");
 206            for (Int_t s = 0;s<kNSpecies;s++)
 207             {
 208               msg+="\n    ";
 209               msg+=s;
 210               msg+="(";
 211               msg+=charge*GetSpeciesPdgCode((ESpecies)s);
 212               msg+="): ";
 213               msg+=w[s];
 214               msg+=" (";
 215               msg+=pidtable[s];
 216               msg+=")";
 217             }
 218            Info("ReadNext","%s",msg.Data());
 219          }//if (AliHBTParticle::GetDebug()>4)
 220
 221         for (Int_t s = 0; s<kNSpecies; s++)
 222          {
 223            Float_t pp = w[s];
 224            if (pp == 0.0) continue;
 225
 226            Int_t pdgcode = charge*GetSpeciesPdgCode((ESpecies)s);
 227            if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type
 228
 229            Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
 230            Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
 231
 232            AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i,
 233                                                       mom[0], mom[1], mom[2], tEtot,
 234                                                       pos[0], pos[1], pos[2], 0.);
 235            //copy probabilitis of other species (if not zero)
 236            for (Int_t k = 0; k<kNSpecies; k++)
 237             {
 238               if (k == s) continue;
 239               if (w[k] == 0.0) continue;
 240               track->SetPIDprobability(charge*GetSpeciesPdgCode( (ESpecies)k ),w[k]);
 241             }
 242            if(Pass(track))//check if meets all criteria of any of our cuts
 243                           //if it does not delete it and take next good track
 244             {
 245               delete track;
 246               continue;
 247             }
 248
 249            fTracksEvent->AddParticle(track);
 250            if (particle) fParticlesEvent->AddParticle(particle);
 251            keeppart = kTRUE;
 252
 253            if (AliHBTParticle::GetDebug() > 4 )
 254             {
 255               Info("ReadNext","\n\nAdding Particle with incarnation %d",pdgcode);
 256               track->Print();
 257             }
 258          }//for (Int_t s = 0; s<kNSpecies; s++)
 259
 260         if (keeppart == kFALSE) delete particle;//particle was not stored in event
 261
 262       }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
 263
 264      Info("ReadNext","Read %d tracks and %d particles from event %d (event %d in dir %d).",
 265             fTracksEvent->GetNumberOfParticles(), fParticlesEvent->GetNumberOfParticles(),
 266             fNEventsRead,fCurrentEvent,fCurrentDir);
 267
 268      fCurrentEvent++;
 269      fNEventsRead++;
 270      delete esd;
 271      return 0;//success -> read one event
 272    }while(fCurrentDir < GetNumberOfDirs());//end of loop over directories specified in fDirs Obj Array
 273
 274   return 1; //no more directories to read
 275 }
 276 /**********************************************************/
 277 void AliHBTReaderESD::Rewind()
 278 {
 279   //rewinds reading
 280   delete fFile;
 281   fFile = 0;
 282   delete fRunLoader;
 283   fCurrentDir = 0;
 284   fNEventsRead = 0;
 285   fCurrentEvent++;
 286 }
 287 /**********************************************************/
 288
 289 TFile* AliHBTReaderESD::OpenFile(Int_t n)
 290 {
 291 //opens fFile with kine tree
 292
 293  const TString& dirname = GetDirName(n);
 294  if (dirname == "")
 295   {
 296    Error("OpenFiles","Can not get directory name");
 297    return 0x0;
 298   }
 299  TString filename = dirname +"/"+ fESDFileName;
 300  TFile *ret = TFile::Open(filename.Data());
 301
 302  if ( ret == 0x0)
 303   {
 304     Error("OpenFiles","Can't open fFile %s",filename.Data());
 305     return 0x0;
 306   }
 307  if (!ret->IsOpen())
 308   {
 309     Error("OpenFiles","Can't open fFile  %s",filename.Data());
 310     return 0x0;
 311   }
 312
 313  if (fReadParticles)
 314   {
 315    fRunLoader = AliRunLoader::Open(dirname +"/"+ fGAlFileName);
 316    if (fRunLoader == 0x0)
 317     {
 318       Error("OpenFiles","Can't get RunLoader for directory %s",dirname.Data());
 319       delete ret;
 320       return 0x0;
 321     }
 322    fRunLoader->LoadHeader();
 323    if (fRunLoader->LoadKinematics())
 324     {
 325       Error("Next","Error occured while loading kinematics.");
 326       delete fRunLoader;
 327       delete ret;
 328       return 0x0;
 329     }
 330   }
 331  return ret;
 332 }
 333 /**********************************************************/
 334
 335 Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
 336 {
 337   //returns pdg code from the PID index
 338   //ask jura about charge
 339   switch (spec)
 340    {
 341      case kESDElectron:
 342        return kPositron;
 343        break;
 344      case kESDMuon:
 345        return kMuonPlus;
 346        break;
 347      case kESDPion:
 348        return kPiPlus;
 349        break;
 350      case kESDKaon:
 351        return kKPlus;
 352        break;
 353      case kESDProton:
 354        return kProton;
 355        break;
 356      default:
 357        ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
 358        break;
 359    }
 360   return 0;
 361 }