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