Reader for ESD
authorskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Aug 2003 15:34:55 +0000 (15:34 +0000)
committerskowron <skowron@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 4 Aug 2003 15:34:55 +0000 (15:34 +0000)
HBTAN/AliHBTReaderESD.cxx [new file with mode: 0644]
HBTAN/AliHBTReaderESD.h [new file with mode: 0644]
HBTAN/HBTAnalysisLinkDef.h
HBTAN/libHBTAN.pkg

diff --git a/HBTAN/AliHBTReaderESD.cxx b/HBTAN/AliHBTReaderESD.cxx
new file mode 100644 (file)
index 0000000..0dbffd2
--- /dev/null
@@ -0,0 +1,332 @@
+#include "AliHBTReaderESD.h"
+
+#include <TPDGCode.h>
+#include <TString.h>
+#include <TObjString.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TParticle.h>
+
+#include <AliESDtrack.h>
+#include <AliESD.h>
+
+#include "AliHBTRun.h"
+#include "AliHBTEvent.h"
+#include "AliHBTParticle.h"
+#include "AliHBTParticleCut.h"
+
+ClassImp(AliHBTReaderESD)
+
+AliHBTReaderESD::AliHBTReaderESD(const Char_t* esdfilename):
+ fParticles(new AliHBTRun()),
+ fTracks(new AliHBTRun()),
+ fESDFileName(esdfilename),
+ fIsRead(kFALSE)
+{
+  //cosntructor
+  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
+    Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
+}
+/********************************************************************/
+  
+AliHBTReaderESD::AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename):
+ AliHBTReader(dirs), 
+ fParticles(new AliHBTRun()),
+ fTracks(new AliHBTRun()),
+ fESDFileName(esdfilename),
+ fIsRead(kFALSE)
+{
+  //cosntructor
+  if ( ((Int_t)kNSpecies) != ((Int_t)AliESDtrack::kSPECIES))
+    Fatal("AliHBTReaderESD","ESD defintions probobly changed. Ask Youra.");
+}
+/********************************************************************/
+
+AliHBTReaderESD::~AliHBTReaderESD()
+{
+ //desctructor
+  delete fParticles;
+  delete fTracks;
+}
+/********************************************************************/
+
+AliHBTEvent* AliHBTReaderESD::GetParticleEvent(Int_t n)
+ {
+ //returns Nth event with simulated particles
+   if (!fIsRead) 
+    if(Read(fParticles,fTracks))
+     {
+       Error("GetParticleEvent","Error in reading");
+       return 0x0;
+     }
+   return fParticles->GetEvent(n);
+ }
+/********************************************************************/
+
+AliHBTEvent* AliHBTReaderESD::GetTrackEvent(Int_t n)
+ {
+ //returns Nth event with reconstructed tracks
+   if (!fIsRead) 
+    if(Read(fParticles,fTracks))
+     {
+       Error("GetTrackEvent","Error in reading");
+       return 0x0;
+     }
+   return fTracks->GetEvent(n);
+ }
+/********************************************************************/
+
+Int_t AliHBTReaderESD::GetNumberOfPartEvents()
+ {
+ //returns number of events of particles
+   if (!fIsRead) 
+    if ( Read(fParticles,fTracks))
+     {
+       Error("GetNumberOfPartEvents","Error in reading");
+       return 0;
+     }
+   return fParticles->GetNumberOfEvents();
+ }
+
+/********************************************************************/
+Int_t AliHBTReaderESD::GetNumberOfTrackEvents()
+ {
+ //returns number of events of tracks
+  if (!fIsRead)
+    if(Read(fParticles,fTracks))
+     {
+       Error("GetNumberOfTrackEvents","Error in reading");
+       return 0;
+     }
+  return fTracks->GetNumberOfEvents();
+ }
+/********************************************************************/
+
+
+Int_t AliHBTReaderESD::Read(AliHBTRun* particles, AliHBTRun *tracks)
+{
+ //reads data and puts put to the particles and tracks objects
+ //reurns 0 if everything is OK
+ //
+  Info("Read","");
+  Int_t totalNevents = 0;
+  Int_t currentdir = 0; //number of current directory name is fDirs array
+  
+  Double_t pidtable[kNSpecies];//array used for reading pid probabilities from ESD track
+  Double_t w[kNSpecies];
+  Double_t mom[3];//momentum
+  Double_t pos[3];//position
+   //****** Tentative particle type "concentrations"
+  const Double_t concentr[5]={0.05, 0., 0.85, 0.10, 0.05};
+  
+  TDatabasePDG* pdgdb = TDatabasePDG::Instance();
+  if (pdgdb == 0x0)
+   {
+     Error("Read","Can not get PDG Database Instance.");
+     return 1;
+   }
+  if (!particles) //check if an object is instatiated
+   {
+     Error("Read"," particles object must instatiated before passing it to the reader");
+     return 1;
+   }
+  if (!tracks)  //check if an object is instatiated
+   {
+     Error("Read"," tracks object must instatiated before passing it to the reader");
+     return 1;
+   }
+  particles->Reset();//clear runs == delete all old events
+  tracks->Reset();
+
+  Int_t Ndirs;
+  if (fDirs) //if array with directories is supplied by user
+   {
+     Ndirs = fDirs->GetEntries(); //get the number if directories
+   }
+  else
+   {
+     Ndirs = 0; //if the array is not supplied read only from current directory
+   }
+
+  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
+   {
+     TFile* file = OpenFile(currentdir);
+     if (file == 0x0)
+       {
+         Error("Read","Cannot get File for dir no. %d",currentdir);
+         currentdir++;
+         continue;
+       } 
+       
+     for(Int_t currentEvent = 0; ;currentEvent++)//loop over all events
+      {
+       TString esdname;
+       esdname+=currentEvent;
+       AliESD* esd = dynamic_cast<AliESD*>(file->Get(esdname));
+       if (esd == 0x0)
+         {
+           if (AliHBTParticle::fgDebug > 2 )
+            {
+              Info("Read","Can not find ESD object named %s.",esdname.Data());
+            }
+           currentdir++;
+           break;//we have to assume there is no more ESD objects in the file
+         } 
+       
+       Info("Read","Reading Event %d",currentEvent);
+  
+       Int_t ntr = esd->GetNumberOfTracks();
+       Info("Read","Found %d tracks.",ntr);
+       for (Int_t i = 0;i<ntr; i++)
+        {
+          AliESDtrack *esdtrack = esd->GetTrack(i);
+          if (esdtrack == 0x0)
+           {
+             Error("Read","Can not get track %d", i);
+             continue;
+           }
+          if ((esdtrack->GetStatus()&AliESDtrack::kESDpid) == kFALSE) 
+           {
+             if (AliHBTParticle::fgDebug > 2) 
+               Info("Read","Particle skipped: PID BIT is not set.");
+             continue;
+           }
+
+          esdtrack->GetESDpid(pidtable);
+          esdtrack->GetPxPyPz(mom);
+          esdtrack->GetXYZ(pos);
+          Double_t rc=0.;
+
+           //Here we apply Bayes' formula
+          for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) rc+=concentr[s]*pidtable[s];
+          if (rc==0.0) 
+           {
+             if (AliHBTParticle::fgDebug > 2) 
+               Info("Read","Particle rejected since total bayessian PID probab. is zero.");
+             continue;
+           }
+
+          for (Int_t s=0; s<AliESDtrack::kSPECIES; s++) w[s]=concentr[s]*pidtable[s]/rc;
+
+          if (AliHBTParticle::fgDebug > 4)
+           { 
+             Info("Read","###########################################################################");
+             Info("Read","Momentum: %f %f %f",mom[0],mom[1],mom[2]);
+             Info("Read","Position: %f %f %f",pos[0],pos[1],pos[2]);
+             Info("Read","Radius: %f ,rc: %f",TMath::Hypot(pos[0],pos[1]),rc);
+             TString msg("Pid list got from track:");
+             for (Int_t s = 0;s<kNSpecies;s++)
+              {
+                msg+="\n    ";
+                msg+=s;
+                msg+="(";
+                msg+=GetSpeciesPdgCode((ESpecies)s);
+                msg+="): ";
+                msg+=w[s];
+                msg+=" (";
+                msg+=pidtable[s];
+                msg+=")";
+              }
+             Info("Read","%s",msg.Data());
+           }//if (AliHBTParticle::fgDebug>4)
+           
+          for (Int_t s = 0; s<kNSpecies; s++)
+           {
+             if (w[s] == 0.0) continue;
+
+             Int_t pdgcode = GetSpeciesPdgCode((ESpecies)s);
+             if(Pass(pdgcode)) continue; //check if we are intersted with particles of this type 
+
+             Double_t mass = pdgdb->GetParticle(pdgcode)->Mass();
+             Double_t tEtot = TMath::Sqrt( mom[0]*mom[0] + mom[1]*mom[1] + mom[2]*mom[2] + mass*mass);//total energy of the track
+
+             AliHBTParticle* track = new AliHBTParticle(pdgcode, w[s],i, 
+                                                        mom[0], mom[1], mom[2], tEtot,
+                                                        pos[0], pos[1], pos[2], 0.);
+             //copy probabilitis of other species (if not zero)
+             for (Int_t k = 0; k<kNSpecies; k++)
+              {
+                if (k == s) continue;
+                if (w[k] == 0.0) continue;
+                track->SetPIDprobability(GetSpeciesPdgCode( (ESpecies)k ),w[k]);
+              }
+             if(Pass(track))//check if meets all criteria of any of our cuts
+                            //if it does not delete it and take next good track
+               { 
+                 delete track;
+                 continue;
+               }
+             tracks->AddParticle(totalNevents,track);
+             if (AliHBTParticle::fgDebug > 4 )
+              {
+                Info("Read","\n\nAdding Particle with incarnation %d",pdgcode);
+                track->Print();
+              }
+           }//for (Int_t s = 0; s<kNSpecies; s++)
+
+        }//for (Int_t i = 0;i<ntr; i++)  -- loop over tracks
+       totalNevents++;
+      }//for(Int_t currentEvent = 0; ;currentEvent++) -- for over ESD in single file
+      
+   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
+  
+  fIsRead = kTRUE;
+  return 0;
+}      
+/**********************************************************/
+   
+TFile* AliHBTReaderESD::OpenFile(Int_t n)
+{
+//opens file with kine tree
+
+ const TString& dirname = GetDirName(n);
+ if (dirname == "")
+  {
+   Error("OpenFiles","Can not get directory name");
+   return 0x0;
+  }
+ TString filename = dirname +"/"+ fESDFileName;
+ TFile *ret = TFile::Open(filename.Data()); 
+
+ if ( ret == 0x0)
+  {
+    Error("OpenFiles","Can't open file %s",filename.Data());
+    return 0x0;
+  }
+ if (!ret->IsOpen())
+  {
+    Error("OpenFiles","Can't open file  %s",filename.Data());
+    return 0x0;
+  }
+ return ret;
+}
+/**********************************************************/
+
+Int_t AliHBTReaderESD::GetSpeciesPdgCode(ESpecies spec)//skowron
+{
+  //returns pdg code from the PID index
+  //ask jura about charge
+  switch (spec)
+   {
+     case kESDElectron:
+       return kPositron;
+       break;
+     case kESDMuon:
+       return kMuonPlus;
+       break;
+     case kESDPion:
+       return kPiPlus;
+       break;
+     case kESDKaon:
+       return kKPlus;
+       break;
+     case kESDProton:
+       return kProton;
+       break;
+     default:
+       ::Warning("GetSpeciesPdgCode","Specie with number %d is not defined.",(Int_t)spec);
+       break;
+   }
+  return 0;
+}
diff --git a/HBTAN/AliHBTReaderESD.h b/HBTAN/AliHBTReaderESD.h
new file mode 100644 (file)
index 0000000..ca1c86a
--- /dev/null
@@ -0,0 +1,54 @@
+#ifndef ALIHBTREADERESD_H
+#define ALIHBTREADERESD_H
+
+#include "AliHBTReader.h"
+//___________________________________________________________________________
+/////////////////////////////////////////////////////////////////////////////
+//                                                                         //
+// Multi file reader for ESD                                               //
+//                                                                         //
+// This reader reads tracks from Event Summary Data                        //
+// do not read particles                                                   //
+// Piotr.Skowronski@cern.ch                                                //
+// more info: http://alisoft.cern.ch/people/skowron/analyzer/index.html    //
+//                                                                         //
+/////////////////////////////////////////////////////////////////////////////
+
+#include <TString.h>
+class TFile;
+
+class AliHBTReaderESD: public AliHBTReader
+{
+  public:
+    AliHBTReaderESD(const Char_t* esdfilename = "AliESDs.root");
+
+    AliHBTReaderESD(TObjArray* dirs,const Char_t* esdfilename = "AliESDs.root");
+
+    virtual ~AliHBTReaderESD();
+    
+    Int_t Read(AliHBTRun* particles, AliHBTRun *tracks);//reads tracks and particles and puts them in runs
+    
+    AliHBTEvent* GetParticleEvent(Int_t);//returns pointer to event with particles
+    AliHBTEvent* GetTrackEvent(Int_t);//returns pointer to event with particles 
+    Int_t        GetNumberOfPartEvents();//returns number of particle events
+    Int_t        GetNumberOfTrackEvents();//returns number of track events
+    
+    enum ESpecies {kESDElectron = 0, kESDMuon, kESDPion, kESDKaon, kESDProton, kNSpecies};
+    static Int_t GetSpeciesPdgCode(ESpecies spec);//skowron
+  protected:
+      
+    TFile*       OpenFile(Int_t evno);//opens files to be read for given event
+    void         CloseFiles(TFile*);//close files
+    
+    AliHBTRun*   fParticles; //!simulated particles
+    AliHBTRun*   fTracks; //!reconstructed tracks (particles)
+
+    TString      fESDFileName;//name of the file with tracks
+    Bool_t       fIsRead;//!flag indicating if the data are already read
+    
+  private:
+    ClassDef(AliHBTReaderESD,1)
+};
+
+
+#endif
index 16388e2c3bc724f0b865b90e7ec3bf6640b8e4dd..c4b51afd2f3d00c5089b3e357f22bc40fd07e077 100644 (file)
@@ -67,6 +67,7 @@
 #pragma link C++ class AliHBTOrCut+;
 
 #pragma link C++ class AliHBTReader+;
