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