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 AliHBTReaderInternal test;
  16
  17 ClassImp(AliHBTReaderInternal)
  18 /********************************************************************/
  19
  20 AliHBTReaderInternal::AliHBTReaderInternal():
  21  fFileName(),
  22  fPartBranch(0x0),
  23  fTrackBranch(0x0),
  24  fTree(0x0),
  25  fFile(0x0)
  26 {
  27 //Defalut constructor
  28 }
  29 /********************************************************************/
  30
  31 AliHBTReaderInternal::AliHBTReaderInternal(const char *filename):
  32  fFileName(filename),
  33  fPartBranch(0x0),
  34  fTrackBranch(0x0),
  35  fTree(0x0),
  36  fFile(0x0)
  37 {
  38 //constructor
  39 //filename - name of file to open
  40 }
  41 /********************************************************************/
  42
  43 AliHBTReaderInternal::AliHBTReaderInternal(TObjArray* dirs, const char *filename):
  44  AliHBTReader(dirs),
  45  fFileName(filename),
  46  fPartBranch(0x0),
  47  fTrackBranch(0x0),
  48  fTree(0x0),
  49  fFile(0x0)
  50 {
  51 //ctor
  52 //dirs contains strings with directories to look data in
  53 //filename - name of file to open
  54 }
  55 /********************************************************************/
  56
  57 AliHBTReaderInternal::~AliHBTReaderInternal()
  58  {
  59  //desctructor
  60    delete fTree;
  61    delete fFile;
  62  }
  63 /********************************************************************/
  64
  65 void AliHBTReaderInternal::Rewind()
  66 {
  67   delete fTree;
  68   fTree = 0x0;
  69   delete fFile;
  70   fFile = 0x0;
  71   fCurrentDir = 0;
  72   fNEventsRead= 0;
  73 }
  74 /********************************************************************/
  75
  76 Int_t AliHBTReaderInternal::ReadNext()
  77  {
  78  //reads data and puts put to the particles and tracks objects
  79  //reurns 0 if everything is OK
  80  //
  81   Info("ReadNext","");
  82   Int_t i; //iterator and some temprary values
  83   Int_t counter;
  84   AliHBTParticle* tpart = 0x0, *ttrack = 0x0;
  85
  86   TDatabasePDG* pdgdb = TDatabasePDG::Instance();
  87   if (pdgdb == 0x0)
  88    {
  89      Error("ReadNext","Can not get PDG Particles Data Base");
  90      return 1;
  91    }
  92
  93   if (fParticlesEvent == 0x0)  fParticlesEvent = new AliHBTEvent();
  94   if (fTracksEvent == 0x0)  fTracksEvent = new AliHBTEvent();
  95
  96   fParticlesEvent->Reset();
  97   fTracksEvent->Reset();
  98
  99   do  //do{}while; is OK even if 0 dirs specified. In that case we try to read from "./"
 100    {
 101     if (fTree == 0x0)
 102      if( (i=OpenNextFile()) )
 103       {
 104         Error("ReadNext","Skipping directory due to problems with opening files. Errorcode %d",i);
 105         fCurrentDir++;
 106         continue;
 107       }
 108     if (fCurrentEvent == (Int_t)fTree->GetEntries())
 109      {
 110        delete fTree;
 111        fTree = 0x0;
 112        delete fFile;
 113        fFile = 0x0;
 114        fPartBranch = 0x0;
 115        fTrackBranch= 0x0;
 116        fCurrentDir++;
 117        continue;
 118      }
 119    /***************************/
 120    /***************************/
 121    /***************************/
 122
 123
 124     Info("ReadNext","Event %d",fCurrentEvent);
 125     fTree->GetEvent(fCurrentEvent);
 126
 127     counter = 0;
 128     if (fPartBranch && fTrackBranch)
 129      {
 130        Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
 131        Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
 132        for(i = 0; i < fPartBuffer->GetEntries(); i++)
 133          {
 134            tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
 135            ttrack =  dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
 136
 137            if( tpart == 0x0 ) continue; //if returned pointer is NULL
 138            if( ttrack == 0x0 ) continue; //if returned pointer is NULL
 139
 140            if (AliHBTParticle::GetDebug() > 9)
 141             {
 142               Info("ReadNext","Particle:");
 143               tpart->Print();
 144               Info("ReadNext","Track:");
 145               ttrack->Print();
 146             }
 147            if (ttrack->GetUID() != tpart->GetUID())
 148              {
 149                Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
 150                Error("ReadNext","They probobly do not correspond to each other.");
 151              }
 152
 153            for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
 154             {
 155               //check if we are intersted with particles of this type
 156               //if not take next partilce
 157               if( Pass(ttrack->GetNthPid(s)) )
 158                {
 159                  if (AliHBTParticle::GetDebug() > 9)
 160                   Info("ReadNext","Track Incarnation %d did not pass PID cut.",ttrack->GetNthPid(s));
 161                  continue;
 162                }
 163               TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
 164               if (pdgp == 0x0)//PDG part corresponding to new incarnation
 165                {
 166                  Error("ReadNext","Particle code unknown to PDG DB.");
 167                  continue;
 168                }
 169
 170               AliHBTParticle* track = new AliHBTParticle(*ttrack);
 171
 172               //apart of setting PDG code of an incarnation
 173               //it is necessary tu recalculate energy on the basis of
 174               //new PDG code (mass) hypothesis
 175               Double_t mass = pdgp->Mass();//mass of new incarnation
 176               Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
 177                                             ttrack->Py()*ttrack->Py() +
 178                                             ttrack->Pz()*ttrack->Pz() +
 179                                             mass*mass);//total energy of the new incarnation
 180               //update Energy and Calc Mass
 181               track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
 182               track->SetCalcMass(mass);
 183               track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
 184
 185               if( Pass(track) )
 186                 {
 187                   if (AliHBTParticle::GetDebug() > 9)
 188                    Info("ReadNext","Track Incarnation %d did not pass cut.",ttrack->GetNthPid(s));
 189                   delete track;
 190                   continue;
 191                 }
 192               AliHBTParticle* part = new AliHBTParticle(*tpart);//particle has only 1 incarnation (real)
 193
 194               fParticlesEvent->AddParticle(part);
 195               fTracksEvent->AddParticle(track);
 196
 197               counter++;
 198             }
 199          }
 200         Info("ReadNext","   Read: %d particles and tracks.",counter);
 201     }
 202    else
 203     {
 204      if (fPartBranch)
 205       {
 206         Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
 207         for(i = 0; i < fPartBuffer->GetEntries(); i++)
 208          {
 209            tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
 210            if(tpart == 0x0) continue; //if returned pointer is NULL
 211
 212            for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
 213             {
 214               if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
 215               if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
 216               AliHBTParticle* part = new AliHBTParticle(*tpart);
 217               part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
 218               if( Pass(part) )
 219                 {
 220                   delete part;
 221                   continue;
 222                 }
 223               fParticlesEvent->AddParticle(part);
 224               counter++;
 225             }
 226          }
 227         Info("ReadNext","   Read: %d particles.",counter);
 228       }
 229      else Info("ReadNext","   Read: 0 particles.");
 230
 231      if (fTrackBranch)
 232       {
 233         Info("ReadNext","Found %d tracks in total.",fTrackBuffer->GetEntries());
 234         for(i = 0; i < fTrackBuffer->GetEntries(); i++)
 235          {
 236            ttrack =  dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
 237            if( ttrack == 0x0 ) continue; //if returned pointer is NULL
 238
 239            for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
 240             {
 241               if( Pass(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
 242                                                          //if not take next partilce
 243               TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
 244               if (pdgp == 0x0)//PDG part corresponding to new incarnation
 245                {
 246                  Error("ReadNext","Particle code unknown to PDG DB.");
 247                  continue;
 248                }
 249               AliHBTParticle* track = new AliHBTParticle(*ttrack);
 250
 251               //apart of setting PDG code of an incarnation
 252               //it is necessary tu recalculate energy on the basis of
 253               //new PDG code (mass) hypothesis
 254               Double_t mass = pdgp->Mass();//mass of new incarnation
 255               Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
 256                                             ttrack->Py()*ttrack->Py() +
 257                                             ttrack->Pz()*ttrack->Pz() +
 258                                             mass*mass);//total energy of the new incarnation
 259               //update Energy and Calc Mass
 260               track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
 261               track->SetCalcMass(mass);
 262               track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
 263
 264               if( Pass(track) )
 265                 {
 266                   delete track;
 267                   continue;
 268                 }
 269               fTracksEvent->AddParticle(track);
 270
 271               counter++;
 272             }
 273          }
 274         Info("ReadNext","   Read: %d tracks",counter);
 275       }
 276      else Info("ReadNext","   Read: 0 tracks.");
 277     }
 278
 279     if (fPartBuffer) fPartBuffer->Delete();
 280     if (fTrackBuffer) fTrackBuffer->Delete();
 281
 282     fCurrentEvent++;
 283     fNEventsRead++;
 284
 285     return 0;
 286
 287    }while(fCurrentDir < GetNumberOfDirs());
 288
 289   return 1;//no more directories to read
 290 }
 291 /********************************************************************/
 292
 293 Int_t AliHBTReaderInternal::OpenNextFile()
 294 {
 295   //open file in current directory
 296    const TString& dirname = GetDirName(fCurrentDir);
 297    if (dirname == "")
 298     {
 299       Error("OpenNextFile","Can not get directory name");
 300       return 4;
 301     }
 302
 303    TString filename = dirname +"/"+ fFileName;
 304    fFile = TFile::Open(filename.Data());
 305    if ( fFile  == 0x0 )
 306      {
 307        Error("OpenNextFile","Can't open file named %s",filename.Data());
 308        return 1;
 309      }
 310    if (fFile->IsOpen() == kFALSE)
 311      {
 312        Error("OpenNextFile","Can't open filenamed %s",filename.Data());
 313        return 1;
 314      }
 315
 316    fTree = (TTree*)fFile->Get("data");
 317    if (fTree == 0x0)
 318     {
 319      Error("OpenNextFile","Can not get the tree.");
 320      return 1;
 321     }
 322
 323
 324    fTrackBranch = fTree->GetBranch("tracks");//get the branch with tracks
 325    if (fTrackBranch == 0x0) ////check if we got the branch
 326      {//if not return with error
 327        Info("OpenNextFile","Can't find a branch with tracks !\n");
 328      }
 329    else
 330     {
 331       if (fTrackBuffer == 0x0)
 332         fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
 333       fTrackBranch->SetAddress(&fTrackBuffer);
 334     }
 335
 336    fPartBranch = fTree->GetBranch("particles");//get the branch with particles
 337    if (fPartBranch == 0x0) ////check if we got the branch
 338      {//if not return with error
 339        Info("OpenNextFile","Can't find a branch with particles !\n");
 340      }
 341    else
 342     {
 343       if (fPartBuffer == 0x0)
 344         fPartBuffer = new TClonesArray("AliHBTParticle",15000);
 345       fPartBranch->SetAddress(&fPartBuffer);
 346     }
 347
 348    Info("OpenNextFile","________________________________________________________");
 349    Info("OpenNextFile","Found %d event(s) in directory %s",
 350          (Int_t)fTree->GetEntries(),GetDirName(fCurrentDir).Data());
 351
 352    fCurrentEvent = 0;
 353
 354    return 0;
 355 }
 356 /********************************************************************/
 357
 358 Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile, Bool_t multcheck)
 359  {
 360   //reads tracks from reader and writes runs to file
 361   //reader - provides data for writing in internal format
 362   //name of output file
 363   //multcheck - switches of checking if particle was stored with other incarnation
 364   // usefull e.g. when using kine data, where all particles have 100% pid prob.and saves a lot of time
 365
 366   Int_t i,j;
 367
 368   ::Info("AliHBTReaderInternal::Write","________________________________________________________");
 369   ::Info("AliHBTReaderInternal::Write","________________________________________________________");
 370   ::Info("AliHBTReaderInternal::Write","________________________________________________________");
 371
 372   TFile *histoOutput = TFile::Open(outfile,"recreate");
 373
 374   if (!histoOutput->IsOpen())
 375    {
 376      ::Error("AliHBTReaderInternal::Write","File is not opened");
 377      return 1;
 378    }
 379
 380   TTree *tracktree = new TTree("data","Tree with tracks");
 381
 382   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
 383   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
 384
 385   TClonesArray &particles = *pbuffer;
 386   TClonesArray &tracks = *tbuffer;
 387
 388   TString name("Tracks");
 389
 390   Int_t nt = reader->GetNumberOfTrackEvents();
 391   Int_t np = reader->GetNumberOfPartEvents();
 392
 393   if (AliHBTParticle::GetDebug() > 0)
 394    ::Info("Write","Reader has %d track events and %d particles events.",nt,np);
 395
 396   Bool_t trck = (nt > 0) ? kTRUE : kFALSE;
 397   Bool_t part = (np > 0) ? kTRUE : kFALSE;
 398
 399   TBranch *trackbranch = 0x0, *partbranch = 0x0;
 400
 401   if (trck) trackbranch = tracktree->Branch("tracks",&tbuffer,32000,0);
 402   if (part) partbranch = tracktree->Branch("particles",&pbuffer,32000,0);
 403
 404   if ( (trck) && (part) && (np != nt))
 405    {
 406      ::Warning("AliHBTReaderInternal::Write","Number of track and particle events is different");
 407    }
 408
 409   Int_t n;
 410   if (nt >= np ) n = nt; else n = np;
 411
 412   if (AliHBTParticle::GetDebug() > 0)
 413    ::Info("Write","Will loop over %d events",n);
 414
 415   for ( i =0;i< n; i++)
 416     {
 417       ::Info("AliHBTReaderInternal::Write","Event %d",i+1);
 418       Int_t counter = 0;
 419       if (trck && (i<=nt))
 420        {
 421          AliHBTEvent* trackev = reader->GetTrackEvent(i);
 422          for ( j = 0; j< trackev->GetNumberOfParticles();j++)
 423           {
 424             const AliHBTParticle& t = *(trackev->GetParticle(j));
 425             if (multcheck)
 426              {
 427               if (FindIndex(tbuffer,t.GetUID()))
 428                {
 429                  if (AliHBTParticle::GetDebug()>4)
 430                   {
 431                    ::Info("Write","Track with Event UID %d already stored",t.GetUID());
 432                   }
 433                  continue; //not to write the same particles with other incarnations
 434                }
 435              }
 436             new (tracks[counter++]) AliHBTParticle(t);
 437           }
 438          ::Info("AliHBTReaderInternal::Write","    Tracks: %d",tracks.GetEntries());
 439        }else ::Info("AliHBTReaderInternal::Write","NO TRACKS");
 440
 441       counter = 0;
 442       if (part && (i<=np))
 443        {
 444 //        ::Warning("AliHBTReaderInternal::Write","Find index switched off!!!");
 445
 446         AliHBTEvent* partev = reader->GetParticleEvent(i);
 447         for ( j = 0; j< partev->GetNumberOfParticles();j++)
 448          {
 449            const AliHBTParticle& part= *(partev->GetParticle(j));
 450             if (multcheck)
 451              {
 452               if (FindIndex(pbuffer,part.GetUID()))
 453                {
 454                  if (AliHBTParticle::GetDebug()>4)
 455                   {
 456                    ::Info("Write","Particle with Event UID %d already stored",part.GetUID());
 457                   }
 458                  continue; //not to write the same particles with other incarnations
 459                }
 460              }
 461            new (particles[counter++]) AliHBTParticle(part);
 462          }
 463          ::Info("AliHBTReaderInternal::Write","    Particles: %d",particles.GetEntries());
 464        }else ::Info("AliHBTReaderInternal::Write","NO PARTICLES");
 465
 466       histoOutput->cd();
 467       tracktree->Fill();
 468       tracktree->AutoSave();
 469       tbuffer->Delete();
 470       pbuffer->Delete();
 471      }
 472
 473   histoOutput->cd();
 474   tracktree->Write(0,TObject::kOverwrite);
 475   delete tracktree;
 476
 477   tbuffer->SetOwner();
 478   pbuffer->SetOwner();
 479   delete pbuffer;
 480   delete tbuffer;
 481
 482   histoOutput->Close();
 483   return 0;
 484  }
 485 /********************************************************************/
 486
 487 Bool_t AliHBTReaderInternal::FindIndex(TClonesArray* arr,Int_t idx)
 488 {
 489 //Checks if in the array exists already partilce with Unique ID idx
 490   if (arr == 0x0)
 491    {
 492      ::Error("FindIndex","Array is 0x0");
 493      return kTRUE;
 494    }
 495   TIter next(arr);
 496   AliHBTParticle* p;
 497   while (( p = (AliHBTParticle*)next()))
 498    {
 499      if (p->GetUID() == idx) return kTRUE;
 500    }
 501   return kFALSE;
 502 }