HBTAN/AliHBTReaderInternal.cxx

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