HBTAN/AliHBTReaderKineTree.cxx

   1 #include "AliHBTReaderKineTree.h"
   2
   3 #include <Riostream.h>
   4 //#include <Riostream.h>
   5 #include <TString.h>
   6 #include <TObjString.h>
   7 #include <TTree.h>
   8 #include <TFile.h>
   9 #include <TParticle.h>
  10
  11 #include <AliRun.h>
  12 #include <AliStack.h>
  13 #include <AliHeader.h>
  14
  15 #include "AliHBTRun.h"
  16 #include "AliHBTEvent.h"
  17 #include "AliHBTParticle.h"
  18 #include "AliHBTParticleCut.h"
  19
  20 ClassImp(AliHBTReaderKineTree)
  21 /**********************************************************/
  22
  23 AliHBTReaderKineTree::
  24 AliHBTReaderKineTree():fFileName("galice.root")
  25 {
  26   fParticles = new AliHBTRun();
  27   fIsRead = kFALSE;
  28 }
  29
  30 AliHBTReaderKineTree::
  31 AliHBTReaderKineTree(TString& fname):fFileName(fname)
  32 {
  33   fParticles = new AliHBTRun();
  34   fIsRead = kFALSE;
  35 }
  36
  37
  38 /**********************************************************/
  39 AliHBTReaderKineTree::
  40 AliHBTReaderKineTree(TObjArray* dirs,const Char_t *filename):AliHBTReader(dirs),fFileName(filename)
  41 {
  42   fParticles = new AliHBTRun();
  43   fIsRead = kFALSE;
  44 }
  45 /**********************************************************/
  46
  47 AliHBTEvent* AliHBTReaderKineTree::GetParticleEvent(Int_t n)
  48  {
  49  //returns Nth event with simulated particles
  50    if (!fIsRead)
  51     if(Read(fParticles,0x0))
  52      {
  53        Error("GetParticleEvent","Error in reading");
  54        return 0x0;
  55      }
  56
  57    return fParticles->GetEvent(n);
  58  }
  59
  60 /********************************************************************/
  61
  62 Int_t AliHBTReaderKineTree::GetNumberOfPartEvents()
  63  {
  64  //returns number of events of particles
  65    if (!fIsRead)
  66     if(Read(fParticles,0x0))
  67      {
  68        Error("GetNumberOfPartEvents","Error in reading");
  69        return 0;
  70      }
  71    return fParticles->GetNumberOfEvents();
  72  }
  73
  74
  75 /**********************************************************/
  76 Int_t AliHBTReaderKineTree::
  77 Read(AliHBTRun* particles, AliHBTRun *tracks)
  78  {
  79  cout<<"AliHBTReaderKineTree::Read()"<<endl;
  80  if (!particles) //check if an object is instatiated
  81    {
  82      Error("Read"," particles object must instatiated before passing it to the reader");
  83      return 1;
  84    }
  85  particles->Reset();//clear runs == delete all old events
  86
  87  Int_t Ndirs;//total number of directories specified in fDirs array
  88  Int_t Nevents; //number of events read in current directory
  89  Int_t totalNevents = 0; //total number of read events
  90  Int_t currentdir = 0; //number of current directory name is fDirs array
  91
  92  if (fDirs) //if array with directories is supplied by user
  93   {
  94     Ndirs = fDirs->GetEntries(); //get the number if directories
  95   }
  96  else
  97   {
  98     Ndirs = 0; //if the array is not supplied read only from current directory
  99   }
 100
 101  do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 102   {
 103     cout<<"________________________________________________________\n";
 104     TFile* fileK = OpenFile(currentdir);
 105     TTree* treeE = (TTree*)fileK->Get("TE");
 106     if(treeE)
 107      {
 108       Nevents = (Int_t)treeE->GetEntries();
 109       cout<<"Found "<<Nevents<<" in directory "<<GetDirName(currentdir)<<endl;
 110       fileK->Delete("TE");
 111      }
 112     else
 113      {
 114       Error("Read","Cannot find Header Tree (TE)");
 115       Nevents = 0;
 116      }
 117
 118     for(Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)//loop over all events
 119      {
 120 //      cout<<"processing event "<<currentEvent<<" in current dir."<<endl;
 121       AliStack* stack = GetStack(currentEvent,fileK);
 122       if (!stack)
 123        {
 124          Error("Read","Can not get stack for event %d",currentEvent);
 125          continue;
 126        }
 127       Int_t npart = stack->GetNtrack();
 128       Int_t nnn = 0;
 129       for (Int_t i = 0;i<npart; i++)
 130        {
 131
 132          TParticle * p = stack->Particle(i);
 133          if (p->GetFirstMother() >= 0) continue;
 134
 135          if(Pass(p->GetPdgCode())) continue; //check if we are intersted with particles of this type
 136                                               //if not take next partilce
 137
 138          AliHBTParticle* part = new AliHBTParticle(*p);
 139          if(Pass(part)) { delete part; continue;}//check if meets all criteria of any of our cuts
 140                                                   //if it does not delete it and take next good track
 141          particles->AddParticle(totalNevents,part);//put track and particle on the run
 142
 143          if ( (nnn%100) == 0)
 144           {
 145             cout<<nnn<<"\r";
 146             fflush(0);
 147           }
 148          nnn++;
 149        }
 150       cout<<"Total read "<<nnn<<endl;
 151       delete stack;
 152       totalNevents++;
 153      }
 154     delete gAlice;
 155     gAlice = 0;
 156     if (fileK)
 157      {
 158       fileK->Close();
 159       delete fileK;
 160       fileK = 0;
 161      }
 162
 163     currentdir++;
 164   }while(currentdir < Ndirs);//end of loop over directories specified in fDirs Obj Array
 165  fIsRead = kTRUE;
 166  return 0;
 167 }
 168
 169 /**********************************************************/
 170 AliStack* AliHBTReaderKineTree::GetStack(Int_t n, TFile* file)
 171 {
 172   AliHeader *header =  new AliHeader();
 173   TTree* treeE = (TTree*)file->Get("TE");
 174   if(treeE)
 175    {
 176     treeE->SetBranchAddress("Header", &header);
 177     if (!treeE->GetEntry(n))
 178      {
 179        Error("GetEvent","Cannot find event:%dn",n);
 180        delete header;
 181        treeE->Delete("TE");
 182        return 0x0;
 183      }
 184    }
 185   else
 186    {
 187      Error("GetStack","Cannot find Header Tree (TE)n");
 188      delete header;
 189      return 0x0;
 190    }
 191
 192   AliStack* stack = header->Stack();
 193   if (stack)
 194    {
 195     if (!stack->GetEvent(n))
 196      {
 197        delete header;
 198        treeE->Delete("TE");
 199        return 0x0;
 200      }
 201    }
 202   delete header;
 203   treeE->Delete("TE");
 204   return stack;
 205 }
 206
 207 TFile* AliHBTReaderKineTree::OpenFile(Int_t n)
 208 {
 209 //opens file with kine tree
 210
 211  const TString& dirname = GetDirName(n);
 212  if (dirname == "")
 213   {
 214    Error("OpenFiles","Can not get directory name");
 215    return 0x0;
 216   }
 217  TString filename = dirname +"/"+ fFileName;
 218  TFile *ret = TFile::Open(filename.Data());
 219
 220  if ( ret == 0x0)
 221   {
 222     Error("OpenFiles","Can't open file %s",filename.Data());
 223     return 0x0;
 224   }
 225  if (!ret->IsOpen())
 226   {
 227     Error("OpenFiles","Can't open file  %s",filename.Data());
 228     return 0x0;
 229   }
 230
 231  return ret;
 232 }