Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / THbtp / AliGenHBTprocessor.cxx
index b192097..6c7b8bf 100644 (file)
 //         A)  HBT PROCESSOR NEEDS MORE THAN ONE EVENT TO WORK
 //             AS MORE AS IT BETTER WORKS
 //         B)  IT IS ABLE TO "ADD" CORRELATIONS ONLY UP TO TWO PARTICLE TYPES AT ONES
+//
+//
+// Artificial particle dennsity enhancment feature
+// HBT Processor is unable to process correctly low multiplicity particles
+// even if the high statistics (more events) is supplied.
+// For that reason was implemented artificial multiplicity enhancement 
+// - see also comments in HBT Processor source file (HBTP/hbt_event_processor.f)
+//   or web page (http://alisoft.cern.ch/people/skowron/hbtprocessor/hbteventprocessor.html) 
+//   concerning trk_accep or  SetTrackRejectionFactor method in this class
+// Fortran is cheated by masking several events a single one. Number is defined by fEventMerge
+// variable. In case it is equal to 1, there is no masking and the feature is off.
+// But, for example if it is 5, and we have 100 events, number of events passed to fortran is 100/5=20
+// When fortran asks about multiplcity of event, let sey, 1 - multiplicity of events 5 to 9 is summed up 
+// and returned. 
+//
+//
 //////////////////////////////////////////////////////////////////////////////////
 
 // 11.11.2001 Piotr Skowronski
 // 
 
 #include "AliGenHBTprocessor.h"
-#include "TROOT.h"
-#include <Riostream.h>
-#include <TFile.h>
-#include <TTree.h>
-#include "AliRun.h"
-#include "AliStack.h"
-#include "TParticle.h"
+
+#include <TParticle.h>
 #include "THBTprocessor.h"
+
+#include "AliStack.h"
+#include "AliMC.h"
+#include "AliRun.h"
 #include "AliGenCocktailAfterBurner.h"
+#include "AliLog.h"
 
 
 
 ClassImp(AliGenHBTprocessor)
 
-
+Int_t AliGenHBTprocessor::fgDebug = 0;
+static TRandom* gAliRandom;//RNG to be used by the fortran code
 
 AliGenCocktailAfterBurner*  GetGenerator();
 /*******************************************************************/
-AliGenHBTprocessor::AliGenHBTprocessor(const AliGenHBTprocessor& in)
-{
-//copy contructor
-  // AliGenHBTprocessor::AliGenHBTprocessor();
-}
 
