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