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