-AliGenHBTprocessor::AliGenHBTprocessor() : AliGenerator() 
+AliGenHBTprocessor::AliGenHBTprocessor(): 
+  AliGenerator(), 
+  fHBTprocessor(0x0),
+  fHbtPStatCodes(0x0),
+  fNPDGCodes(0),
+  fTrackRejectionFactor(0),
+  fReferenceControl(0),
+  fPrintFull(0),
+  fPrintSectorData(0),
+  fNPidTypes(0),
+  fNevents(0),
+  fSwitch1d(0),
+  fSwitch3d(0),
+  fSwitchType(0),
+  fSwitchCoherence(0),
+  fSwitchCoulomb(0),
+  fSwitchFermiBose(0),
+  fEventLineCounter(0),
+  fMaxit(0),
+  fIrand(0),
+  fLambda(0),
+  fR1d(0),
+  fRside(0),
+  fRout(0),
+  fRlong(0),
+  fRperp(0),
+  fRparallel(0),
+  fR0(0),
+  fQ0(0),
+  fDeltap(0),
+  fDelchi(0),
+  fNPtBins(0),
+  fNPhiBins(0),
+  fNEtaBins(0),
+  fN1dFine(0),
+  fN1dCoarse(0),
+  fN1dTotal(0),
+  fN3dFine(0),
+  fN3dCoarse(0),
+  fN3dTotal(0),
+  fN3dFineProject(0),
+  fNPxBins(0),
+  fNPyBins(0),
+  fNPzBins(0),
+  fNSectors(0),
+  fPtBinSize(0),
+  fPhiBinSize(0),
+  fEtaBinSize(0),
+  fEtaMin(0),
+  fEtaMax(0),
+  fBinsize1dFine(0),
+  fBinsize1dCoarse(0),
+  fQmid1d(0),
+  fQmax1d(0),
+  fBinsize3dFine(0),
+  fBinsize3dCoarse(0),
+  fQmid3d(0),
+  fQmax3d(0),
+  fPxMin(0),
+  fPxMax(0),
+  fDelpx(0),
+  fPyMin(0),
+  fPyMax(0),
+  fDelpy(0),
+  fPzMin(0),
+  fPzMax(0),
+  fDelpz(0),
+  fEventMerge(1),
+  fActiveStack(0)
 {
   //
   // Standard constructor
   // Sets default veues of all parameters
-  fHbtPStatCodes = 0x0;
-  fHBTprocessor = 0x0;
 
   SetName("AliGenHBTprocessor");
   SetTitle("AliGenHBTprocessor");
   
-  sRandom = fRandom;
+  gAliRandom = fRandom;
   fHBTprocessor = new THBTprocessor();
 
   fNPDGCodes = 0;
@@ -163,7 +241,10 @@ void AliGenHBTprocessor::InitStatusCodes()
  //creates and inits status codes array to zero
   AliGenCocktailAfterBurner *cab = GetGenerator();
 
-  if(!cab) Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
+  if(!cab) {
+    Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
+    return;
+  }
 
   Int_t nev = cab->GetNumberOfEvents();
 
@@ -187,7 +268,7 @@ void AliGenHBTprocessor::CleanStatusCodes()
   {
     for (Int_t i =0; i<GetGenerator()->GetNumberOfEvents(); i++)
       delete [] fHbtPStatCodes[i];  
-    delete fHbtPStatCodes;
+    delete [] fHbtPStatCodes;
     fHbtPStatCodes = 0;
   }
 
@@ -248,7 +329,8 @@ void AliGenHBTprocessor::Init()
    thbtp->SetPxRange(fPxMin,fPxMax);
    thbtp->SetPyRange(fPyMin,fPyMax);
    thbtp->SetPzRange(fPzMin,fPzMax);
-   thbtp->SetPhiRange(fPhiMin*180./TMath::Pi(),fPhiMax*180./TMath::Pi());
+   thbtp->SetPhiRange(fPhiMin*180./((Float_t)TMath::Pi())+180.0, //casting is because if fPhiMin = 180.0 then
+                      fPhiMax*180./((Float_t)TMath::Pi())+180.0);//TMath::Pi() != TMath::Pi()*fPhiMin/180.0,
    thbtp->SetEtaRange(fEtaMin,fEtaMax);
    thbtp->SetNPtBins(fNPtBins);
    thbtp->SetNPhiBins(fNPhiBins);
@@ -265,6 +347,8 @@ void AliGenHBTprocessor::Init()
    thbtp->SetNBins3DCoarseMesh(fN3dCoarse);
    thbtp->SetBinSize3DCoarseMesh(fBinsize3dCoarse);
    thbtp->SetNBins3DFineProjectMesh(fN3dFineProject);
+   
+   thbtp->SetPrintFull(fPrintFull);
        
  }
 /*******************************************************************/
@@ -276,6 +360,7 @@ void AliGenHBTprocessor::Generate()
    if (cab == 0x0)
     {
       Fatal("Generate()","AliGenHBTprocessor needs AliGenCocktailAfterBurner to be main generator");
+      return;
     }
    if (cab->GetNumberOfEvents() <2)
     {
@@ -307,7 +392,7 @@ void AliGenHBTprocessor::GetParticles(TClonesArray * particles) const
  //practically dumm
    if (!particles)
     {
-      cout<<"Particles has to be initialized"<<endl;
+      Error("GetParticles","Particles has to be initialized");
       return;
     } 
    fHBTprocessor->ImportParticles(particles);
@@ -320,16 +405,29 @@ Int_t AliGenHBTprocessor::GetHbtPStatusCode(Int_t part) const
 //returns the status code of the given particle in the active event
 //see SetActiveEvent in the bottom of AliGenHBTprocessor.cxx
 //and in AliCocktailAfterBurner
-  Int_t activeEvent = GetGenerator()->GetActiveEventNumber();
-  return fHbtPStatCodes[activeEvent][part];
+ Int_t ev, idx;
+ GetTrackEventIndex(part,ev,idx);
+ if ( (ev<0) || (idx<0) )
+  {
+    Error("GetHbtPStatusCode","GetTrackEventIndex returned error");
+    return 0;
+  }
+ return fHbtPStatCodes[ev][idx];
+  
 }
 
 /*******************************************************************/
 void  AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part)
 {
  //Sets the given status code to given particle number (part) in the active event
-  Int_t activeEvent = GetGenerator()->GetActiveEventNumber();
-  fHbtPStatCodes[activeEvent][part] = hbtstatcode;
+ Int_t ev, idx;
+ GetTrackEventIndex(part,ev,idx);
+ if ( (ev<0) || (idx<0) )
+  {
+    Error("GetHbtPStatusCode","GetTrackEventIndex returned error");
+    return;
+  }
+ else fHbtPStatCodes[ev][idx] = hbtstatcode;
 }
 
 /*******************************************************************/
