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