HBTAN/AliHBTReaderInternal.cxx

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