@@ -629,7 +727,7 @@ void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax)
    fPzMax =pzmax; 
    fHBTprocessor->SetPzRange(pzmin,pzmax);
  }
-void AliGenHBTprocessor::SetMomentumRange(Float_t pmin, Float_t pmax)
+void AliGenHBTprocessor::SetMomentumRange(Float_t /*pmin*/, Float_t /*pmax*/)
  {
  //default pmin=0, pmax=0
  //Do not use this method!
@@ -643,7 +741,7 @@ void AliGenHBTprocessor::SetPhiRange(Float_t phimin, Float_t phimax)
 //Sets \\Phi range  
   AliGenerator::SetPhiRange(phimin,phimax);
   
-  fHBTprocessor->SetPhiRange(phimin,phimax);
+  fHBTprocessor->SetPhiRange(phimin+180.0,phimax+180.0);
  }
 /*******************************************************************/
 void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax)
@@ -803,10 +901,139 @@ void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n )
   fHBTprocessor->SetNBins3DFineProjectMesh(n);
 }
 /*******************************************************************/
+void AliGenHBTprocessor::SetPrintFull(Int_t flag)
+{
+//sets the print full flag
+ fPrintFull = flag;
+ fHBTprocessor->SetPrintFull(flag);
+}
+
+
+/*******************************************************************/
 
+Int_t  AliGenHBTprocessor::GetNumberOfEvents()
+{
+//returns number of available events
+  AliGenerator* g = gAlice->GetMCApp()->Generator();
+  AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
+  if (cab == 0x0)
+   {
+     Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner");
+     return 0;
+   }
+  return (Int_t)TMath::Ceil(cab->GetNumberOfEvents()/((Float_t)fEventMerge));
+}
 
 /*******************************************************************/
 
