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        for(i = 0; i < fPartBuffer->GetEntries(); i++)
 131          {
 132            tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
 133            ttrack =  dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
 134
 135            if( tpart == 0x0 ) continue; //if returned pointer is NULL
 136            if( ttrack == 0x0 ) continue; //if returned pointer is NULL
 137
 138            if (ttrack->GetUID() != tpart->GetUID())
 139              {
 140                Error("ReadNext","Sth. is wrong: Track and Particle has different UID.");
 141                Error("ReadNext","They probobly do not correspond to each other.");
 142              }
 143
 144            for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
 145             {
 146               if( Pass(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
 147                                                          //if not take next partilce
 148               TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
 149               if (pdgp == 0x0)//PDG part corresponding to new incarnation
 150                {
 151                  Error("ReadNext","Particle code unknown to PDG DB.");
 152                  continue;
 153                }
 154
 155               AliHBTParticle* track = new AliHBTParticle(*ttrack);
 156
 157               //apart of setting PDG code of an incarnation
 158               //it is necessary tu recalculate energy on the basis of
 159               //new PDG code (mass) hypothesis
 160               Double_t mass = pdgp->Mass();//mass of new incarnation
 161               Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
 162                                             ttrack->Py()*ttrack->Py() +
 163                                             ttrack->Pz()*ttrack->Pz() +
 164                                             mass*mass);//total energy of the new incarnation
 165               //update Energy and Calc Mass
 166               track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
 167               track->SetCalcMass(mass);
 168               track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
 169
 170               if( Pass(track) )
 171                 {
 172                   delete track;
 173                   continue;
 174                 }
 175               AliHBTParticle* part = new AliHBTParticle(*tpart);//particle has only 1 incarnation (real)
 176
 177               fParticlesEvent->AddParticle(part);
 178               fTracksEvent->AddParticle(track);
 179
 180               counter++;
 181             }
 182          }
 183         Info("ReadNext","   Read: %d particles and tracks.",counter);
 184     }
 185    else
 186     {
 187      if (fPartBranch)
 188       {
 189         Info("ReadNext","Found %d particles in total.",fPartBuffer->GetEntries());
 190         for(i = 0; i < fPartBuffer->GetEntries(); i++)
 191          {
 192            tpart = dynamic_cast<AliHBTParticle*>(fPartBuffer->At(i));
 193            if(tpart == 0x0) continue; //if returned pointer is NULL
 194
 195            for (Int_t s = 0; s < tpart->GetNumberOfPids(); s++)
 196             {
 197               if( pdgdb->GetParticle(tpart->GetNthPid(s)) == 0x0 ) continue; //if particle has crazy PDG code (not known to our database)
 198               if( Pass(tpart->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
 199               AliHBTParticle* part = new AliHBTParticle(*tpart);
 200               part->SetPdgCode(tpart->GetNthPid(s),tpart->GetNthPidProb(s));
 201               if( Pass(part) )
 202                 {
 203                   delete part;
 204                   continue;
 205                 }
 206               fParticlesEvent->AddParticle(part);
 207               counter++;
 208             }
 209          }
 210         Info("ReadNext","   Read: %d particles.",counter);
 211       }
 212      else Info("ReadNext","   Read: 0 particles.");
 213
 214      if (fTrackBranch)
 215       {
 216         for(i = 0; i < fTrackBuffer->GetEntries(); i++)
 217          {
 218            ttrack =  dynamic_cast<AliHBTParticle*>(fTrackBuffer->At(i));
 219            if( ttrack == 0x0 ) continue; //if returned pointer is NULL
 220
 221            for (Int_t s = 0; s < ttrack->GetNumberOfPids(); s++)
 222             {
 223               if( Pass(ttrack->GetNthPid(s)) ) continue; //check if we are intersted with particles of this type
 224                                                          //if not take next partilce
 225               TParticlePDG* pdgp = pdgdb->GetParticle(ttrack->GetNthPid(s));
 226               if (pdgp == 0x0)//PDG part corresponding to new incarnation
 227                {
 228                  Error("ReadNext","Particle code unknown to PDG DB.");
 229                  continue;
 230                }
 231               AliHBTParticle* track = new AliHBTParticle(*ttrack);
 232
 233               //apart of setting PDG code of an incarnation
 234               //it is necessary tu recalculate energy on the basis of
 235               //new PDG code (mass) hypothesis
 236               Double_t mass = pdgp->Mass();//mass of new incarnation
 237               Double_t tEtot = TMath::Sqrt( ttrack->Px()*ttrack->Px() +
 238                                             ttrack->Py()*ttrack->Py() +
 239                                             ttrack->Pz()*ttrack->Pz() +
 240                                             mass*mass);//total energy of the new incarnation
 241               //update Energy and Calc Mass
 242               track->SetMomentum(ttrack->Px(),ttrack->Py(),ttrack->Pz(),tEtot);
 243               track->SetCalcMass(mass);
 244               track->SetPdgCode(ttrack->GetNthPid(s),ttrack->GetNthPidProb(s));
 245
 246               if( Pass(track) )
 247                 {
 248                   delete track;
 249                   continue;
 250                 }
 251               fTracksEvent->AddParticle(track);
 252
 253               counter++;
 254             }
 255          }
 256         Info("ReadNext","   Read: %d tracks",counter);
 257       }
 258      else Info("ReadNext","   Read: 0 tracks.");
 259     }
 260
 261     if (fPartBuffer) fPartBuffer->Delete();
 262     if (fTrackBuffer) fTrackBuffer->Delete();
 263
 264     fCurrentEvent++;
 265     fNEventsRead++;
 266
 267     return 0;
 268
 269    }while(fCurrentDir < GetNumberOfDirs());
 270
 271   return 1;//no more directories to read
 272 }
 273 /********************************************************************/
 274
 275 Int_t AliHBTReaderInternal::OpenNextFile()
 276 {
 277   //open file in current directory
 278    const TString& dirname = GetDirName(fCurrentDir);
 279    if (dirname == "")
 280     {
 281       Error("OpenNextFile","Can not get directory name");
 282       return 4;
 283     }
 284
 285    TString filename = dirname +"/"+ fFileName;
 286    fFile = TFile::Open(filename.Data());
 287    if ( fFile  == 0x0 )
 288      {
 289        Error("OpenNextFile","Can't open file named %s",filename.Data());
 290        return 1;
 291      }
 292    if (fFile->IsOpen() == kFALSE)
 293      {
 294        Error("OpenNextFile","Can't open filenamed %s",filename.Data());
 295        return 1;
 296      }
 297
 298    fTree = (TTree*)fFile->Get("data");
 299    if (fTree == 0x0)
 300     {
 301      Error("OpenNextFile","Can not get the tree.");
 302      return 1;
 303     }
 304
 305
 306    fTrackBranch = fTree->GetBranch("tracks");//get the branch with tracks
 307    if (fTrackBranch == 0x0) ////check if we got the branch
 308      {//if not return with error
 309        Info("OpenNextFile","Can't find a branch with tracks !\n");
 310      }
 311    else
 312     {
 313       if (fTrackBuffer == 0x0)
 314         fTrackBuffer = new TClonesArray("AliHBTParticle",15000);
 315       fTrackBranch->SetAddress(&fTrackBuffer);
 316     }
 317
 318    fPartBranch = fTree->GetBranch("particles");//get the branch with particles
 319    if (fPartBranch == 0x0) ////check if we got the branch
 320      {//if not return with error
 321        Info("OpenNextFile","Can't find a branch with particles !\n");
 322      }
 323    else
 324     {
 325       if (fPartBuffer == 0x0)
 326         fPartBuffer = new TClonesArray("AliHBTParticle",15000);
 327       fPartBranch->SetAddress(&fPartBuffer);
 328     }
 329
 330    Info("OpenNextFile","________________________________________________________");
 331    Info("OpenNextFile","Found %d event(s) in directory %s",
 332          (Int_t)fTree->GetEntries(),GetDirName(fCurrentEvent).Data());
 333
 334    fCurrentEvent = 0;
 335
 336    return 0;
 337 }
 338 /********************************************************************/
 339
 340 Int_t AliHBTReaderInternal::Write(AliHBTReader* reader,const char* outfile, Bool_t multcheck)
 341  {
 342   //reads tracks from reader and writes runs to file
 343   //reader - provides data for writing in internal format
 344   //name of output file
 345   //multcheck - switches of checking if particle was stored with other incarnation
 346   // usefull e.g. when using kine data, where all particles have 100% pid prob.and saves a lot of time
 347
 348   Int_t i,j;
 349
 350   ::Info("AliHBTReaderInternal::Write","________________________________________________________");
 351   ::Info("AliHBTReaderInternal::Write","________________________________________________________");
 352   ::Info("AliHBTReaderInternal::Write","________________________________________________________");
 353
 354   TFile *histoOutput = TFile::Open(outfile,"recreate");
 355
 356   if (!histoOutput->IsOpen())
 357    {
 358      ::Error("AliHBTReaderInternal::Write","File is not opened");
 359      return 1;
 360    }
 361
 362   TTree *tracktree = new TTree("data","Tree with tracks");
 363
 364   TClonesArray* pbuffer = new TClonesArray("AliHBTParticle",15000);
 365   TClonesArray* tbuffer = new TClonesArray("AliHBTParticle",15000);
 366
 367   TClonesArray &particles = *pbuffer;
 368   TClonesArray &tracks = *tbuffer;
 369
 370   TString name("Tracks");
 371
 372   Int_t nt = reader->GetNumberOfTrackEvents();
 373   Int_t np = reader->GetNumberOfPartEvents();
 374
 375   if (AliHBTParticle::GetDebug() > 0)
 376    ::Info("Write","Reader has %d track events and %d particles events.",nt,np);
 377
 378   Bool_t trck = (nt > 0) ? kTRUE : kFALSE;
 379   Bool_t part = (np > 0) ? kTRUE : kFALSE;
 380
 381   TBranch *trackbranch = 0x0, *partbranch = 0x0;
 382
 383   if (trck) trackbranch = tracktree->Branch("tracks",&tbuffer,32000,0);
 384   if (part) partbranch = tracktree->Branch("particles",&pbuffer,32000,0);
 385
 386   if ( (trck) && (part) && (np != nt))
 387    {
 388      ::Warning("AliHBTReaderInternal::Write","Number of track and particle events is different");
 389    }
 390
 391   Int_t n;
 392   if (nt >= np ) n = nt; else n = np;
 393
 394   if (AliHBTParticle::GetDebug() > 0)
 395    ::Info("Write","Will loop over %d events",n);
 396
 397   for ( i =0;i< n; i++)
 398     {
 399       ::Info("AliHBTReaderInternal::Write","Event %d",i+1);
 400       Int_t counter = 0;
 401       if (trck && (i<=nt))
 402        {
 403          AliHBTEvent* trackev = reader->GetTrackEvent(i);
 404          for ( j = 0; j< trackev->GetNumberOfParticles();j++)
 405           {
 406             const AliHBTParticle& t = *(trackev->GetParticle(j));
 407             if (multcheck)
 408              {
 409               if (FindIndex(tbuffer,t.GetUID()))
 410                {
 411                  if (AliHBTParticle::GetDebug()>4)
 412                   {
 413                    ::Info("Write","Track with Event UID %d already stored",t.GetUID());
 414                   }
 415                  continue; //not to write the same particles with other incarnations
 416                }
 417              }
 418             new (tracks[counter++]) AliHBTParticle(t);
 419           }
 420          ::Info("AliHBTReaderInternal::Write","    Tracks: %d",tracks.GetEntries());
 421        }else ::Info("AliHBTReaderInternal::Write","NO TRACKS");
 422
 423       counter = 0;
 424       if (part && (i<=np))
 425        {
 426 //        ::Warning("AliHBTReaderInternal::Write","Find index switched off!!!");
 427
 428         AliHBTEvent* partev = reader->GetParticleEvent(i);
 429         for ( j = 0; j< partev->GetNumberOfParticles();j++)
 430          {
 431            const AliHBTParticle& part= *(partev->GetParticle(j));
 432             if (multcheck)
 433              {
 434               if (FindIndex(pbuffer,part.GetUID()))
 435                {
 436                  if (AliHBTParticle::GetDebug()>4)
 437                   {
 438                    ::Info("Write","Particle with Event UID %d already stored",part.GetUID());
 439                   }
 440                  continue; //not to write the same particles with other incarnations
 441                }
 442              }
 443            new (particles[counter++]) AliHBTParticle(part);
 444          }
 445          ::Info("AliHBTReaderInternal::Write","    Particles: %d",particles.GetEntries());
 446        }else ::Info("AliHBTReaderInternal::Write","NO PARTICLES");
 447
 448       histoOutput->cd();
 449       tracktree->Fill();
 450       tracktree->AutoSave();
 451       tbuffer->Delete();
 452       pbuffer->Delete();
 453      }
 454
 455   histoOutput->cd();
 456   tracktree->Write(0,TObject::kOverwrite);
 457   delete tracktree;
 458
 459   tbuffer->SetOwner();
 460   pbuffer->SetOwner();
 461   delete pbuffer;
 462   delete tbuffer;
 463
 464   histoOutput->Close();
 465   return 0;
 466  }
 467 /********************************************************************/
 468
 469 Bool_t AliHBTReaderInternal::FindIndex(TClonesArray* arr,Int_t idx)
 470 {
 471 //Checks if in the array exists already partilce with Unique ID idx
 472   if (arr == 0x0)
 473    {
 474      ::Error("FindIndex","Array is 0x0");
 475      return kTRUE;
 476    }
 477   TIter next(arr);
 478   AliHBTParticle* p;
 479   while (( p = (AliHBTParticle*)next()))
 480    {
 481      if (p->GetUID() == idx) return kTRUE;
 482    }
 483   return kFALSE;
 484 }