HBTAN/AliHBTReaderInternal.cxx

   1 #include "AliHBTReaderInternal.h"
   2
   3 #include <iostream.h>
   4 //#include <fstream.h>
   5 #include <TTree.h>
   6 #include <TFile.h>
   7 #include <TParticle.h>
   8
   9 #include <AliRun.h>
  10 #include <AliMagF.h>
  11
  12 #include "AliHBTRun.h"
  13 #include "AliHBTEvent.h"
  14 #include "AliHBTParticle.h"
  15 #include "AliHBTParticleCut.h"
  16
  17
  18 AliHBTReaderInternal nnn;
  19
  20 ClassImp(AliHBTReaderInternal)
  21 /********************************************************************/
  22
  23 AliHBTReaderInternal::AliHBTReaderInternal()
  24 {
  25   fParticles = 0x0;
  26   fTracks = 0x0;
  27   fIsRead = kFALSE;
  28 }
  29 /********************************************************************/
  30
  31 AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):fFileName(filename)
  32 {
  33  fParticles = new AliHBTRun();
  34  fTracks    = new AliHBTRun();
  35  fIsRead = kFALSE;
  36 }
  37 /********************************************************************/
  38 AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
  39   AliHBTReader(dirs),fFileName(filename)
  40 {
  41  fParticles = new AliHBTRun();
  42  fTracks    = new AliHBTRun();
  43  fIsRead = kFALSE;
  44 }
  45 AliHBTReaderInternal::~AliHBTReaderInternal()
  46  {
  47  //desctructor
  48    delete fParticles;
  49    delete fTracks;
  50  }
  51 /********************************************************************/
  52 AliHBTEvent* AliHBTReaderInternal::GetParticleEvent(Int_t n)
  53  {
  54  //returns Nth event with simulated particles
  55    if (!fIsRead)
  56     if(Read(fParticles,fTracks))
  57      {
  58        Error("GetParticleEvent","Error in reading");
  59        return 0x0;
  60      }
  61    return fParticles->GetEvent(n);
  62  }
  63 /********************************************************************/
  64 AliHBTEvent* AliHBTReaderInternal::GetTrackEvent(Int_t n)
  65  {
  66  //returns Nth event with reconstructed tracks
  67    if (!fIsRead)
  68     if(Read(fParticles,fTracks))
  69      {
  70        Error("GetTrackEvent","Error in reading");
  71        return 0x0;
  72      }
  73    return fTracks->GetEvent(n);
  74  }
  75 /********************************************************************/
  76
  77 Int_t AliHBTReaderInternal::GetNumberOfPartEvents()
  78  {
  79  //returns number of events of particles
  80    if (!fIsRead)
  81     if ( Read(fParticles,fTracks))
  82      {
  83        Error("GetNumberOfPartEvents","Error in reading");
  84        return 0;
  85      }
  86    return fParticles->GetNumberOfEvents();
  87  }
  88
  89 /********************************************************************/
  90 Int_t AliHBTReaderInternal::GetNumberOfTrackEvents()
  91  {
  92  //returns number of events of tracks
  93   if (!fIsRead)
  94     if(Read(fParticles,fTracks))
  95      {
  96        Error("GetNumberOfTrackEvents","Error in reading");
  97        return 0;
  98      }
  99   return fTracks->GetNumberOfEvents();
 100  }
 101 /********************************************************************/
 102
 103
 104 Int_t AliHBTReaderInternal::Read(AliHBTRun* particles, AliHBTRun *tracks)
 105  {
 106  //reads data and puts put to the particles and tracks objects
 107  //reurns 0 if everything is OK
 108  //
 109   cout<<"AliHBTReaderInternal::Read()"<<endl;
 110   Int_t i; //iterator and some temprary values
 111   Int_t Nevents = 0;
 112   TFile *aFile;//file with tracks
 113   AliHBTParticle* p = 0x0;
 114
 115   if (!particles) //check if an object is instatiated
 116    {
 117      Error("Read"," particles object must instatiated before passing it to the reader");
 118    }
 119   if (!tracks)  //check if an object is instatiated
 120    {
 121      Error("Read"," tracks object must instatiated before passing it to the reader");
 122    }
 123   particles->Reset();//clear runs == delete all old events
 124   tracks->Reset();
 125
 126   Int_t currentdir = 0;
 127
 128   Int_t Ndirs;
 129   if (fDirs) //if array with directories is supplied by user
 130    {
 131      Ndirs = fDirs->GetEntries(); //get the number if directories
 132    }
 133   else
 134    {
 135      Ndirs = 0; //if the array is not supplied read only from current directory
 136    }
 137
 138   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
 139   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
 140
 141   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 142    {
 143
 144     if( (i=OpenFile(aFile, currentdir)) )
 145      {
 146        Error("Read","Exiting due to problems with opening files. Errorcode %d",i);
 147        return i;
 148      }
 149    /***************************/
 150    /***************************/
 151    /***************************/
 152
 153      TTree* tree = (TTree*)aFile->Get("data");
 154      if (tree == 0x0)
 155       {
 156        Error("Read","Can not get the tree");
 157        return 1;
 158       }
 159
 160      TBranch *trackbranch=tree->GetBranch("tracks");//get the branch with tracks
 161      if (trackbranch == 0x0) ////check if we got the branch
 162        {//if not return with error
 163          Warning("Read","Can't find a branch with tracks !\n");
 164        }
 165      else
 166       {
 167         trackbranch->SetAddress(&tbuffer);
 168       }
 169
 170      TBranch *partbranch=tree->GetBranch("particles");//get the branch with particles
 171      if (partbranch == 0x0) ////check if we got the branch
 172        {//if not return with error
 173          Warning("Read","Can't find a branch with particles !\n");
 174        }
 175      else
 176       {
 177         partbranch->SetAddress(&pbuffer);
 178       }
 179
 180      Nevents = (Int_t)tree->GetEntries();
 181      cout<<"________________________________________________________\n";
 182      cout<<"Found "<<Nevents<<" event(s) in directory "<<GetDirName(currentdir)<<endl;
 183
 184      for (Int_t currentEvent =0; currentEvent<Nevents;currentEvent++)
 185        {
 186          tree->GetEvent(currentEvent);
 187          if (partbranch)
 188           {
 189             for(i = 0; i < pbuffer->GetEntries(); i++)
 190              {
 191                p = dynamic_cast<AliHBTParticle*>(pbuffer->At(i));
 192                if(p == 0x0) continue; //if returned pointer is NULL
 193                if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 194                if(Pass(p)) continue; //check if we are intersted with particles of this type
 195                                                    //if not take next partilce
 196                AliHBTParticle* part = new AliHBTParticle(*p);
 197                particles->AddParticle(currentEvent,part);//put track and particle on the run
 198              }
 199             cout<<"Read: "<<particles->GetNumberOfParticlesInEvent(currentEvent)<<" particles  ";
 200           }
 201          else cout<<"Read: 0 particles  ";
 202
 203          if (trackbranch)
 204           {
 205             for(i = 0; i < tbuffer->GetEntries(); i++)
 206              {
 207                p = dynamic_cast<AliHBTParticle*>(tbuffer->At(i));
 208                if(p == 0x0) continue; //if returned pointer is NULL
 209                if(p->GetPDG() == 0x0) continue; //if particle has crezy PDG code (not known to our database)
 210                if(Pass(p)) continue; //check if we are intersted with particles of this type
 211                                                    //if not take next partilce
 212                AliHBTParticle* part = new AliHBTParticle(*p);
 213                tracks->AddParticle(currentEvent,part);//put track and particle on the run
 214              }
 215             cout<<tracks->GetNumberOfParticlesInEvent(currentEvent)<<" tracks"<<endl;
 216           }
 217          else cout<<" 0 tracks"<<endl;
 218
 219
 220
 221        }
 222
 223
 224    /***************************/
 225    /***************************/
 226    /***************************/
 227    currentdir++;
 228    aFile->Close();
 229    aFile = 0x0;
 230    }while(currentdir < Ndirs);
 231
 232   fIsRead = kTRUE;
 233   return 0;
 234 }
 235
 236 /********************************************************************/
 237
 238 Int_t AliHBTReaderInternal::OpenFile(TFile*& aFile,Int_t event)
 239 {
 240
 241    const TString& dirname = GetDirName(event);
 242    if (dirname == "")
 243     {
 244       Error("OpenFile","Can not get directory name");
 245       return 4;
 246     }
 247
 248    TString filename = dirname +"/"+ fFileName;
 249    aFile = TFile::Open(filename.Data());
 250    if ( aFile  == 0x0 )
 251      {
 252        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
 253        return 1;
 254      }
 255    if (!aFile->IsOpen())
 256      {
 257        Error("OpenFiles","Can't open file with tacks named %s",filename.Data());
 258        return 1;
 259      }
 260    return 0;
 261 }
 262
 263
 264
 265
 266 Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile)
 267  {
 268   //reads tracks from runs and writes them to file
 269     Int_t i,j;
 270
 271
 272   TFile *histoOutput = TFile::Open(outfile,"recreate");
 273
 274   if (!histoOutput->IsOpen())
 275    {
 276      cout<<"File is not opened"<<endl;
 277      return 1;
 278    }
 279
 280
 281   TTree *tracktree = new TTree("data","Tree with tracks");
 282
 283   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
 284   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
 285   tbuffer->SetOwner();
 286   pbuffer->SetOwner();
 287
 288   TClonesArray &particles = *pbuffer;
 289   TClonesArray &tracks = *tbuffer;
 290
 291   TString name("Tracks");
 292
 293
 294   Int_t NT = reader->GetNumberOfTrackEvents();
 295   Int_t NP = reader->GetNumberOfPartEvents();
 296
 297   Bool_t trck = (NT > 0) ? kTRUE : kFALSE;
 298   Bool_t part = (NP > 0) ? kTRUE : kFALSE;
 299
 300   TBranch *trackbranch = 0x0, *partbranch = 0x0;
 301
 302
 303   if (trck) trackbranch = tracktree->Branch("tracks","TClonesArray",&tbuffer);
 304   if (part) partbranch = tracktree->Branch("particles","TClonesArray",&pbuffer);
 305
 306
 307
 308   if ( (trck) && (part) && (NP != NT))
 309    {
 310      cerr<<"Warning number of track and particle events is different"<<endl;
 311    }
 312
 313   Int_t N;
 314   if (NT >= NP ) N = NT; else N = NP;
 315
 316   for ( i =0;i< N; i++)
 317     {
 318       if (trck && (i<=NT))
 319        {
 320          AliHBTEvent* trackev = reader->GetTrackEvent(i);
 321          for ( j = 0; j< trackev->GetNumberOfParticles();j++)
 322           {
 323             cout<<j<<"\r";
 324             new (tracks[j]) AliHBTParticle(*(trackev->GetParticle(j)));
 325           }
 326
 327        }
 328       cout<<endl;
 329
 330       if (part && (i<=NP))
 331        {
 332         AliHBTEvent* partev = reader->GetParticleEvent(i);
 333         for ( j = 0; j< partev->GetNumberOfParticles();j++)
 334          {
 335            cout<<j<<"\r";
 336            new (particles[j]) AliHBTParticle(*(partev->GetParticle(j)));
 337          }
 338
 339        }
 340
 341       histoOutput->cd();
 342       tracktree->Fill();
 343       tbuffer->Delete();
 344       pbuffer->Delete();
 345
 346      }
 347
 348
 349
 350   histoOutput->cd();
 351   tracktree->Write(0,TObject::kOverwrite);
 352   histoOutput->Close();
 353
 354   return 0;
 355  }