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