THbtp/AliGenHBTprocessor.cxx

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