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