TUHKMgen/UHKM/InitialStateHydjet.cxx

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