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