TUHKMgen/UHKM/InitialStateHydjet.cxx

   1 //expanding localy equilibated fireball with volume hadron radiation
   2 //thermal part: Blast wave model, Bjorken-like parametrization
   3 //hyght-pt: PYTHIA + jet quenching model PYQUEN
   4
   5 /*
   6          HYDJET++
   7          version 1.0:
   8          InitialStateHydjet is the modified InitialStateBjorken
   9          The high-pt part related with PYTHIA-PYQUEN is included
  10          InitialStateBjorken (FASTMC) was used.
  11
  12
  13
  14          InitialStateBjorken
  15          version 2.0:
  16          Ludmila Malinina  malinina@lav01.sinp.msu.ru,   SINP MSU/Moscow and JINR/Dubna
  17          Ionut Arsene  i.c.arsene@fys.uio.no,            Oslo University
  18                      June 2007
  19
  20          version 1.0:
  21          Nikolai Amelin, Ludmila Malinina, Timur Pocheptsov (C) JINR/Dubna
  22          amelin@sunhe.jinr.ru, malinina@sunhe.jinr.ru, pocheptsov@sunhe.jinr.ru
  23                            November. 2, 2005
  24
  25
  26 */
  27
  28 #include <TLorentzVector.h>
  29 #include <TVector3.h>
  30 #include <TMath.h>
  31
  32 #ifndef INITIALSTATEHYDJET_H
  33 #include "InitialStateHydjet.h"
  34 #endif
  35 #ifndef RANDARRAYFUNCTION_INCLUDED
  36 #include "RandArrayFunction.h"
  37 #endif
  38 #ifndef GRANDCANONICAL_INCLUDED
  39 #include "GrandCanonical.h"
  40 #endif
  41 #ifndef NAStrangePotential_h
  42 #include "StrangePotential.h"
  43 #endif
  44 #ifndef PARTICLE_INCLUDED
  45 #include "Particle.h"
  46 #endif
  47 #ifndef PARTICLE_PDG
  48 #include "ParticlePDG.h"
  49 #endif
  50 #include <iostream>
  51 #include <fstream>
  52 #include "HYJET_COMMONS.h"
  53 extern "C" void  hyevnt_();
  54 extern "C" void  myini_();
  55 extern HYIPARCommon HYIPAR;
  56 extern HYFPARCommon HYFPAR;
  57 extern HYJPARCommon HYJPAR;
  58 extern HYPARTCommon HYPART;
  59 extern SERVICECommon SERVICE;
  60
  61 using std::cout;
  62 using std::endl;
  63
  64 class ParticleAllocator;
  65 class TRandom3;
  66
  67 // declaration of the static member fLastIndex
  68 Int_t Particle::fLastIndex;
  69
  70 void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocator) {
  71   // Generate initial particles from the soft and hard components
  72
  73   // Initialize the static "last index variable"
  74   Particle::InitIndexing();
  75
  76   //----- high-pt part------------------------------
  77   TLorentzVector partJMom, partJPos, zeroVec;
  78
  79   // run a HYDJET event
  80   hyevnt_();
  81
  82   if(fParams.fNhsel != 0) {
  83     //get number of particles in jets
  84     Int_t numbJetPart = HYPART.njp;
  85
  86     for(Int_t i = 0; i <numbJetPart; i++) {
  87       Int_t pdg = Int_t(HYPART.ppart[i][1]);
  88       Double_t px = HYPART.ppart[i][2];
  89       Double_t py = HYPART.ppart[i][3];
  90       Double_t pz = HYPART.ppart[i][4];
  91       Double_t e =  HYPART.ppart[i][5];
  92       Double_t vx = HYPART.ppart[i][6];
  93       Double_t vy = HYPART.ppart[i][7];
  94       Double_t vz = HYPART.ppart[i][8];
  95       Double_t vt = HYPART.ppart[i][9];
  96       ParticlePDG *partDef = fDatabase->GetPDGParticle(pdg);
  97       Int_t type =1;                //from jet
  98       if(partDef) {
  99         partJMom.SetXYZT(px, py, pz, e);
 100         partJPos.SetXYZT(vx, vy, vz, vt);
 101         Particle *particle=new Particle(partDef, partJPos, partJMom, 0, 0, type, -1, zeroVec, zeroVec);
 102         particle->SetIndex();
 103         allocator.AddParticle(*particle, source);
 104         delete particle;
 105       }
 106     }
 107   }       //nhsel !=0 not only hydro!
 108
 109
 110   //----------HYDRO part------------------------------------------------
 111   if(fParams.fNhsel < 3) {
 112     const Double_t  weightMax = 2*TMath::CosH(fParams.fUmax);
 113     const Int_t nBins = 100;
 114     Double_t probList[nBins];
 115     RandArrayFunction arrayFunctDistE(nBins);
 116     RandArrayFunction arrayFunctDistR(nBins);
 117
 118     TLorentzVector partPos, partMom, n1, p0;
 119     TVector3 vec3;
 120     //set maximal hadron energy
 121     const Double_t eMax = 5.;
 122     //-------------------------------------
 123     // get impact parameter
 124
 125     //effective volume for central
 126     double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]
 127     if(fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ;  //Gaussian distr.
 128     Double_t volEffcent = 2 * TMath::Pi() * fParams.fTau * dYl *
 129     (fParams.fR * fParams.fR)/TMath::Power((fParams.fUmax),2)*
 130     ((fParams.fUmax)*TMath::SinH((fParams.fUmax))-TMath::CosH((fParams.fUmax))+ 1);
 131
 132     //effective volume for non-central Simpson2
 133     Double_t volEffnoncent = fParams.fTau * dYl * SimpsonIntegrator2(0., 2.*TMath::Pi());
 134     fVolEff = volEffcent * HYFPAR.npart/HYFPAR.npart0;
 135
 136     Double_t coeffRB = TMath::Sqrt(volEffcent * HYFPAR.npart/HYFPAR.npart0/volEffnoncent);
 137     Double_t coeffR1 = HYFPAR.npart/HYFPAR.npart0;
 138     coeffR1 = TMath::Power(coeffR1, 0.333333);
 139
 140     double veff=fVolEff;
 141
 142     //------------------------------------
 143     //cycle on particles types
 144     for(Int_t i = 0; i < fParams.fNPartTypes; ++i) {
 145       Double_t mparam = fParams.fPartMult[2 * i] * veff;
 146       Int_t multiplicity = gRandom->Poisson(mparam);
 147       const Int_t encoding = fParams.fPartEnc[i];
 148
 149       if(multiplicity > 0) {
 150         ParticlePDG *partDef = fDatabase->GetPDGParticle(encoding);
 151         if(!partDef) {
 152           Error("InitialStateHydjet::Initialize", "No particle with encoding %d", encoding);
 153           continue;
 154         }
 155         //no charm now !
 156         if(partDef->GetCharmQNumber()!=0 || partDef->GetCharmAQNumber()!=0){
 157           continue;
 158         }
 159
 160         //compute chemical potential for single f.o. mu==mu_ch
 161         //compute chemical potential for thermal f.o.
 162         Double_t mu = fParams.fPartMu[2 * i];
 163
 164         //choose Bose-Einstein or Fermi-Dirac statistics
 165         const Double_t d    = !(Int_t(2*partDef->GetSpin()) & 1) ? -1 : 1;
 166         const Double_t mass = partDef->GetMass();
 167
 168         //prepare histogram to sample hadron energy:
 169         Double_t h = (eMax - mass) / nBins;
 170         Double_t x = mass + 0.5 * h;
 171         Int_t ii;
 172         for(ii = 0; ii < nBins; ++ii) {
 173           if(x>=mu && fParams.fThFO>0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) /
 174                                         (TMath::Exp((x - mu) / (fParams.fThFO)) + d);
 175           if(x>=mu && fParams.fThFO<=0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) /
 176                                          (TMath::Exp((x - mu) / (fParams.fT)) + d);
 177           if(x<mu)probList[ii] = 0.;
 178           x += h;
 179         }
 180         arrayFunctDistE.PrepareTable(probList);
 181
 182         //prepare histogram to sample hadron transverse radius:
 183         h = (fParams.fR) / nBins;
 184         x =  0.5 * h;
 185         Double_t param = (fParams.fUmax) / (fParams.fR);
 186         for(ii = 0; ii < nBins; ++ii) {
 187           probList[ii] = x * TMath::CosH(param*x);
 188           x += h;
 189         }
 190         arrayFunctDistR.PrepareTable(probList);
 191
 192         //loop over hadrons, assign hadron coordinates and momenta
 193         Double_t weight = 0., yy = 0., px0 = 0., py0 = 0., pz0 = 0.;
 194         Double_t e = 0., x0 = 0., y0 = 0., z0 = 0., t0 = 0., etaF = 0.;
 195         Double_t r, rB, phiF;
 196
 197         for(Int_t j = 0; j < multiplicity; ++j) {
 198           do {
 199             fParams.fEtaType <=0 ? etaF = fParams.fYlmax * (2. * gRandom->Rndm() - 1.)
 200               : etaF = (fParams.fYlmax) * (gRandom->Gaus());
 201             n1.SetXYZT(0.,0.,TMath::SinH(etaF),TMath::CosH(etaF));
 202             if(TMath::Abs(etaF)>5.)continue;
 203
 204             rB = fParams.fR * coeffRB * coeffR1;
 205
 206             Double_t rho = TMath::Sqrt(gRandom->Rndm());
 207             Double_t phi = TMath::TwoPi() * gRandom->Rndm();
 208             Double_t rx =  TMath::Sqrt(1-fParams.fEpsilon)*rB;
 209             Double_t ry =  TMath::Sqrt(1+fParams.fEpsilon)*rB;
 210
 211             x0 = rx * rho * TMath::Cos(phi);
 212             y0 = ry * rho * TMath::Sin(phi);
 213             r = TMath::Sqrt(x0*x0+y0*y0);
 214             phiF = TMath::Abs(TMath::ATan(y0/x0));
 215
 216             if(x0<0&&y0>0)phiF = TMath::Pi()-phiF;
 217             if(x0<0&&y0<0)phiF = TMath::Pi()+phiF;
 218             if(x0>0&&y0<0)phiF = 2.*TMath::Pi()-phiF;
 219
 220             //proper time with emission duration
 221             Double_t tau = coeffR1 * fParams.fTau +  sqrt(2.) * fParams.fSigmaTau * coeffR1 * (gRandom->Gaus());
 222             z0 = tau  * TMath::SinH(etaF);
 223             t0 = tau  * TMath::CosH(etaF);
 224
 225             Double_t rhou = fParams.fUmax * r / rB;
 226
 227
 228             Double_t uxf = TMath::SinH(rhou)*TMath::Sqrt(1+fParams.fDelta)*TMath::Cos(phiF);
 229             Double_t uyf = TMath::SinH(rhou)*TMath::Sqrt(1-fParams.fDelta)*TMath::Sin(phiF);
 230             Double_t utf = TMath::CosH(etaF) * TMath::CosH(rhou) *
 231               TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*phiF)*TMath::TanH(rhou)*TMath::TanH(rhou));
 232             Double_t uzf = TMath::SinH(etaF) * TMath::CosH(rhou) *
 233               TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*phiF)*TMath::TanH(rhou)*TMath::TanH(rhou));
 234
 235             vec3.SetXYZ(uxf / utf, uyf / utf, uzf / utf);
 236             n1.Boost(-vec3);
 237
 238             yy = weightMax * gRandom->Rndm();
 239
 240             Double_t php0 = TMath::TwoPi() * gRandom->Rndm();
 241             Double_t ctp0 = 2. * gRandom->Rndm() - 1.;
 242             Double_t stp0 = TMath::Sqrt(1. - ctp0 * ctp0);
 243             e = mass + (eMax - mass) * arrayFunctDistE();
 244             Double_t pp0 = TMath::Sqrt(e * e - mass * mass);
 245             px0 = pp0 * stp0 * TMath::Sin(php0);
 246             py0 = pp0 * stp0 * TMath::Cos(php0);
 247             pz0 = pp0 * ctp0;
 248             p0.SetXYZT(px0, py0, pz0, e);
 249
 250             //weight for rdr
 251             weight = (n1 * p0) /e;  // weight for rdr gammar: weight = (n1 * p0) / n1[3] / e;
 252           } while(yy >= weight);
 253
 254           partMom.SetXYZT(px0, py0, pz0, e);
 255           partPos.SetXYZT(x0, y0, z0, t0);
 256           partMom.Boost(vec3);
 257           Int_t type =0; //hydro
 258
 259           Particle *particle=new Particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec);
 260           particle->SetIndex();
 261           allocator.AddParticle(*particle, source);
 262           delete particle;
 263         } //nhsel==4 , no hydro part
 264       }
 265     }
 266   }
 267
 268 }
 269
 270 Bool_t InitialStateHydjet::ReadParams() {
 271   // Read parameters from an input file in ascii
 272
 273   Float_t par[200] = {0.};
 274   Int_t i = 0;
 275   std::string s(40,' ');
 276   std::ifstream input("RunInputHydjet");
 277   if (!input) {
 278     Error("Ukm::ReadParams", "Cannot open RunInputHydjet");
 279     return kFALSE;
 280   }
 281
 282   while (std::getline(input, s)) {
 283     input>>par[i];
 284     if (i < 140)
 285       std::cout<<s<<"     =  "<<par[i]<<std::endl;
 286     ++i;
 287     std::getline(input,s);
 288   }
 289
 290   std::cout<<"\nFor output use the files RunOutput.root  \n\n"<< std::endl;
 291
 292   fParams.fNevnt  = Int_t(par[0]); //number of events
 293   fParams.fSqrtS  = par[1];        //cms energy per nucleon
 294   fParams.fAw     = par[2];        // atomic number of colliding nuclei
 295   fParams.fIfb    = Int_t(par[3]);      // flag of type of centrality generation (=0 is fixed by fBfix, not 0
 296                                          //impact parameter is generated in each event between fBfmin
 297                                           //and fBmax according with Glauber model (f-la 30)
 298   fParams.fBmin = par[4];         //minimum impact parameter in units of nuclear radius RA
 299   fParams.fBmax = par[5];         //maximum impact parameter in units of nuclear radius RA
 300   fParams.fBfix = par[6];         //fix impact parameter in units of nuclear radius RA
 301
 302   fParams.fSeed = Int_t(par[7]);         //parameter to set the random nuber seed (=0 the current time is used
 303                                    //to set the random generator seed, !=0 the value fSeed is
 304                                    //used to set the random generator seed and then the state of random
 305                                    //number generator in PYTHIA MRPY(1)=fSeed
 306
 307   fParams.fT         = par[8];     //chemical freeze-out temperature in GeV
 308   fParams.fMuB       = par[9];     //baryon potential
 309   fParams.fMuS       = par[10];    //strangeness potential
 310   fParams.fMuI3      = par[11];    //isospin potential
 311   fParams.fThFO      = par[12];    //thermal freeze-out temperature T^th in GeV
 312   fParams.fMu_th_pip = par[13];    // effective chemical potential of positivly charged pions at thermal in GeV
 313
 314
 315   fParams.fTau       = par[14];     //proper time value
 316   fParams.fSigmaTau  = par[15];     //its standart deviation (emission duration)
 317   fParams.fR         = par[16];     //maximal transverse radius
 318   fParams.fYlmax     = par[17];     //maximal longitudinal rapidity
 319   fParams.fUmax      = par[18];     //maximal transverse velocity multiplaed on \gamma_r
 320   fParams.fDelta     = par[19];     //momentum asymmetry parameter
 321   fParams.fEpsilon   = par[20];     //coordinate asymmetry parameter
 322
 323   fParams.fDecay      = Int_t(par[21]);    // flag to switch on/off hadron decays<0: decays off,>=0: decays on, (default: 0)
 324   fParams.fWeakDecay  = Int_t(par[22]);    //flag to switch on/off weak hadron decays <0: decays off, >0: decays on, (default: 0)
 325
 326   fParams.fEtaType   = Int_t(par[23]);     // flag to choose rapidity distribution, if fEtaType<=0,
 327                                      //then uniform rapidity distribution in [-fYlmax,fYlmax] if fEtaType>0,
 328                                      //then Gaussian with dispertion = fYlmax
 329
 330   fParams.fTMuType   = Int_t(par[24]);     // flag to use calculated chemical freeze-out temperature,
 331                                      //baryon potential and strangeness potential as a function of fSqrtS
 332
 333   fParams.fCorrS     = par[25];     // flag and value to include strangeness supression factor
 334   fParams.fNhsel = Int_t(par[26]);         //flag to switch on/off jet and hydro-state production (0: jet
 335                                      // production off and hydro on, 1: jet production on and jet quenching
 336                                      // off and hydro on, 2: jet production on and jet quenching on and
 337                                      // hydro on, 3: jet production on and jet quenching off and hydro
 338                                      // off, 4: jet production on and jet quenching on and hydro off
 339
 340   fParams.fIshad= Int_t(par[27]);         //flag to switch on/off impact parameter dependent nuclear
 341                                     // shadowing for gluons and light sea quarks (u,d,s) (0: shadowing off,
 342                                     // 1: shadowing on for fAw=207, 197, 110, 40, default: 1
 343
 344   fParams.fPtmin = par[28];       //minimal transverse momentum transfer p_T of hard
 345                                    // parton-parton scatterings in GeV (the PYTHIA parameter ckin(3)=fPtmin)
 346
 347   //  PYQUEN energy loss model parameters:
 348
 349   fParams.fT0 = par[29];          // initial temperature (in GeV) of QGP for
 350                                    //central Pb+Pb collisions at mid-rapidity (initial temperature for other
 351                                   //centralities and atomic numbers will be calculated automatically) (allowed range is 0.2<fT0<2)
 352
 353   fParams.fTau0= par[30];        //proper QGP formation time in fm/c (0.01<fTau0<10)
 354   fParams.fNf= Int_t(par[31]);          //number of active quark flavours N_f in QGP fNf=0, 1,2 or 3
 355   fParams.fIenglu= Int_t(par[32]);      // flag to fix type of in-medium partonic energy loss
 356                                            //(0: radiative and collisional loss, 1: radiative loss only, 2:
 357                                   //collisional loss only) (default: 0);
 358   fParams.fIanglu= Int_t(par[33]);      //flag to fix type of angular distribution of in-medium emitted
 359                                   // gluons (0: small-angular, 1: wide-angular, 2:collinear) (default: 0).
 360
 361
 362   //PYTHIA parameters:
 363   Int_t jj;
 364   for (Int_t j = 0; j <25; ++j) {
 365     jj= 35+j;
 366     SERVICE.parPYTH[j]=par[jj];
 367   }
 368
 369   // Set Random Number seed
 370
 371   gRandom->SetSeed(fParams.fSeed); //Set 0 to use the current time
 372   //to send seed in PYTHIA
 373   SERVICE.iseed_fromC=gRandom->GetSeed();
 374   std::cout<<"Seed for random number generation= "<<gRandom->GetSeed()<<std::endl;
 375
 376   fParams.fNPartTypes = 0;         //counter of hadron species
 377   return kTRUE;
 378 }
 379
 380 Bool_t InitialStateHydjet::MultIni() {
 381   // Calculate average multiplicities, chemical potentials (if necessary),
 382   // initialize pyquen
 383
 384   //check and redefine input parameters
 385   if(fParams.fTMuType>0 &&  fParams.fSqrtS > 2.24) {
 386     if(fParams.fSqrtS < 2.24){
 387       Error("InitialStateHydjet::MultIni", "SqrtS<2.24 not allowed with fParams.fTMuType>0");
 388       return 0;
 389     }
 390
 391     //sqrt(s) = 2.24 ==> T_kin = 0.8 GeV
 392     //see J. Cleymans, H. Oeschler, K. Redlich,S. Wheaton, Phys Rev. C73 034905 (2006)
 393     fParams.fMuB = 1.308/(1. + fParams.fSqrtS*0.273);
 394     fParams.fT = 0.166 - 0.139*fParams.fMuB*fParams.fMuB - 0.053*fParams.fMuB*fParams.fMuB*
 395       fParams.fMuB*fParams.fMuB;
 396     fParams.fMuI3 = 0.;
 397     fParams.fMuS = 0.;
 398     //create strange potential object and set strangeness density 0
 399     NAStrangePotential* psp = new NAStrangePotential(0., fDatabase);
 400     psp->SetBaryonPotential(fParams.fMuB);
 401     psp->SetTemperature(fParams.fT);
 402     //compute strangeness potential
 403     if(fParams.fMuB > 0.01)
 404       fParams.fMuS = psp->CalculateStrangePotential();
 405
 406     //if user choose fYlmax larger then allowed by kinematics at the specified beam energy sqrt(s)
 407     if(fParams.fYlmax > TMath::Log(fParams.fSqrtS/0.94)){
 408       Error("InitialStateHydjet::MultIni", "fParams.fYlmax > TMath::Log(fParams.fSqrtS/0.94)!!! ");
 409       return 0;
 410     }
 411
 412
 413     if(fParams.fCorrS <= 0.) {
 414       //see F. Becattini, J. Mannien, M. Gazdzicki, Phys Rev. C73 044905 (2006)
 415       fParams.fCorrS = 1. - 0.386* TMath::Exp(-1.23*fParams.fT/fParams.fMuB);
 416
 417     }
 418     std::cout<<"The phenomenological f-la J. Cleymans et al. PRC73 034905 (2006) for Tch mu_B was used." << std::endl;
 419     std::cout<<"The simulation will be done with the calculated parameters:" << std::endl;
 420     std::cout<<"Baryon chemical potential = "<<fParams.fMuB<< " [GeV]" << std::endl;
 421     std::cout<<"Strangeness chemical potential = "<<fParams.fMuS<< " [GeV]" << std::endl;
 422     std::cout<<"Isospin chemical potential = "<<fParams.fMuI3<< " [GeV]" << std::endl;
 423     std::cout<<"Strangeness suppression parameter = "<<fParams.fCorrS << std::endl;
 424     std::cout<<"Eta_max = "<<fParams.fYlmax<<  std::endl;
 425     std::cout << std::endl;
 426
 427   }
 428
 429
 430   //initialisation of high-pt part
 431
 432   HYJPAR.nhsel = fParams.fNhsel;
 433   HYJPAR.ptmin = fParams.fPtmin;
 434   HYJPAR.ishad = fParams.fIshad;
 435   HYIPAR.bminh = fParams.fBmin;
 436   HYIPAR.bmaxh = fParams.fBmax;
 437   HYIPAR.AW = fParams.fAw;
 438
 439   HYPYIN.ifb = fParams.fIfb;
 440   HYPYIN.bfix = fParams.fBfix;
 441   HYPYIN.ene = fParams.fSqrtS;
 442
 443   PYQPAR.T0 = fParams.fT0;
 444   PYQPAR.tau0 = fParams.fTau0;
 445   PYQPAR.nf = fParams.fNf;
 446   PYQPAR.ienglu = fParams.fIenglu;
 447   PYQPAR.ianglu = fParams.fIanglu;
 448
 449
 450   myini_();  //
 451
 452
 453   // calculation of  multiplicities of different particle species
 454   // according to the grand canonical approach
 455   GrandCanonical gc(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
 456   GrandCanonical gcCh(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
 457   GrandCanonical gcPiTh(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
 458   GrandCanonical gcTh0(15, fParams.fThFO, 0., 0., 0.);
 459
 460   //effective volume for central
 461   double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]
 462   if (fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ;  //Gaussian distr.
 463   fVolEff = 2 * TMath::Pi() * fParams.fTau * dYl * (fParams.fR * fParams.fR)/TMath::Power((fParams.fUmax),2) *
 464     ((fParams.fUmax)*TMath::SinH((fParams.fUmax))-TMath::CosH((fParams.fUmax))+ 1);
 465
 466   Double_t particleDensityPiCh=0;
 467   Double_t particleDensityPiTh=0;
 468
 469   if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
 470     particleDensityPiCh = gcCh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
 471     particleDensityPiTh = gcPiTh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
 472   }
 473
 474   for(Int_t particleIndex = 0; particleIndex < fDatabase->GetNParticles(); particleIndex++) {
 475     ParticlePDG *currParticle = fDatabase->GetPDGParticleByIndex(particleIndex);
 476     Int_t encoding = currParticle->GetPDG();
 477     //strangeness supression
 478     Double_t gammaS = 1;
 479     Int_t s = Int_t(currParticle->GetStrangeness());
 480     if(encoding == 333)
 481       s = 2;
 482     if(fParams.fCorrS < 1. && s != 0)
 483       gammaS = TMath::Power(fParams.fCorrS,-TMath::Abs(s));
 484     //average densities
 485     Double_t particleDensity = gc.ParticleNumberDensity(currParticle)/gammaS;
 486
 487     //compute chemical potential for single f.o. mu==mu_ch
 488     Double_t mu = fParams.fMuB  * Int_t(currParticle->GetBaryonNumber()) +
 489       fParams.fMuS  * Int_t(currParticle->GetStrangeness()) +
 490       fParams.fMuI3 * Int_t(currParticle->GetElectricCharge());
 491
 492     //thermal f.o.
 493     if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
 494       Double_t particleDensityCh = gcCh.ParticleNumberDensity(currParticle);
 495       Double_t particleDensityTh0 = gcTh0.ParticleNumberDensity(currParticle);
 496       Double_t numbDensBolt = particleDensityPiTh*particleDensityCh/particleDensityPiCh;
 497       mu = fParams.fThFO*TMath::Log(numbDensBolt/particleDensityTh0);
 498       if(abs(encoding)==211 || encoding==111)mu= fParams.fMu_th_pip;
 499       particleDensity = numbDensBolt;
 500     }
 501
 502     // set particle number density to zero for some species
 503     // photons
 504     if(abs(encoding)==22)
 505       particleDensity=0;
 506     // K0L and K0S
 507     if(abs(encoding)==130 || abs(encoding)==310) {
 508       particleDensity=0;
 509     }
 510
 511     if(particleDensity > 0.) {
 512       //      outMult<<encoding<< "         " <<particleDensity<< "      "<<mu<<std::endl;
 513       fParams.fPartEnc[fParams.fNPartTypes] = encoding;
 514       fParams.fPartMult[2 * fParams.fNPartTypes] = particleDensity;
 515       fParams.fPartMu[2 * fParams.fNPartTypes] = mu;
 516       ++fParams.fNPartTypes;
 517       if(fParams.fNPartTypes > 1000)
 518         Error("in Bool_t MultIni:", "fNPartTypes is too large %d", fParams.fNPartTypes);
 519     }
 520   }
 521   return kTRUE;
 522 }
 523
 524 Double_t InitialStateHydjet::SimpsonIntegrator2(Double_t a, Double_t b) {
 525   // Simpson integration
 526   Int_t nsubIntervals=10000;
 527   Double_t h = (b - a)/nsubIntervals; //0-pi, phi
 528   Double_t s=0;
 529   Double_t x = 0; //phi
 530   for(Int_t j = 1; j < nsubIntervals; j++) {
 531     x += h; // phi
 532     Double_t e = fParams.fEpsilon;
 533     Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
 534     Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
 535     Double_t sr = SimpsonIntegrator(0,rB,x);
 536     s += sr;
 537   }
 538   return s*h;
 539
 540 }
 541
 542 Double_t InitialStateHydjet::SimpsonIntegrator(Double_t a, Double_t b, Double_t phi) {
 543   // Simpson integration
 544   Int_t nsubIntervals=100;
 545   Double_t h = (b - a)/nsubIntervals;
 546   Double_t s = f2(phi,a + 0.5*h);
 547   Double_t t = 0.5*(f2(phi,a) + f2(phi,b));
 548   Double_t x = a;
 549   Double_t y = a + 0.5*h;
 550   for(Int_t i = 1; i < nsubIntervals; i++) {
 551     x += h;
 552     y += h;
 553     s += f2(phi,y);
 554     t += f2(phi,x);
 555   }
 556   t += 2.0*s;
 557   return t*h/3.0;
 558 }
 559
 560
 561 //f2=f(phi,r)
 562 Double_t InitialStateHydjet::f2(Double_t x, Double_t y) {
 563   // formula
 564   Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
 565   Double_t rhou =  fParams.fUmax * y / rSB;
 566   Double_t ff = y*TMath::CosH(rhou)*
 567     TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*x)*TMath::TanH(rhou)*TMath::TanH(rhou));
 568 //n_mu u^mu f-la 20
 569   return ff;
 570 }
 571
 572
 573 Double_t InitialStateHydjet::MidpointIntegrator2(Double_t a, Double_t b) {
 574   // Perform integration through the mid-point method
 575   Int_t nsubIntervals=2000;
 576   Int_t nsubIntervals2=1;
 577   Double_t h = (b - a)/nsubIntervals; //0-pi , phi
 578   Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal rB ?
 579   Double_t x = a + 0.5*h;
 580   Double_t y = 0;
 581   Double_t t = f2(x,y);
 582   Double_t e = fParams.fEpsilon;
 583   for(Int_t j = 1; j < nsubIntervals; j++) {
 584     x += h; // integr  phi
 585     Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
 586     Double_t  rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
 587     nsubIntervals2 = Int_t(rB / h2)+1;
 588     // integr R
 589     y=0;
 590     for(Int_t i = 1; i < nsubIntervals2; i++)
 591       t += f2(x,(y += h2));
 592   }
 593   return t*h*h2;
 594 }
 595