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