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