THbtp/AliGenHBTprocessor.cxx

   1 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   2  * See cxx source for full Copyright notice                               */
   3 // Implementation of the interface for THBTprocessor
   4 // Author: Piotr Krzysztof Skowronski <Piotr.Skowronski@cern.ch>
   5 //////////////////////////////////////////////////////////////////////////////////
   6 //Wrapper class for "hbt processor" after burner
   7 //The origibal code is written in fortran by Lanny Ray
   8 //and is put in the directory $ALICE_ROOT/HBTprocessor
   9 //Detailed description is on the top of the file hbt_event_processor.f
  10 //
  11 //This class can be used ONLY with AliGenCocktailAfterBurner wrapper generator
  12 //i.e. (in Config.C)
  13 // ....
  14 // AliGenCocktailAfterBurner = gener = new AliGenCocktailAfterBurner();
  15 // gener->SetPhiRange(0, 360); //Set global parameters
  16 // gener->SetThetaRange(35., 145.);
  17 // AliGenHIJINGpara *hijing = new AliGenHIJINGpara(10000); //Add some generator
  18 // hijing->SetMomentumRange(0, 999);
  19 // gener->AddGenerator(hijing,"HIJING PARAMETRIZATION",1);
  20 //
  21 // AliGenHBTprocessor *hbtp = new AliGenHBTprocessor(); //create object
  22 // hbtp->SetRefControl(2); //set parameters
  23 // hbtp->SetSwitch1D(1);
  24 // hbtp->SetR0(6);//fm - radius
  25 // hbtp->SetLambda(0.7);//chaocity parameter
  26 // gener->AddAfterBurner(hbtp,"HBT PROCESSOR",1); //add to the main generator
  27 //
  28 // gener->Init();
  29 //
  30 //CAUTIONS:
  31 //         A)  HBT PROCESSOR NEEDS MORE THAN ONE EVENT TO WORK
  32 //             AS MORE AS IT BETTER WORKS
  33 //         B)  IT IS ABLE TO "ADD" CORRELATIONS ONLY UP TO TWO PARTICLE TYPES AT ONES
  34 //
  35 //
  36 // Artificial particle dennsity enhancment feature
  37 // HBT Processor is unable to process correctly low multiplicity particles
  38 // even if the high statistics (more events) is supplied.
  39 // For that reason was implemented artificial multiplicity enhancement
  40 // - see also comments in HBT Processor source file (HBTP/hbt_event_processor.f)
  41 //   or web page (http://alisoft.cern.ch/people/skowron/hbtprocessor/hbteventprocessor.html)
  42 //   concerning trk_accep or  SetTrackRejectionFactor method in this class
  43 // Fortran is cheated by masking several events a single one. Number is defined by fEventMerge
  44 // variable. In case it is equal to 1, there is no masking and the feature is off.
  45 // But, for example if it is 5, and we have 100 events, number of events passed to fortran is 100/5=20
  46 // When fortran asks about multiplcity of event, let sey, 1 - multiplicity of events 5 to 9 is summed up
  47 // and returned.
  48 //
  49 //
  50 //////////////////////////////////////////////////////////////////////////////////
  51
  52 // 11.11.2001 Piotr Skowronski
  53 // Setting seed (date) in RNG in the constructor
  54
  55 // 09.10.2001 Piotr Skowronski
  56 //
  57 // Theta - Eta cohernecy facilities added:
  58 //    AliGenerator::SetThetaRange method overriden
  59 //    Static methods added
  60 //    EtaToTheta
  61 //    ThetaToEta
  62 //    DegreesToRadians
  63 //    RadiansToDegrees
  64 //
  65 // Class description comments put on proper place
  66
  67 // 27.09.2001 Piotr Skowronski
  68 // removing of redefinition of defaults velues
  69 // in method's implementation.
  70 //
  71 //
  72
  73 #include "AliGenHBTprocessor.h"
  74
  75 #include <TParticle.h>
  76 #include "THBTprocessor.h"
  77
  78 #include "AliStack.h"
  79 #include "AliGenCocktailAfterBurner.h"
  80 #include "AliLog.h"
  81
  82
  83
  84 ClassImp(AliGenHBTprocessor)
  85
  86 Int_t AliGenHBTprocessor::fgDebug = 0;
  87 static TRandom* gAliRandom;//RNG to be used by the fortran code
  88
  89 AliGenCocktailAfterBurner*  GetGenerator();
  90 /*******************************************************************/
  91
  92 AliGenHBTprocessor::AliGenHBTprocessor():
  93   AliGenerator(),
  94   fHBTprocessor(0x0),
  95   fHbtPStatCodes(0x0),
  96   fEventMerge(1)
  97 {
  98   //
  99   // Standard constructor
 100   // Sets default veues of all parameters
 101
 102   SetName("AliGenHBTprocessor");
 103   SetTitle("AliGenHBTprocessor");
 104
 105   gAliRandom = fRandom;
 106   fHBTprocessor = new THBTprocessor();
 107
 108   fNPDGCodes = 0;
 109   DefineParticles();
 110
 111   SetTrackRejectionFactor();
 112   SetRefControl();
 113   SetPIDs();
 114   SetNPIDtypes();
 115   SetDeltap();
 116   SetMaxIterations();
 117   SetDelChi();
 118   SetIRand();
 119   SetLambda();
 120   SetR1d() ;
 121   SetRSide();
 122   SetROut() ;
 123   SetRLong() ;
 124   SetRPerp();
 125   SetRParallel();
 126   SetR0();
 127   SetQ0();
 128   SetSwitch1D();
 129   SetSwitch3D();
 130   SetSwitchType();
 131   SetSwitchCoherence();
 132   SetSwitchCoulomb();
 133   SetSwitchFermiBose();
 134   //SetMomentumRange();
 135   SetPtRange();
 136   SetPxRange();
 137   SetPyRange();
 138   SetPzRange();
 139   SetPhiRange();
 140   SetEtaRange();
 141   SetNPtBins();
 142   SetNPhiBins();
 143   SetNEtaBins();
 144   SetNPxBins();
 145   SetNPyBins();
 146   SetNPzBins();
 147   SetNBins1DFineMesh();
 148   SetBinSize1DFineMesh();
 149   SetNBins1DCoarseMesh();
 150   SetBinSize1DCoarseMesh();
 151   SetNBins3DFineMesh();
 152   SetBinSize3DFineMesh();
 153   SetNBins3DCoarseMesh();
 154   SetBinSize3DCoarseMesh();
 155   SetNBins3DFineProjectMesh();
 156 }
 157
 158 /*******************************************************************/
 159
 160
 161 /*******************************************************************/
 162
 163 AliGenHBTprocessor::~AliGenHBTprocessor()
 164 {
 165 //destructor
 166   CleanStatusCodes();
 167   if (fHBTprocessor) delete fHBTprocessor; //delete generator
 168
 169 }
 170
 171 /*******************************************************************/
 172
 173 void AliGenHBTprocessor::InitStatusCodes()
 174 {
 175  //creates and inits status codes array to zero
 176   AliGenCocktailAfterBurner *cab = GetGenerator();
 177
 178   if(!cab) Fatal("InitStatusCodes()","Can not find AliGenCocktailAfterBurner generator");
 179
 180   Int_t nev = cab->GetNumberOfEvents();
 181
 182   fHbtPStatCodes = new Int_t* [nev];
 183   for( Int_t i =0; i<nev;i++)
 184   {
 185     Int_t nprim = cab->GetStack(i)->GetNprimary();
 186     fHbtPStatCodes[i] = new Int_t[nprim];
 187     for (Int_t k =0 ;k<nprim;k++)
 188       fHbtPStatCodes[i][k] =0;
 189
 190   }
 191
 192 }
 193 /*******************************************************************/
 194
 195 void AliGenHBTprocessor::CleanStatusCodes()
 196 {
 197  //Cleans up status codes
 198   if (fHbtPStatCodes)
 199   {
 200     for (Int_t i =0; i<GetGenerator()->GetNumberOfEvents(); i++)
 201       delete [] fHbtPStatCodes[i];
 202     delete fHbtPStatCodes;
 203     fHbtPStatCodes = 0;
 204   }
 205
 206 }
 207 /*******************************************************************/
 208
 209 void AliGenHBTprocessor::Init()
 210   {
 211   //sets up parameters in generator
 212
 213    THBTprocessor *thbtp = fHBTprocessor;
 214
 215
 216    thbtp->SetTrackRejectionFactor(fTrackRejectionFactor);
 217    thbtp->SetRefControl(fReferenceControl);
 218
 219    if ((fPid[0] == fPid[1]) || (fPid[0] == 0) || (fPid[1] == 0))
 220     {
 221        if (fPid[0] == 0)
 222          thbtp->SetPIDs(IdFromPDG(fPid[1]) ,0);
 223        else
 224          thbtp->SetPIDs(IdFromPDG(fPid[0]) ,0);
 225        thbtp->SetNPIDtypes(1);
 226
 227        if (fSwitchType !=1)
 228           Warning("AliGenHBTprocessor::Init","\nThere is only one particle type set,\n\
 229                    and Switch_Type differnt then 1,\n which does not make sense.\n\
 230                    Setting it to 1.\n");
 231
 232        SetSwitchType(1);
 233     }
 234    else
 235     {
 236        thbtp->SetPIDs(IdFromPDG(fPid[0]) ,IdFromPDG(fPid[1]));
 237        SetNPIDtypes(2);
 238        thbtp->SetSwitchType(fSwitchType);
 239     }
 240
 241    thbtp->SetMaxIterations(fMaxit);
 242    thbtp->SetDelChi(fDelchi);
 243    thbtp->SetIRand(fIrand);
 244    thbtp->SetLambda(fLambda);
 245    thbtp->SetR1d(fR1d);
 246    thbtp->SetRSide(fRside);
 247    thbtp->SetROut(fRout);
 248    thbtp->SetRLong(fRlong);
 249    thbtp->SetRPerp(fRperp);
 250    thbtp->SetRParallel(fRparallel);
 251    thbtp->SetR0(fR0);
 252    thbtp->SetQ0(fQ0);
 253    thbtp->SetSwitch1D(fSwitch1d);
 254    thbtp->SetSwitch3D(fSwitch3d);
 255    thbtp->SetSwitchType(fSwitchType);
 256    thbtp->SetSwitchCoherence(fSwitchCoherence);
 257    thbtp->SetSwitchCoulomb(fSwitchCoulomb);
 258    thbtp->SetSwitchFermiBose(fSwitchFermiBose);
 259    thbtp->SetPtRange(fPtMin,fPtMax);
 260    thbtp->SetPxRange(fPxMin,fPxMax);
 261    thbtp->SetPyRange(fPyMin,fPyMax);
 262    thbtp->SetPzRange(fPzMin,fPzMax);
 263    thbtp->SetPhiRange(fPhiMin*180./((Float_t)TMath::Pi())+180.0, //casting is because if fPhiMin = 180.0 then
 264                       fPhiMax*180./((Float_t)TMath::Pi())+180.0);//TMath::Pi() != TMath::Pi()*fPhiMin/180.0,
 265    thbtp->SetEtaRange(fEtaMin,fEtaMax);
 266    thbtp->SetNPtBins(fNPtBins);
 267    thbtp->SetNPhiBins(fNPhiBins);
 268    thbtp->SetNEtaBins(fNEtaBins);
 269    thbtp->SetNPxBins(fNPxBins);
 270    thbtp->SetNPyBins(fNPyBins);
 271    thbtp->SetNPzBins(fNPzBins);
 272    thbtp->SetNBins1DFineMesh(fN1dFine);
 273    thbtp->SetBinSize1DFineMesh(fBinsize1dFine);
 274    thbtp->SetNBins1DCoarseMesh(fN1dCoarse);
 275    thbtp->SetBinSize1DCoarseMesh(fBinsize1dCoarse);
 276    thbtp->SetNBins3DFineMesh(fN3dFine);
 277    thbtp->SetBinSize3DFineMesh(fBinsize3dFine);
 278    thbtp->SetNBins3DCoarseMesh(fN3dCoarse);
 279    thbtp->SetBinSize3DCoarseMesh(fBinsize3dCoarse);
 280    thbtp->SetNBins3DFineProjectMesh(fN3dFineProject);
 281
 282    thbtp->SetPrintFull(fPrintFull);
 283
 284  }
 285 /*******************************************************************/
 286
 287 void AliGenHBTprocessor::Generate()
 288  {
 289  //starts processig
 290    AliGenCocktailAfterBurner* cab = GetGenerator();
 291    if (cab == 0x0)
 292     {
 293       Fatal("Generate()","AliGenHBTprocessor needs AliGenCocktailAfterBurner to be main generator");
 294     }
 295    if (cab->GetNumberOfEvents() <2)
 296     {
 297       Fatal("Generate()",
 298             "HBT Processor needs more than 1 event to work properly,\
 299              but there is only %d set", cab->GetNumberOfEvents());
 300     }
 301
 302
 303    fHBTprocessor->Initialize(); //reset all fields of common blocks
 304                                    //in case there are many HBT processors
 305                                        //run one after one (i.e. in AliCocktailAfterBurner)
 306                                    //between Init() called and Generate there might
 307    Init();                         //be different instance running - be sure that we have our settings
 308
 309    InitStatusCodes(); //Init status codes
 310
 311    fHBTprocessor->GenerateEvent(); //Generates event
 312
 313    CleanStatusCodes(); //Clean Status codes - thet are not needed anymore
 314  }
 315
 316 /*******************************************************************/
 317
 318
 319 /*******************************************************************/
 320 void AliGenHBTprocessor::GetParticles(TClonesArray * particles) const
 321  {
 322  //practically dumm
 323    if (!particles)
 324     {
 325       Error("GetParticles","Particles has to be initialized");
 326       return;
 327     }
 328    fHBTprocessor->ImportParticles(particles);
 329  }
 330
 331 /*******************************************************************/
 332
 333 Int_t AliGenHBTprocessor::GetHbtPStatusCode(Int_t part) const
 334 {
 335 //returns the status code of the given particle in the active event
 336 //see SetActiveEvent in the bottom of AliGenHBTprocessor.cxx
 337 //and in AliCocktailAfterBurner
 338  Int_t ev, idx;
 339  GetTrackEventIndex(part,ev,idx);
 340  if ( (ev<0) || (idx<0) )
 341   {
 342     Error("GetHbtPStatusCode","GetTrackEventIndex returned error");
 343     return 0;
 344   }
 345  return fHbtPStatCodes[ev][idx];
 346
 347 }
 348
 349 /*******************************************************************/
 350 void  AliGenHBTprocessor::SetHbtPStatusCode(Int_t hbtstatcode, Int_t part)
 351 {
 352  //Sets the given status code to given particle number (part) in the active event
 353  Int_t ev, idx;
 354  GetTrackEventIndex(part,ev,idx);
 355  if ( (ev<0) || (idx<0) )
 356   {
 357     Error("GetHbtPStatusCode","GetTrackEventIndex returned error");
 358     return;
 359   }
 360  else fHbtPStatCodes[ev][idx] = hbtstatcode;
 361 }
 362
 363 /*******************************************************************/
 364
 365 void AliGenHBTprocessor::SetTrackRejectionFactor(Float_t trf) //def 1.0
 366   {
 367    //Sets the Track Rejection Factor
 368    //variates in range 0.0 <-> 1.0
 369    //Describes the factor of particles rejected from the output.
 370    //Used only in case of low muliplicity particles e.g. lambdas.
 371    //Processor generates addisional particles and builds the
 372    //correletions on such a statistics.
 373    //At the end these particels are left in the event according
 374    //to this factor: 1==all particles are left
 375    //                0==all are removed
 376
 377     fTrackRejectionFactor=trf;
 378     fHBTprocessor->SetTrackRejectionFactor(trf);
 379   }
 380 /*******************************************************************/
 381
 382 void AliGenHBTprocessor::SetRefControl(Int_t rc) //default 2
 383  {
 384  //Sets the Refernce Control Switch
 385  //switch wether read reference histograms from file =1
 386  //              compute from input events =2 - default
 387    fReferenceControl=rc;
 388    fHBTprocessor->SetRefControl(rc);
 389  }
 390 /*******************************************************************/
 391
 392 void AliGenHBTprocessor::SetPIDs(Int_t pid1,Int_t pid2)
 393   {
 394    //default pi+ pi-
 395    //Sets PDG Codes of particles to be processed, default \\Pi^{+} and \\Pi^{-}
 396    //This method accepts PDG codes which is
 397    //in opposite to THBBProcessor which accepts GEANT PID
 398     if ( (pid1 == 0) && (pid2 == 0) )
 399     {
 400       Error("AliGenHBTprocessor::SetPIDs","Sensless Particle Types setting: 0 0, Ignoring\n");
 401     }
 402
 403     fPid[0]=pid1;
 404     fPid[1]=pid2;
 405
 406     if(pid1 == pid2)
 407        {
 408         fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,0);
 409         SetNPIDtypes(1);
 410         SetSwitchType(1);
 411        }
 412     else
 413        {
 414         fHBTprocessor->SetPIDs(IdFromPDG(pid1) ,IdFromPDG(pid2));
 415         SetNPIDtypes(2);
 416        }
 417   }
 418 /*******************************************************************/
 419
 420 void AliGenHBTprocessor::SetNPIDtypes(Int_t npidt)
 421   {
 422     //Number ofparticle types to be processed - default 2
 423     //see AliGenHBTprocessor::SetPIDs
 424
 425     fNPidTypes = npidt;
 426     fHBTprocessor->SetNPIDtypes(npidt);
 427   }
 428 /*******************************************************************/
 429
 430 void AliGenHBTprocessor::SetDeltap(Float_t deltp)
 431   {
 432   //default = 0.1 GeV
 433   //maximum range for random momentum shifts in GeV/c;
 434   //px,py,pz independent; Default = 0.1 GeV/c.
 435     fDeltap=deltp;
 436     fHBTprocessor->SetDeltap(deltp);
 437   }
 438 /*******************************************************************/
 439
 440 void AliGenHBTprocessor::SetMaxIterations(Int_t maxiter)
 441   {
 442   //maximum # allowed iterations thru all the
 443   //tracks for each event. Default = 50.
 444   //If maxit=0, then calculate the correlations
 445   //for the input set of events without doing the
 446   //track adjustment procedure.
 447
 448     fMaxit=maxiter;
 449     fHBTprocessor->SetMaxIterations(maxiter);
 450   }
 451
 452 /*******************************************************************/
 453 void AliGenHBTprocessor::SetDelChi(Float_t dc)
 454   {
 455     //min % change in total chi-square for which
 456     //the track adjustment iterations may stop,
 457     //Default = 0.1%.
 458
 459     fDelchi=dc;
 460     fHBTprocessor->SetDelChi(dc);
 461   }
 462 /*******************************************************************/
 463
 464 void AliGenHBTprocessor::SetIRand(Int_t irnd)
 465   {
 466     //if fact dummy - used only for compatibility
 467     //we are using AliRoot built in RNG
 468     //dummy in fact since we are using aliroot build-in RNG
 469     //Sets renaodom number generator
 470     fIrand=irnd;
 471     fHBTprocessor->SetIRand(irnd);
 472   }
 473 /*******************************************************************/
 474
 475 void AliGenHBTprocessor::SetLambda(Float_t lam)
 476   {
 477   //default = 0.6
 478   // Sets Chaoticity Parameter
 479     fLambda = lam;
 480     fHBTprocessor->SetLambda(lam);
 481   }
 482 /*******************************************************************/
 483 void AliGenHBTprocessor::SetR1d(Float_t r)
 484   {
 485     //default 7.0
 486     //Sets Spherical source model radius (fm)
 487     fR1d = r;
 488     fHBTprocessor->SetR1d(r);
 489   }
 490 /*******************************************************************/
 491 void AliGenHBTprocessor::SetRSide(Float_t rs)
 492   {
 493    //default rs = 6.0
 494    //Rside,Rout,Rlong  = Non-spherical Bertsch-Pratt source model (fm)
 495
 496     fRside = rs;
 497     fHBTprocessor->SetRSide(rs);
 498   }
 499 /*******************************************************************/
 500 void AliGenHBTprocessor::SetROut(Float_t ro)
 501   {
 502     //default ro = 7.0
 503     //Rside,Rout,Rlong  = Non-spherical Bertsch-Pratt source model (fm)
 504     fRout = ro;
 505     fHBTprocessor->SetROut(ro);
 506   }
 507 /*******************************************************************/
 508 void AliGenHBTprocessor::SetRLong(Float_t rl)
 509   {
 510     //default rl = 4.0
 511     //Rside,Rout,Rlong  = Non-spherical Bertsch-Pratt source model (fm)
 512     fRlong = rl;
 513     fHBTprocessor->SetRLong(rl);
 514   }
 515 /*******************************************************************/
 516 void AliGenHBTprocessor::SetRPerp(Float_t rp)
 517   {
 518    //default rp = 6.0
 519    //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
 520     fRperp = rp;
 521     fHBTprocessor->SetRPerp(rp);
 522   }
 523 /*******************************************************************/
 524 void AliGenHBTprocessor::SetRParallel(Float_t rprl)
 525   {
 526    //default rprl = 4.0
 527    //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
 528     fRparallel = rprl;
 529     fHBTprocessor->SetRParallel(rprl);
 530   }
 531 /*******************************************************************/
 532 void AliGenHBTprocessor::SetR0(Float_t r0)
 533   {
 534   //default r0 = 4.0
 535   //Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source model (fm).
 536     fR0 = r0;
 537     fHBTprocessor->SetR0(r0);
 538   }
 539 /*******************************************************************/
 540 void AliGenHBTprocessor::SetQ0(Float_t q0)
 541   {
 542   //default q0 = 9.0
 543   //Sets Q0                = NA35 Coulomb parameter for finite source size in (GeV/c)
 544   //                         if fSwitchCoulomb = 2
 545   //                       = Spherical Coulomb source radius in (fm)
 546   //                         if switchCoulomb = 3, used to interpolate the
 547   //                         input Pratt/Cramer discrete Coulomb source
 548   //                         radii tables.
 549     fQ0 = q0;
 550     fHBTprocessor->SetQ0(q0);
 551   }
 552
 553 /*******************************************************************/
 554 void AliGenHBTprocessor::SetSwitch1D(Int_t s1d)
 555   {
 556 //default s1d = 3
 557 // Sets fSwitch1d
 558 //                          = 0 to not compute the 1D two-body //orrelations.
 559 //                          = 1 to compute this using Q-invariant
 560 //                          = 2 to compute this using Q-total
 561 //                          = 3 to compute this using Q-3-ve//tor
 562
 563     fSwitch1d = s1d;
 564     fHBTprocessor->SetSwitch1D(s1d);
 565   }
 566 /*******************************************************************/
 567 void AliGenHBTprocessor::SetSwitch3D(Int_t s3d)
 568   {
 569 //default s3d = 0
 570 // Sets fSwitch3d
 571 //                         = 0 to not compute the 3D two-body correlations.
 572 //                         = 1 to compute this using the side-out-long form
 573 //                         = 2 to compute this using the Yanno-Koonin-Pogoredskij form.
 574
 575     fSwitch3d = s3d;
 576     fHBTprocessor->SetSwitch3D(s3d);
 577   }
 578 /*******************************************************************/
 579 void AliGenHBTprocessor::SetSwitchType(Int_t st)
 580   {
 581 //default st = 3
 582 //  Sets  switch_type       = 1 to fit only the like pair correlations
 583 //                          = 2 to fit only the unlike pair correlations
 584 //                          = 3 to fit both the like and unlike pair correl.
 585 //See SetPIDs and Init
 586 //If only one particle type is set, unly==1 makes sens
 587
 588     fSwitchType = st;
 589     fHBTprocessor->SetSwitchType(st);
 590   }
 591 /*******************************************************************/
 592 void AliGenHBTprocessor::SetSwitchCoherence(Int_t sc)
 593   {
 594 // default  sc = 0
 595 //        switchCoherence  = 0 for purely incoherent source (but can have
 596 //                              lambda < 1.0)
 597 //                          = 1 for mixed incoherent and coherent source
 598
 599     fSwitchCoherence = sc;
 600     fHBTprocessor->SetSwitchCoherence(sc);
 601   }
 602 /*******************************************************************/
 603 void AliGenHBTprocessor::SetSwitchCoulomb(Int_t scol)
 604   {
 605 //default scol = 2
 606 //        switchCoulomb    = 0 no Coulomb correction
 607 //                          = 1 Point source Gamow correction only
 608 //                          = 2 NA35 finite source size correction
 609 //                          = 3 Pratt/Cramer finite source size correction;
 610 //                              interpolated from input tables.
 611     fSwitchCoulomb =scol;
 612     fHBTprocessor->SetSwitchCoulomb(scol);
 613   }
 614 /*******************************************************************/
 615 void AliGenHBTprocessor::SetSwitchFermiBose(Int_t sfb)
 616   {
 617 //default sfb = 1
 618 //        switchFermiBose =  1 Boson pairs
 619 //                          = -1 Fermion pairs
 620
 621     fSwitchFermiBose = sfb;
 622     fHBTprocessor->SetSwitchFermiBose(sfb);
 623   }
 624 /*******************************************************************/
 625 void AliGenHBTprocessor::SetPtRange(Float_t ptmin, Float_t ptmax)
 626  {
 627 // default ptmin = 0.1, ptmax = 0.98
 628 //Sets Pt range (GeV)
 629    AliGenerator::SetPtRange(ptmin,ptmax);
 630    fHBTprocessor->SetPtRange(ptmin,ptmax);
 631  }
 632
 633 /*******************************************************************/
 634 void AliGenHBTprocessor::SetPxRange(Float_t pxmin, Float_t pxmax)
 635  {
 636 //default pxmin = -1.0, pxmax = 1.0
 637 //Sets Px range
 638   fPxMin =pxmin;
 639   fPxMax =pxmax;
 640   fHBTprocessor->SetPxRange(pxmin,pxmax);
 641  }
 642 /*******************************************************************/
 643 void AliGenHBTprocessor::SetPyRange(Float_t pymin, Float_t pymax)
 644  {
 645 //default  pymin = -1.0, pymax = 1.0
 646 //Sets Py range
 647   fPyMin =pymin;
 648   fPyMax =pymax;
 649    fHBTprocessor->SetPyRange(pymin,pymax);
 650  }
 651 /*******************************************************************/
 652 void AliGenHBTprocessor::SetPzRange(Float_t pzmin, Float_t pzmax)
 653  {
 654 //default pzmin = -3.6, pzmax = 3.6
 655 //Sets Py range
 656    fPzMin =pzmin;
 657    fPzMax =pzmax;
 658    fHBTprocessor->SetPzRange(pzmin,pzmax);
 659  }
 660 void AliGenHBTprocessor::SetMomentumRange(Float_t /*pmin*/, Float_t /*pmax*/)
 661  {
 662  //default pmin=0, pmax=0
 663  //Do not use this method!
 664     MayNotUse("AliGenHBTprocessor::SetMomentumRange Method is Dummy");
 665  }
 666
 667  /*******************************************************************/
 668 void AliGenHBTprocessor::SetPhiRange(Float_t phimin, Float_t phimax)
 669  {
 670 //
 671 //Sets \\Phi range
 672   AliGenerator::SetPhiRange(phimin,phimax);
 673
 674   fHBTprocessor->SetPhiRange(phimin+180.0,phimax+180.0);
 675  }
 676 /*******************************************************************/
 677 void AliGenHBTprocessor::SetEtaRange(Float_t etamin, Float_t etamax)
 678  {
 679 //default etamin = -1.5, etamax = 1.5
 680 //Sets \\Eta range
 681    fEtaMin= etamin;
 682    fEtaMax =etamax;
 683    fHBTprocessor->SetEtaRange(etamin,etamax);
 684
 685    //set the azimothal angle range in the AliGeneraor -
 686    //to keep coherency between azimuthal angle and pseudorapidity
 687    //DO NOT CALL this->SetThetaRange, because it calls this method (where we are)
 688    //which must cause INFINITE LOOP
 689    AliGenerator::SetThetaRange(RadiansToDegrees(EtaToTheta(fEtaMin)),
 690                                RadiansToDegrees(EtaToTheta(fEtaMax)));
 691
 692  }
 693 /*******************************************************************/
 694 void AliGenHBTprocessor::SetThetaRange(Float_t thetamin, Float_t thetamax)
 695 {
 696   //default thetamin = 0, thetamax = 180
 697   //Azimuthal angle, override AliGenerator method which uses widely (i.e. wrapper generators)
 698   //core fortran HBTProcessor uses Eta (pseudorapidity)
 699   //so these methods has to be synchronized
 700
 701   AliGenerator::SetThetaRange(thetamin,thetamax);
 702
 703   SetEtaRange( ThetaToEta(fThetaMin) , ThetaToEta(fThetaMax) );
 704
 705 }
 706
 707 /*******************************************************************/
 708 void AliGenHBTprocessor::SetNPtBins(Int_t nptbin)
 709  {
 710   //default nptbin = 50
 711   //set number of Pt bins
 712    fNPtBins= nptbin;
 713    fHBTprocessor->SetNPtBins(nptbin);
 714  }
 715 /*******************************************************************/
 716 void AliGenHBTprocessor::SetNPhiBins(Int_t nphibin)
 717 {
 718   //default nphibin = 50
 719   //set number of Phi bins
 720   fNPhiBins=nphibin;
 721   fHBTprocessor->SetNPhiBins(nphibin);
 722 }
 723 /*******************************************************************/
 724 void AliGenHBTprocessor::SetNEtaBins(Int_t netabin)
 725 {
 726   //default netabin = 50
 727   //set number of Eta bins
 728   fNEtaBins = netabin;
 729   fHBTprocessor->SetNEtaBins(netabin);
 730 }
 731 /*******************************************************************/
 732 void AliGenHBTprocessor::SetNPxBins(Int_t npxbin)
 733 {
 734   //default  npxbin = 20
 735   //set number of Px bins
 736   fNPxBins = npxbin;
 737   fHBTprocessor->SetNPxBins(npxbin);
 738 }
 739 /*******************************************************************/
 740 void AliGenHBTprocessor::SetNPyBins(Int_t npybin)
 741 {
 742   //default  npybin = 20
 743   //set number of Py bins
 744   fNPyBins = npybin;
 745   fHBTprocessor->SetNPyBins(npybin);
 746 }
 747 /*******************************************************************/
 748 void AliGenHBTprocessor::SetNPzBins(Int_t npzbin)
 749 {
 750   //default npzbin = 70
 751   //set number of Pz bins
 752   fNPzBins = npzbin;
 753   fHBTprocessor->SetNPzBins(npzbin);
 754 }
 755 /*******************************************************************/
 756 void AliGenHBTprocessor::SetNBins1DFineMesh(Int_t n)
 757 {
 758 //default n = 10
 759 //Sets the number of bins in the 1D mesh
 760    fN1dFine =n;
 761    fHBTprocessor->SetNBins1DFineMesh(n);
 762
 763 }
 764 /*******************************************************************/
 765 void AliGenHBTprocessor::SetBinSize1DFineMesh(Float_t x)
 766 {
 767 //default x=0.01
 768 //Sets the bin size in the 1D mesh
 769    fBinsize1dFine = x;
 770    fHBTprocessor->SetBinSize1DFineMesh(x);
 771 }
 772 /*******************************************************************/
 773
 774 void AliGenHBTprocessor::SetNBins1DCoarseMesh(Int_t n)
 775 {
 776 //default n =2
 777 //Sets the number of bins in the coarse 1D mesh
 778   fN1dCoarse =n;
 779   fHBTprocessor->SetNBins1DCoarseMesh(n);
 780 }
 781 /*******************************************************************/
 782 void AliGenHBTprocessor::SetBinSize1DCoarseMesh(Float_t x)
 783 {
 784 //default x=0.05
 785 //Sets the bin size in the coarse 1D mesh
 786   fBinsize1dCoarse =x;
 787   fHBTprocessor->SetBinSize1DCoarseMesh(x);
 788 }
 789 /*******************************************************************/
 790
 791 void AliGenHBTprocessor::SetNBins3DFineMesh(Int_t n)
 792 {
 793 //default n = 8
 794 //Sets the number of bins in the 3D mesh
 795   fN3dFine =n;
 796   fHBTprocessor->SetNBins3DFineMesh(n);
 797 }
 798 /*******************************************************************/
 799 void AliGenHBTprocessor::SetBinSize3DFineMesh(Float_t x)
 800 {
 801 //default x=0.01
 802 //Sets the bin size in the 3D mesh
 803   fBinsize3dFine =x;
 804   fHBTprocessor->SetBinSize3DFineMesh(x);
 805 }
 806 /*******************************************************************/
 807
 808 void AliGenHBTprocessor::SetNBins3DCoarseMesh(Int_t n)
 809 {
 810 //default n = 2
 811 //Sets the number of bins in the coarse 3D mesh
 812
 813   fN3dCoarse = n;
 814   fHBTprocessor->SetNBins3DCoarseMesh(n);
 815 }
 816 /*******************************************************************/
 817 void AliGenHBTprocessor::SetBinSize3DCoarseMesh(Float_t x)
 818 {
 819 //default x=0.08
 820 //Sets the bin size in the coarse 3D mesh
 821   fBinsize3dCoarse = x;
 822   fHBTprocessor->SetBinSize3DCoarseMesh(x);
 823 }
 824 /*******************************************************************/
 825
 826 void AliGenHBTprocessor::SetNBins3DFineProjectMesh(Int_t n )
 827 {
 828 //default n =3
 829 //Sets the number of bins in the fine project mesh
 830   fN3dFineProject = n;
 831   fHBTprocessor->SetNBins3DFineProjectMesh(n);
 832 }
 833 /*******************************************************************/
 834 void AliGenHBTprocessor::SetPrintFull(Int_t flag)
 835 {
 836 //sets the print full flag
 837  fPrintFull = flag;
 838  fHBTprocessor->SetPrintFull(flag);
 839 }
 840
 841
 842 /*******************************************************************/
 843
 844 Int_t  AliGenHBTprocessor::GetNumberOfEvents()
 845 {
 846 //returns number of available events
 847   AliGenerator* g = gAlice->Generator();
 848   AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
 849   if (cab == 0x0)
 850    {
 851      Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner");
 852      return 0;
 853    }
 854   return (Int_t)TMath::Ceil(cab->GetNumberOfEvents()/((Float_t)fEventMerge));
 855 }
 856
 857 /*******************************************************************/
 858
 859 void AliGenHBTprocessor::SetActiveEventNumber(Int_t n)
 860 {
 861 //sets the active event
 862  fActiveStack =  n*fEventMerge;
 863  AliDebug(1,Form("Settimg active event %d passed %d",fActiveStack,n));
 864 }
 865 /*******************************************************************/
 866
 867 Int_t  AliGenHBTprocessor::GetNumberOfTracks()
 868 {
 869 //returns number of tracks in active event
 870   AliGenerator* g = gAlice->Generator();
 871   AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
 872   if (cab == 0x0)
 873    {
 874      Fatal("GetNumberOfEvents","Master Generator is not an AliGenCocktailAfterBurner");
 875      return 0;
 876    }
 877  Int_t n = 0;
 878  for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++)
 879   {
 880     if (i >= GetNumberOfEvents()) break; //protection not to overshoot nb of events
 881     AliStack* stack = cab->GetStack(i);
 882     if (stack == 0x0)
 883      Error("GetNumberOfTracks","There is no stack %d",i);
 884
 885     n+=stack->GetNprimary();
 886   }
 887  return n;
 888 }
 889 /*******************************************************************/
 890
 891 void AliGenHBTprocessor::SetNEventsToMerge(Int_t nev)
 892 {
 893  //sets number of events to merge
 894  if (nev > 0 ) fEventMerge = nev;
 895 }
 896 /*******************************************************************/
 897
 898 TParticle* AliGenHBTprocessor::GetTrack(Int_t n)
 899 {
 900 //returns track that hbtp thinks is n in active event
 901   AliDebug(5,Form("n = %d",n));
 902   AliGenerator* g = gAlice->Generator();
 903   AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
 904   if (cab == 0x0)
 905    {
 906      Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner");
 907      return 0;
 908    }
 909
 910  Int_t ev, idx;
 911  GetTrackEventIndex(n,ev,idx);
 912  AliDebug(5,Form("Event = %d Particle = %d",ev,idx));
 913  if ( (ev<0) || (idx<0) )
 914   {
 915     Error("GetTrack","GetTrackEventIndex returned error");
 916     return 0x0;
 917   }
 918  AliDebug(5,Form("Number of Tracks in Event(%d) = %d",ev,cab->GetStack(ev)->GetNprimary()));
 919  return cab->GetStack(ev)->Particle(idx);//safe - in case stack does not exist
 920 }
 921 /*******************************************************************/
 922
 923 void AliGenHBTprocessor::GetTrackEventIndex(Int_t n, Int_t &evno, Int_t &index) const
 924 {
 925  //returns event(stack) number and particle index
 926   AliGenerator* g = gAlice->Generator();
 927   AliGenCocktailAfterBurner* cab = (g)?dynamic_cast<AliGenCocktailAfterBurner*>(g):0x0;
 928   if (cab == 0x0)
 929    {
 930      Fatal("GetTrackEventIndex","Master Generator is not an AliGenCocktailAfterBurner");
 931      return;
 932    }
 933
 934  evno = -1;
 935  index = -1;
 936  for (Int_t i = fActiveStack;i < fActiveStack+fEventMerge; i++)
 937   {
 938     AliStack* stack = cab->GetStack(i);
 939     if (stack == 0x0)
 940      {
 941        Error("GetTrackEventIndex","There is no stack %d",i);
 942        return;
 943      }
 944
 945     Int_t ntracks = stack->GetNprimary();
 946     AliDebug(10,Form("Event %d has %d tracks. n = %d",i,ntracks,n));
 947
 948     if ( ntracks > n)
 949      {
 950        evno = i;
 951        index = n;
 952        return ;
 953      }
 954     else
 955      {
 956        n-=ntracks;
 957        continue;
 958      }
 959   }
 960  Error("GetTrackEventIndex","Could not find given track");
 961 }
 962
 963 /*******************************************************************/
 964 /*******************************************************************/
 965 /*******************************************************************/
 966
 967
 968
 969
 970
 971 void AliGenHBTprocessor::DefineParticles()
 972 {
 973   //
 974   // Load standard numbers for GEANT particles and PDG conversion
 975   fNPDGCodes = 0; //this is done in the constructor - but in any case ...
 976
 977   fPDGCode[fNPDGCodes++]=-99;   //  0 = unused location
 978   fPDGCode[fNPDGCodes++]=22;    //  1 = photon
 979   fPDGCode[fNPDGCodes++]=-11;   //  2 = positron
 980   fPDGCode[fNPDGCodes++]=11;    //  3 = electron
 981   fPDGCode[fNPDGCodes++]=12;    //  4 = neutrino e
 982   fPDGCode[fNPDGCodes++]=-13;   //  5 = muon +
 983   fPDGCode[fNPDGCodes++]=13;    //  6 = muon -
 984   fPDGCode[fNPDGCodes++]=111;   //  7 = pi0
 985   fPDGCode[fNPDGCodes++]=211;   //  8 = pi+
 986   fPDGCode[fNPDGCodes++]=-211;  //  9 = pi-
 987   fPDGCode[fNPDGCodes++]=130;   // 10 = Kaon Long
 988   fPDGCode[fNPDGCodes++]=321;   // 11 = Kaon +
 989   fPDGCode[fNPDGCodes++]=-321;  // 12 = Kaon -
 990   fPDGCode[fNPDGCodes++]=2112;  // 13 = Neutron
 991   fPDGCode[fNPDGCodes++]=2212;  // 14 = Proton
 992   fPDGCode[fNPDGCodes++]=-2212; // 15 = Anti Proton
 993   fPDGCode[fNPDGCodes++]=310;   // 16 = Kaon Short
 994   fPDGCode[fNPDGCodes++]=221;   // 17 = Eta
 995   fPDGCode[fNPDGCodes++]=3122;  // 18 = Lambda
 996   fPDGCode[fNPDGCodes++]=3222;  // 19 = Sigma +
 997   fPDGCode[fNPDGCodes++]=3212;  // 20 = Sigma 0
 998   fPDGCode[fNPDGCodes++]=3112;  // 21 = Sigma -
 999   fPDGCode[fNPDGCodes++]=3322;  // 22 = Xi0
1000   fPDGCode[fNPDGCodes++]=3312;  // 23 = Xi-
1001   fPDGCode[fNPDGCodes++]=3334;  // 24 = Omega-
1002   fPDGCode[fNPDGCodes++]=-2112; // 25 = Anti Proton
1003   fPDGCode[fNPDGCodes++]=-3122; // 26 = Anti Proton
1004   fPDGCode[fNPDGCodes++]=-3222; // 27 = Anti Sigma -
1005   fPDGCode[fNPDGCodes++]=-3212; // 28 = Anti Sigma 0
1006   fPDGCode[fNPDGCodes++]=-3112; // 29 = Anti Sigma 0
1007   fPDGCode[fNPDGCodes++]=-3322; // 30 = Anti Xi 0
1008   fPDGCode[fNPDGCodes++]=-3312; // 31 = Anti Xi +
1009   fPDGCode[fNPDGCodes++]=-3334; // 32 = Anti Omega +
1010 }
1011
1012 /*******************************************************************/
1013 Int_t AliGenHBTprocessor::IdFromPDG(Int_t pdg) const
1014 {
1015   //
1016   // Return Geant3 code from PDG and pseudo ENDF code
1017   //
1018   for(Int_t i=0;i<fNPDGCodes;++i)
1019     if(pdg==fPDGCode[i]) return i;
1020   return -1;
1021 }
1022 Int_t AliGenHBTprocessor::PDGFromId(Int_t id) const
1023 {
1024   //
1025   // Return PDG code and pseudo ENDF code from Geant3 code
1026   //
1027   if(id>0 && id<fNPDGCodes) return fPDGCode[id];
1028   else return -1;
1029 }
1030 Double_t AliGenHBTprocessor::ThetaToEta(Double_t arg)
1031  {
1032   //converts etha(azimuthal angle) to Eta (pseudorapidity). Argument in radians
1033
1034    if(arg>= TMath::Pi()) return  708.0;//This number is the biggest wich not crashes exp(x)
1035    if(arg<= 0.0) return -708.0;//
1036    if(arg == TMath::Pi()/2.) return 0.0;//
1037
1038    if (arg > 0.0)
1039     {
1040       return TMath::Log( TMath::Tan(arg/2.)) ;
1041     }
1042    else
1043     {
1044       return -TMath::Log( TMath::Tan(-arg/2.)) ;
1045     }
1046  }
1047
1048 /*******************************************************************/
1049 /******      ROUTINES    USED    FOR     COMMUNUCATION      ********/
1050 /********************     WITH      FORTRAN     ********************/
1051 /*******************************************************************/
1052
1053 #ifndef WIN32
1054   # define hbtpran hbtpran_
1055   # define alihbtp_puttrack alihbtp_puttrack_
1056   # define alihbtp_gettrack alihbtp_gettrack_
1057   # define alihbtp_getnumberevents alihbtp_getnumberevents_
1058   # define alihbtp_getnumbertracks  alihbtp_getnumbertracks_
1059   # define alihbtp_initialize alihbtp_initialize_
1060   # define alihbtp_setactiveeventnumber alihbtp_setactiveeventnumber_
1061   # define alihbtp_setparameters alihbtp_setparameters_
1062   # define type_ofCall
1063
1064 #else
1065   # define hbtpran HBTPRAN
1066   # define alihbtp_puttrack ALIHBTP_PUTTRACK
1067   # define alihbtp_gettrack ALIHBTP_GETTRACK
1068   # define alihbtp_getnumberevents ALIHBTP_GETNUMBEREVENTS
1069   # define alihbtp_getnumbertracks  ALIHBTP_GETNUMBERTRACKS
1070   # define alihbtp_initialize ALIHBTP_INITIALIZE
1071   # define alihbtp_setactiveeventnumber ALIHBTP_SETACTIVEEVENTNUMBER
1072   # define alihbtp_setparameters ALIHBTP_SETPARAMETERS
1073   # define type_ofCall  _stdcall
1074 #endif
1075
1076 /*******************************************************************/
1077
1078 AliGenCocktailAfterBurner*  GetGenerator()
1079  {
1080    // This function has two tasks:
1081    // Check if environment is OK (exist gAlice and generator)
1082    // Returns pointer to genrator
1083    //to be changed with TFolders
1084
1085    if(!gAlice)
1086     {
1087       ::Error("AliGenHBTprocessor.cxx: GetGenerator()",
1088               "There is NO gALICE! Check what you are doing!");
1089       ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
1090               "Running HBT Processor without gAlice... Exiting \n");
1091       return 0x0;//pro forma
1092     }
1093    AliGenerator * gen = gAlice->Generator();
1094
1095    if (!gen)
1096     {
1097       ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
1098                "There is no generator in gAlice, exiting\n");
1099       return 0x0;
1100     }
1101
1102    //we do not sure actual type of the genetator
1103    //and simple casting is risky - we use ROOT machinery to do safe cast
1104
1105    TClass* cabclass = AliGenCocktailAfterBurner::Class(); //get AliGenCocktailAfterBurner TClass
1106    TClass* genclass = gen->IsA();//get TClass of the generator we got from galice
1107    //use casting implemented in TClass
1108    //cast gen to cabclass
1109    AliGenCocktailAfterBurner* cab=(AliGenCocktailAfterBurner*)genclass->DynamicCast(cabclass,gen);
1110
1111    if (cab == 0x0)//if generator that we got is not AliGenCocktailAfterBurner or its descendant we get null
1112    {              //then quit with error
1113       ::Fatal("AliGenHBTprocessor.cxx: GetGenerator()",
1114               "The main Generator is not a AliGenCocktailAfterBurner, exiting\n");
1115       return 0x0;
1116    }
1117    return cab;
1118  }
1119 /*******************************************************************/
1120
1121 AliGenHBTprocessor* GetAliGenHBTprocessor()
1122 {
1123 //returns pointer to the current instance of AliGenHBTprocessor in
1124 //AliGenCocktailAfterBurner (can be more than one)
1125 //
1126  AliGenCocktailAfterBurner* gen = GetGenerator();
1127  AliGenerator* g = gen->GetCurrentGenerator();
1128  if(g == 0x0)
1129   {
1130     ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
1131                   "Can not get the current generator. Exiting");
1132     return 0x0;//pro forma
1133   }
1134
1135  TClass* hbtpclass = AliGenHBTprocessor::Class(); //get AliGenCocktailAfterBurner TClass
1136  TClass* gclass = g->IsA();//get TClass of the current generator we got from CAB
1137  AliGenHBTprocessor* hbtp = (AliGenHBTprocessor*)gclass->DynamicCast(hbtpclass,g);//try to cast
1138  if (hbtp == 0x0)
1139    {
1140       ::Fatal("AliGenHBTprocessor.cxx: GetAliGenHBTprocessor()",
1141               "\nCurrernt generator in AliGenCocktailAfterBurner is not a AliGenHBTprocessor, exiting\n");
1142       return 0x0;
1143    }
1144  return hbtp;
1145 }
1146
1147 /*******************************************************************/
1148 extern "C" void type_ofCall alihbtp_setparameters()
1149  {
1150    //dummy
1151  }
1152
1153 extern "C" void type_ofCall  alihbtp_initialize()
1154  {
1155    //dummy
1156  }
1157
1158 /*******************************************************************/
1159
1160 extern "C" void type_ofCall alihbtp_getnumberevents(Int_t &nev)
1161  {
1162    //passes number of events to the fortran
1163    if(AliGenHBTprocessor::GetDebug())
1164     AliInfoGeneral("alihbtp_getnumberevents",Form("(%d) ....",nev));
1165
1166    AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1167    nev = gen->GetNumberOfEvents();
1168
1169    if(AliGenHBTprocessor::GetDebug()>5)
1170     AliInfoGeneral("alihbtp_getnumberevents",Form("EXITED N Ev = %d",nev));
1171  }
1172
1173 /*******************************************************************/
1174
1175 extern "C" void type_ofCall  alihbtp_setactiveeventnumber(Int_t & nev)
1176  {
1177 //sets active event in generator (AliGenCocktailAfterBurner)
1178
1179    if(AliGenHBTprocessor::GetDebug()>5)
1180     ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","(%d)",nev);
1181    if(AliGenHBTprocessor::GetDebug()>0)
1182     ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","Asked for event %d",nev-1);
1183
1184    AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1185
1186    gen->SetActiveEventNumber(nev - 1); //fortran numerates events from 1 to N
1187
1188    if(AliGenHBTprocessor::GetDebug()>5)
1189     ::Info("AliGenHBTpocessor.cxx: alihbtp_setactiveeventnumber","EXITED returned %d",nev);
1190  }
1191 /*******************************************************************/
1192
1193 extern "C" void type_ofCall  alihbtp_getnumbertracks(Int_t &ntracks)
1194  {
1195 //passes number of particles in active event to the fortran
1196    if(AliGenHBTprocessor::GetDebug()>5)
1197     ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","(%d)",ntracks);
1198
1199    AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1200
1201    ntracks = gen->GetNumberOfTracks();
1202    if(AliGenHBTprocessor::GetDebug()>5)
1203     ::Info("AliGenHBTpocessor.cxx: alihbtp_getnumbertracks","EXITED Ntracks = %d",ntracks);
1204  }
1205
1206 /*******************************************************************/
1207
1208 extern "C" void type_ofCall
1209    alihbtp_puttrack(Int_t & n,Int_t& flag, Float_t& px,
1210                     Float_t& py, Float_t& pz, Int_t& geantpid)
1211  {
1212 //sets new parameters (momenta) in track number n
1213 // in the active event
1214 // n - number of the track in active event
1215 // flag - flag of the track
1216 // px,py,pz - momenta
1217 // geantpid - type of the particle - Geant Particle ID
1218
1219    if(AliGenHBTprocessor::GetDebug()>5)
1220     ::Info("AliGenHBTpocessor.cxx: alihbtp_puttrack","(%d)",n);
1221
1222    AliGenHBTprocessor* gen = GetAliGenHBTprocessor();//we dont check because it is done in function
1223
1224    TParticle * track = gen->GetTrack(n-1);
1225    if (track == 0x0)
1226     {
1227       ::Error("AliGenHBTprocessor.cxx","Can not get track from AliGenHBTprocessor");
1228       return;
1229     }
1230
1231    //check to be deleted
1232    if (geantpid != (gen->IdFromPDG( track->GetPdgCode() )))
1233     {
1234       ::Error("AliGenHBTprocessor.cxx: alihbtp_puttrack",
1235               "SOMETHING IS GOING BAD:\n   GEANTPIDS ARE NOT THE SAME");
1236     }
1237
1238    if(AliGenHBTprocessor::GetDebug()>0)
1239      if (px != track->Px())
1240        ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack",
1241               "Px diff. = %f", px - track->Px());
1242
1243    if(AliGenHBTprocessor::GetDebug()>3)
1244      ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack",
1245             "track->GetPdgCode() --> %d",track->GetPdgCode());
1246
1247
1248
1249    Float_t m = track->GetMass();
1250    Float_t e = TMath::Sqrt(m*m+px*px+py*py+pz*pz);
1251    track->SetMomentum(px,py,pz,e);
1252
1253    gen->SetHbtPStatusCode(flag,n-1);
1254
1255    if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_puttrack","EXITED ");
1256  }
1257
1258 /*******************************************************************/
1259
1260 extern "C" void type_ofCall
1261   alihbtp_gettrack(Int_t & n,Int_t & flag, Float_t & px,
1262                    Float_t & py, Float_t & pz, Int_t & geantpid)
1263
1264  {
1265 //passes track parameters to the fortran
1266 // n - number of the track in active event
1267 // flag - flag of the track
1268 // px,py,pz - momenta
1269 // geantpid - type of the particle - Geant Particle ID
1270
1271    if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","(%d)",n);
1272    AliGenHBTprocessor* g = GetAliGenHBTprocessor();
1273
1274    TParticle * track = g->GetTrack(n-1);
1275
1276    flag = g->GetHbtPStatusCode(n-1);
1277
1278    px = (Float_t)track->Px();
1279    py = (Float_t)track->Py();
1280    pz = (Float_t)track->Pz();
1281
1282    geantpid = g->IdFromPDG( track->GetPdgCode() );
1283
1284    if(AliGenHBTprocessor::GetDebug()>5) ::Info("AliGenHBTprocessor.cxx: alihbtp_gettrack","EXITED");
1285  }
1286
1287 /*******************************************************************/
1288 extern "C" Float_t type_ofCall hbtpran(Int_t &)
1289 {
1290 //interface to the random number generator
1291   return gAliRandom->Rndm();
1292 }
1293
1294 /*******************************************************************/
1295
1296
1297 /*******************************************************************/