EVGEN/AliGenHBTprocessor.cxx

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