+void AliGenHBTprocessor::SetActiveEventNumber(Int_t n)
+{
+//sets the active event
+ fActiveStack =  n*fEventMerge;
+ AliDebug(1,Form("Settimg active event %d passed %d",fActiveStack,n));
+}
+/*******************************************************************/
+
+Int_t  AliGenHBTprocessor::GetNumberOfTracks()
+{
+//returns number of tracks in active event
+  AliGenerator* g = gAlice->GetMCApp()->Generator();
+  AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
+  if (cab == 0x0)
+   {
+     Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner");
+     return 0;
+   }
+ Int_t n = 0;
+ for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++) 
+  { 
+    if (i >= GetNumberOfEvents()) break; //protection not to overshoot nb of events
+    AliStack* stack = cab->GetStack(i);
+    if (stack == 0x0) {
+     Error("GetNumberOfTracks","There is no stack %d",i);
+     continue;
+    }
+    n+=stack->GetNprimary();
+  }
+ return n;
+}
+/*******************************************************************/
+
+void AliGenHBTprocessor::SetNEventsToMerge(Int_t nev)
+{
+ //sets number of events to merge
+ if (nev > 0 ) fEventMerge = nev;
+}
+/*******************************************************************/
+
+TParticle* AliGenHBTprocessor::GetTrack(Int_t n)
+{ 
+//returns track that hbtp thinks is n in active event
+  AliDebug(5,Form("n = %d",n));
+  AliGenerator* g = gAlice->GetMCApp()->Generator();
+  AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
+  if (cab == 0x0)
+   {
+     Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner");
+     return 0;
+   }
+
+ Int_t ev, idx;
+ GetTrackEventIndex(n,ev,idx);
+ AliDebug(5,Form("Event = %d Particle = %d",ev,idx));
+ if ( (ev<0) || (idx<0) )
+  {
+    Error("GetTrack","GetTrackEventIndex returned error");
+    return 0x0;
+  }
+ AliDebug(5,Form("Number of Tracks in Event(%d) = %d",ev,cab->GetStack(ev)->GetNprimary()));
+ return cab->GetStack(ev)->Particle(idx);//safe - in case stack does not exist 
+}
+/*******************************************************************/
+
+void AliGenHBTprocessor::GetTrackEventIndex(Int_t n, Int_t &evno, Int_t &index) const
+{
+ //returns event(stack) number and particle index
+  AliGenerator* g = gAlice->GetMCApp()->Generator();
+  AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
+  if (cab == 0x0)
+   {
+     Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner");
+     return;
+   }
+
+ evno = -1;
+ index = -1;
+ for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++) 
+  { 
+    AliStack* stack = cab->GetStack(i);
+    if (stack == 0x0)
+     {
+       Error("GetTrackEventIndex","There is no stack %d",i);
+       return;
+     }
+
+    Int_t ntracks = stack->GetNprimary();
+    AliDebug(10,Form("Event %d has %d tracks. n = %d",i,ntracks,n));
+    
+    if ( ntracks > n) 
+     {
+       evno = i;
+       index = n;
+       return ;
+     }
+    else 
+     {  
+       n-=ntracks;
+       continue;
+     }
+  }
+ Error("GetTrackEventIndex","Could not find given track");
+}
+
+/*******************************************************************/
+/*******************************************************************/
+/*******************************************************************/
 
 
 
@@ -874,18 +1101,18 @@ Int_t AliGenHBTprocessor::PDGFromId(Int_t id) const
 Double_t AliGenHBTprocessor::ThetaToEta(Double_t arg)
  {
   //converts etha(azimuthal angle) to Eta (pseudorapidity). Argument in radians
+
+   if(arg>= TMath::Pi()) return  708.0;//This number is the biggest wich not crashes exp(x)
+   if(arg<= 0.0) return -708.0;//
+   if(arg == TMath::Pi()/2.) return 0.0;//
    
-   if(arg>= TMath::Pi()) return  709.0;//This number is the biggest wich not crashes exp(x)
-   if(arg<= 0.0) return -709.0;//
-   
-   arg -= TMath::Pi()/2.;
    if (arg > 0.0) 
     { 
-      return -TMath::Log( TMath::Tan(arg/2.)) ;
+      return TMath::Log( TMath::Tan(arg/2.)) ;
     }
    else 
-    {
-      return TMath::Log( TMath::Tan(-arg/2.)) ;
+    { 
+      return -TMath::Log( TMath::Tan(-arg/2.)) ;
     }
  }
                                   
@@ -917,8 +1144,6 @@ Double_t AliGenHBTprocessor::ThetaToEta(Double_t arg)
   # define type_ofCall  _stdcall
 #endif    
 
-#include "AliGenCocktailAfterBurner.h"
-#include <string.h>
 /*******************************************************************/
 
 AliGenCocktailAfterBurner*  GetGenerator()
@@ -930,17 +1155,18 @@ AliGenCocktailAfterBurner*  GetGenerator()
 
    if(!gAlice)
     {
-      cout<<endl<<"ERROR: There is NO gALICE! Check what you are doing!"<<endl;
-      gROOT->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
-                    "\nRunning HBT Processor without gAlice... Exiting \n");
-      return 0x0;
+      ::Error("AliGenHBTprocessor.cxx: GetGenerator()",
+              "There is NO gALICE! Check what you are doing!");
+      ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
+              "Running HBT Processor without gAlice... Exiting \n");
+      return 0x0;//pro forma
     }
-   AliGenerator * gen = gAlice->Generator();
+   AliGenerator * gen = gAlice->GetMCApp()->Generator();
    
    if (!gen) 
     {
-      gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
-                   "\nThere is no generator in gAlice, exiting\n");
+      ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
+               "There is no generator in gAlice, exiting\n");
       return 0x0;
     }
 
