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