HBTAN/AliHBTReaderInternal.cxx

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