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