@@ -955,13 +1181,11 @@ AliGenCocktailAfterBurner*  GetGenerator()
                                                                         
    if (cab == 0x0)//if generator that we got is not AliGenCocktailAfterBurner or its descendant we get null
    {              //then quit with error
-      gAlice->Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
-                    "\nThe main Generator is not a AliGenCocktailAfterBurner, exiting\n");
+      ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
+              "The main Generator is not a AliGenCocktailAfterBurner, exiting\n");
       return 0x0;
    }
-   //   cout<<endl<<"Got generator"<<endl;
    return cab;
-   
  }
 /*******************************************************************/
 
@@ -974,9 +1198,9 @@ AliGenHBTprocessor* GetAliGenHBTprocessor()
  AliGenerator* g = gen->GetCurrentGenerator();
  if(g == 0x0)
   {
-    gAlice->Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
+    ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
                   "Can not get the current generator. Exiting");
-    return 0x0;
+    return 0x0;//pro forma
   }
   
  TClass* hbtpclass = AliGenHBTprocessor::Class(); //get AliGenCocktailAfterBurner TClass
@@ -984,8 +1208,8 @@ AliGenHBTprocessor* GetAliGenHBTprocessor()
  AliGenHBTprocessor* hbtp = (AliGenHBTprocessor*)gclass->DynamicCast(hbtpclass,g);//try to cast 
  if (hbtp == 0x0)
    {
-      gAlice->Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
-                    "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
+      ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
+              "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
       return 0x0;
    }
  return hbtp;
@@ -1006,19 +1230,15 @@ extern "C" void type_ofCall  alihbtp_initialize()
 
 extern "C" void type_ofCall alihbtp_getnumberevents(Int_t &nev)
  {
-  //passes number of events to the fortran 
-   if(gDebug) cout<<"alihbtp_getnumberevents("<<nev<<") ....";
-   AliGenCocktailAfterBurner* gen = GetGenerator();
-   if(!gen)
-    {
-     nev = -1;
-     return;
-    } 
-     
-    nev = gen->GetNumberOfEvents();
+   //passes number of events to the fortran 
+   if(AliGenHBTprocessor::GetDebug()) 
+    AliInfoGeneral("alihbtp_getnumberevents",Form("(%d) ....",nev));
+
+   AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
+   nev = gen->GetNumberOfEvents();
     
-   if(gDebug>5) cout<<"EXITED N Ev = "<<nev<<endl; 
-   
+   if(AliGenHBTprocessor::GetDebug()>5) 
+    AliInfoGeneral("alihbtp_getnumberevents",Form("EXITED N Ev = %d",nev));
  }
 
 /*******************************************************************/
@@ -1027,31 +1247,31 @@ extern "C" void type_ofCall  alihbtp_setactiveeventnumber(Int_t & nev)
  {
 //sets active event in generator (AliGenCocktailAfterBurner)
 
-   if(gDebug>5) cout<<"alihbtp_setactiveeventnumber("<<nev<<") ....";
-   if(gDebug>0) cout<<"Asked for event "<<nev-1<<endl;
-   AliGenCocktailAfterBurner* gen = GetGenerator();
-   if(!gen) return;
-   gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N
+   if(AliGenHBTprocessor::GetDebug()>5) 
+    ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","(%d)",nev);
+   if(AliGenHBTprocessor::GetDebug()>0)
+    ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","Asked for event %d",nev-1);
+    
+   AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
    
-   if(gDebug>5) cout<<"EXITED returned "<<nev<<endl; 
+   gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N
    
+   if(AliGenHBTprocessor::GetDebug()>5) 
+    ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","EXITED returned %d",nev);
  }
 /*******************************************************************/
  
 extern "C" void type_ofCall  alihbtp_getnumbertracks(Int_t &ntracks)
  {
 //passes number of particles in active event to the fortran  
-   if(gDebug>5) cout<<"alihbtp_getnumbertracks("<<ntracks<<") ....";
+   if(AliGenHBTprocessor::GetDebug()>5) 
+    ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","(%d)",ntracks);
 
-   AliGenCocktailAfterBurner* gen = GetGenerator();
-   if (!gen) 
-    {
-     ntracks = -1;
-     return;
-    } 
+   AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
 
-   ntracks = gen->GetActiveStack()->GetNprimary();
-   if(gDebug>5) cout<<"EXITED Ntracks = "<<ntracks<<endl; 
+   ntracks = gen->GetNumberOfTracks();
+   if(AliGenHBTprocessor::GetDebug()>5)
+    ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","EXITED Ntracks = %d",ntracks); 
  }
  
 /*******************************************************************/
