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