JETAN/AliJetParticlesReaderESD.cxx

   1 // $Id$
   2
   3 //____________________________________________________________________
   4 //////////////////////////////////////////////////////////////////////
   5 //                                                                  //
   6 // class AliHBTReaderESD                                            //
   7 //                                                                  //
   8 // Reader for ALICE Event Summary Data (ESD).                       //
   9 //                                                                  //
  10 // Piotr.Skowronski@cern.ch                                         //
  11 //                                                                  //
  12 //////////////////////////////////////////////////////////////////////
  13
  14 #include <TMath.h>
  15 #include <TPDGCode.h>
  16 #include <TString.h>
  17 #include <TObjString.h>
  18 #include <TTree.h>
  19 #include <TFile.h>
  20 #include <TKey.h>
  21 #include <TH1.h>
  22 #include <AliESD.h>
  23 #include <AliESDtrack.h>
  24 #include <AliJetEventParticles.h>
  25 #include "AliJetParticlesReaderESD.h"
  26
  27 ClassImp(AliJetParticlesReaderESD)
  28
  29 AliJetParticlesReaderESD::AliJetParticlesReaderESD(const Char_t* esdfilename) :
  30   AliJetParticlesReader(),
  31   fESDFileName(esdfilename),
  32   fFile(0),
  33   fKeyIterator(0),
  34   fPassFlag(AliESDtrack::kTPCrefit)
  35 {
  36   //constructor
  37 }
  38
  39 /********************************************************************/
  40
  41 AliJetParticlesReaderESD::AliJetParticlesReaderESD(
  42                                       TObjArray* dirs,
  43                                       const Char_t* esdfilename) :
  44   AliJetParticlesReader(dirs),
  45   fESDFileName(esdfilename),
  46   fFile(0),
  47   fKeyIterator(0),
  48   fPassFlag(AliESDtrack::kTPCrefit)
  49 {
  50   //constructor
  51 }
  52
  53 /********************************************************************/
  54
  55 AliJetParticlesReaderESD::~AliJetParticlesReaderESD()
  56 {
  57   //desctructor
  58   if(fFile) delete fFile;
  59   if(fKeyIterator) delete fKeyIterator;
  60 }
  61
  62
  63 Int_t AliJetParticlesReaderESD::ReadNext()
  64 {
  65   //reads next event from fFile
  66
  67   do   // is OK even if 0 dirs specified,
  68     {  // in that case we try to read from "./"
  69       if (fFile == 0)
  70         {
  71           fFile = OpenFile(fCurrentDir);
  72           if (fFile == 0)
  73             {
  74               Error("ReadNext","Cannot get fFile for dir no. %d",fCurrentDir);
  75               fCurrentDir++;
  76               continue;
  77             }
  78
  79           fKeyIterator = new TIter(fFile->GetListOfKeys());
  80           fCurrentEvent = 0;
  81           //fFile->Dump();
  82           //fFile->GetListOfKeys()->Print();
  83         }
  84
  85       TKey* key = (TKey*)fKeyIterator->Next();
  86       if (key == 0)
  87         {
  88           fCurrentDir++;
  89           delete fKeyIterator;
  90           fKeyIterator = 0;
  91           delete fFile; //we have to assume there are no more ESD objects in the fFile
  92           fFile = 0;
  93           continue;
  94         }
  95
  96       TString esdname = "ESD";
  97       esdname+=fCurrentEvent;
  98       AliESD* esd = dynamic_cast<AliESD*>(fFile->Get(esdname));
  99       if (esd == 0)
 100       {
 101         //Info("ReadNext","This key is not an AliESD object %s",key->GetName());
 102         Info("ReadNext","Can not find AliESD object named %s",esdname.Data());
 103
 104         fCurrentDir++;
 105         delete fKeyIterator;
 106         fKeyIterator = 0;
 107         delete fFile;//we have to assume there is no more ESD objects in the fFile
 108         fFile = 0;
 109         continue;
 110       }
 111
 112       ReadESD(esd);
 113
 114       fCurrentEvent++;
 115       fNEventsRead++;
 116       delete esd;
 117       return kTRUE;//success -> read one event
 118     }  while(fCurrentDir < GetNumberOfDirs());
 119        //end of loop over directories specified in fDirs Obj Array
 120   return kFALSE; //no more directories to read
 121 }
 122
 123 /**********************************************************/
 124
 125 Int_t AliJetParticlesReaderESD::ReadESD(AliESD* esd)
 126 {
 127   //Reads one ESD
 128
 129   if (esd == 0)
 130    {
 131      Error("ReadESD","ESD is NULL");
 132      return kFALSE;
 133    }
 134
 135   //TDatabasePDG* pdgdb = TDatabasePDG::Instance();
 136   //if (pdgdb == 0)
 137   //{
 138   //   Error("ReadESD","Can not get PDG Database Instance.");
 139   //   return kFALSE;
 140   // }
 141
 142   //Float_t mf = esd->GetMagneticField();
 143   //if ( (mf == 0.0) && (fNTrackPoints > 0) )
 144   //{
 145   //   Error("ReadESD","Magnetic Field is 0 and Track Points demanded. Skipping to next event.");
 146   //   return kFALSE;
 147   //}
 148
 149   Info("ReadESD","Reading Event %d",fCurrentEvent);
 150   if((!fOwner) || (fEventParticles==0))
 151     fEventParticles = new AliJetEventParticles();
 152
 153   Double_t vertexpos[3];//vertex position
 154   const AliESDVertex* kvertex = esd->GetVertex();
 155   if (kvertex == 0)
 156    {
 157      Info("ReadESD","ESD returned NULL pointer to vertex - assuming (0.0,0.0,0.0)");
 158      vertexpos[0] = 0.0;
 159      vertexpos[1] = 0.0;
 160      vertexpos[2] = 0.0;
 161    }
 162   else
 163    {
 164      kvertex->GetXYZ(vertexpos);
 165    }
 166   fEventParticles->SetVertex(vertexpos[0],vertexpos[1],vertexpos[2]);
 167
 168   const Int_t kntr = esd->GetNumberOfTracks();
 169   Info("ReadESD","Found %d tracks.",kntr);
 170   for (Int_t i = 0;i<kntr; i++)
 171    {
 172      const AliESDtrack *kesdtrack = esd->GetTrack(i);
 173      if (kesdtrack == 0)
 174       {
 175         Error("ReadESD","Can not get track %d", i);
 176         continue;
 177       }
 178
 179      if ((kesdtrack->GetStatus() & fPassFlag) == kFALSE)
 180       {
 181         Info("ReadNext","Particle skipped.");
 182         continue;
 183       }
 184
 185      Double_t mom[3];  //momentum
 186      kesdtrack->GetPxPyPz(mom);
 187      //kesdtrack->GetConstrainedPxPyPz(mom);
 188      //Double_t pos[3];//position
 189      //kesdtrack->GetXYZ(pos);
 190      //kesdtrack->GetConstrainedXYZ(pos);
 191      const Float_t kmass=kesdtrack->GetMass();
 192      const Float_t kp2=mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2];
 193      const Float_t ketot=TMath::Sqrt(kmass*kmass+kp2);
 194      const Float_t kpt=TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
 195      const Float_t kp=TMath::Sqrt(kp2);
 196      const Float_t keta=0.5*TMath::Log((kp+mom[2]+1e-30)/(kp-mom[2]+1e-30));
 197      const Float_t kphi=TMath::Pi()+TMath::ATan2(-mom[1],-mom[0]);
 198      if(IsAcceptedParticle(kpt,kphi,keta))
 199        fEventParticles->AddParticle(mom[0],mom[1],mom[2],ketot,i,kesdtrack->GetLabel(),kpt,kphi,keta);
 200
 201    } // loop over tracks
 202
 203   return kTRUE;
 204 }
 205
 206 /**********************************************************/
 207
 208 void AliJetParticlesReaderESD::Rewind()
 209 {
 210   //rewinds reading
 211
 212   if(fFile) delete fFile;
 213   if(fKeyIterator) delete fKeyIterator;
 214   fKeyIterator = 0;
 215   fFile = 0;
 216   fCurrentDir = 0;
 217   fNEventsRead = 0;
 218 }
 219
 220 /**********************************************************/
 221
 222 TFile* AliJetParticlesReaderESD::OpenFile(Int_t n)
 223 {
 224   //opens fFile with kine tree
 225
 226  const TString& dirname = GetDirName(n);
 227  if (dirname == "")
 228   {
 229    Error("OpenFiles","Can not get directory name");
 230    return 0;
 231   }
 232  TString filename = dirname +"/"+ fESDFileName;
 233  TFile *ret = TFile::Open(filename.Data());
 234  if (ret == 0)
 235   {
 236     Error("OpenFiles","Can't open fFile %s",filename.Data());
 237     return 0;
 238   }
 239  if (!ret->IsOpen())
 240   {
 241     Error("OpenFiles","Can't open fFile  %s",filename.Data());
 242     return 0;
 243   }
 244
 245  return ret;
 246 }