@@ -1067,26 +1287,33 @@ extern "C" void type_ofCall
 // px,py,pz - momenta
 // geantpid - type of the particle - Geant Particle ID
  
-   if(gDebug>5) cout<<"alihbtp_puttrack("<<n<<") ....";
+   if(AliGenHBTprocessor::GetDebug()>5)
+    ::Info("AliGenHBTpocessor.cxx: alihbtp_puttrack","(%d)",n);
 
-   AliGenCocktailAfterBurner* gen = GetGenerator();
-   if(!gen) return;
+   AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
    
-   TParticle * track = gen->GetActiveStack()->Particle(n-1);
-   
-   AliGenHBTprocessor* g = GetAliGenHBTprocessor();
+   TParticle * track = gen->GetTrack(n-1);
+   if (track == 0x0)
+    {
+      ::Error("AliGenHBTprocessor.cxx","Can not get track from AliGenHBTprocessor");
+      return;
+    }
        
    //check to be deleted 
-   if (geantpid != (g->IdFromPDG( track->GetPdgCode() )))
+   if (geantpid != (gen->IdFromPDG( track->GetPdgCode() )))
     {
-      cerr<<endl<<" AliGenHBTprocessor.cxx: alihbtp_puttrack: SOMETHING IS GOING BAD:\n   GEANTPIDS ARE NOT THE SAME"<<endl;
+      ::Error("AliGenHBTprocessor.cxx: alihbtp_puttrack",
+              "SOMETHING IS GOING BAD:\n   GEANTPIDS ARE NOT THE SAME");
     }
    
-   if(gDebug>0)
+   if(AliGenHBTprocessor::GetDebug()>0)
      if (px != track->Px()) 
-       cout<<"Px diff. = "<<px - track->Px()<<endl;
+       ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack",
+              "Px diff. = %f", px - track->Px());
    
-   if(gDebug>3) cout<<" track->GetPdgCode() --> "<<track->GetPdgCode()<<endl;
+   if(AliGenHBTprocessor::GetDebug()>3)
+     ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack",
+            "track->GetPdgCode() --> %d",track->GetPdgCode());
    
    
    
@@ -1094,9 +1321,9 @@ extern "C" void type_ofCall
    Float_t e = TMath::Sqrt(m*m+px*px+py*py+pz*pz);
    track->SetMomentum(px,py,pz,e);
    
-   g->SetHbtPStatusCode(flag,n-1);
+   gen->SetHbtPStatusCode(flag,n-1);
    
-   if(gDebug>5) cout<<"EXITED "<<endl; 
+   if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack","EXITED ");
  }
 
 /*******************************************************************/
@@ -1112,20 +1339,10 @@ extern "C" void type_ofCall
 // px,py,pz - momenta
 // geantpid - type of the particle - Geant Particle ID
  
-   if(gDebug>5) cout<<"alihbtp_gettrack("<<n<<") ....";
-   AliGenCocktailAfterBurner* gen = GetGenerator();
-
-   if (!gen) 
-    {
-     n = -1;
-     flag =-1;
-     px = py = pz = -1;
-     geantpid = -1;
-     return;
-    } 
-
-   TParticle * track = gen->GetActiveStack()->Particle(n-1);
+   if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","(%d)",n);
    AliGenHBTprocessor* g = GetAliGenHBTprocessor();
+
+   TParticle * track = g->GetTrack(n-1);
    
    flag = g->GetHbtPStatusCode(n-1);
 
@@ -1135,14 +1352,14 @@ extern "C" void type_ofCall
   
    geantpid = g->IdFromPDG( track->GetPdgCode() );
   
-   if(gDebug>5) cout<<"EXITED "<<endl; 
+   if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","EXITED"); 
  }
 
 /*******************************************************************/
 extern "C" Float_t type_ofCall hbtpran(Int_t &)
 {
 //interface to the random number generator
-  return sRandom->Rndm();
+  return gAliRandom->Rndm();
 }        
 
 /*******************************************************************/