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