+#pragma link C++ class AliHBTReaderESD+;
 #pragma link C++ class AliHBTReaderTPC+;
 #pragma link C++ class AliHBTReaderITSv1+;
 #pragma link C++ class AliHBTReaderITSv2+;
index 5706193471cb84e595b281b02a24bff1ddddf1c1..f64184ed6be1ef97f1a64ca2b3e2e1c33e068f40 100644 (file)
@@ -1,20 +1,22 @@
-SRCS          =   \
-AliHBTAnalysis.cxx    AliHBTPair.cxx           AliHBTReaderPPprod.cxx \
-AliHBTReader.cxx      AliHBTReaderTPC.cxx \
-AliHBTCorrelFctn.cxx  AliHBTPairCut.cxx \
-AliHBTEvent.cxx       AliHBTParticle.cxx       AliHBTRun.cxx \
-AliHBTFunction.cxx    AliHBTMonitorFunction.cxx   AliHBTParticleCut.cxx  \
-AliHBTReaderITSv1.cxx AliHBTReaderITSv2.cxx    AliHBTReaderKineTree.cxx \
-AliHBTTwoTrackEffFctn.cxx       AliHBTReaderInternal.cxx\
-AliHBTQResolutionFctns.cxx      AliHBTQDistributionFctns.cxx    \
-AliHBTLLWeights.cxx             AliHBTLLWeightFctn.cxx\
-AliHBTMonDistributionFctns.cxx     AliHBTMonResolutionFctns.cxx \
-AliHBTLLWeightsPID.cxx AliHBTLLWeightTheorFctn.cxx \
+SRCS          =       AliHBTAnalysis.cxx  \
+AliHBTReader.cxx      AliHBTReaderKineTree.cxx \
+AliHBTReaderESD.cxx   AliHBTReaderInternal.cxx \
+AliHBTReaderTPC.cxx   AliHBTReaderPPprod.cxx \
+AliHBTReaderITSv1.cxx AliHBTReaderITSv2.cxx    \
+AliHBTParticle.cxx    AliHBTParticleCut.cxx \
+AliHBTPair.cxx        AliHBTPairCut.cxx \
+AliHBTEvent.cxx       AliHBTRun.cxx \
+AliHBTFunction.cxx    AliHBTCorrelFctn.cxx  \
+AliHBTMonitorFunction.cxx       AliHBTTwoTrackEffFctn.cxx \
+AliHBTQResolutionFctns.cxx      AliHBTQDistributionFctns.cxx \
+AliHBTMonDistributionFctns.cxx  AliHBTMonResolutionFctns.cxx \
+AliHBTLLWeights.cxx             AliHBTLLWeightFctn.cxx \
+AliHBTLLWeightsPID.cxx          AliHBTLLWeightTheorFctn.cxx \
 AliHBTPositionRandomizer.cxx
 
 FSRCS   = fsiini.F  fsiw.F  led_bldata.F  ltran12.F
 
-HDRS:= $(SRCS:.cxx=.h) 
+HDRS:= $(SRCS:.cxx=.h)
 
 DHDR= HBTAnalysisLinkDef.h