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