+//expanding localy equilibated fireball with volume hadron radiation
+//thermal part: Blast wave model, Bjorken-like parametrization
+//hyght-pt: PYTHIA + jet quenching model PYQUEN
+
/*
HYDJET++
version 1.0:
*/
-
-//expanding localy equilibated fireball with volume hadron radiation
-//thermal part: Blast wave model, Bjorken-like parametrization
-//hyght-pt: PYTHIA + jet quenching model PYQUEN
-
-
#include <TLorentzVector.h>
#include <TVector3.h>
-#include <TRandom.h>
#include <TMath.h>
-#ifndef INITIALSTATEHYDJET_INCLUDED
+#ifndef INITIALSTATEHYDJET_H
#include "InitialStateHydjet.h"
#endif
#ifndef RANDARRAYFUNCTION_INCLUDED
#include "RandArrayFunction.h"
#endif
-#ifndef HADRONDECAYER_INCLUDED
-#include "HadronDecayer.h"
-#endif
#ifndef GRANDCANONICAL_INCLUDED
#include "GrandCanonical.h"
#endif
#ifndef NAStrangePotential_h
#include "StrangePotential.h"
#endif
-#ifndef NAEquationSolver_h
-#include "EquationSolver.h"
-#endif
#ifndef PARTICLE_INCLUDED
#include "Particle.h"
#endif
#ifndef PARTICLE_PDG
#include "ParticlePDG.h"
#endif
-#ifndef UKUTILITY_INCLUDED
-#include "UKUtility.h"
-#endif
#include <iostream>
#include <fstream>
#include "HYJET_COMMONS.h"
extern "C" void hyevnt_();
extern "C" void myini_();
-//extern "C" void hyinit_();
extern HYIPARCommon HYIPAR;
extern HYFPARCommon HYFPAR;
extern HYJPARCommon HYJPAR;
using std::cout;
using std::endl;
-void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocator) {
-
- //----- high-pt part------------------------------
-
- TLorentzVector partJMom, partJPos, zeroVec;
-
- HYJPAR.ishad = fParams.fIshad;
- PYQPAR.T0 = fParams.fT0;
- PYQPAR.tau0 = fParams.fTau0;
- PYQPAR.nf = fParams.fNf;
- PYQPAR.ienglu = fParams.fIenglu;
- PYQPAR.ianglu = fParams.fIanglu;
+class ParticleAllocator;
+class TRandom3;
- //std::cout<<"in InitialStateHydjet::Initialize nhsel"<<fParams.fNhsel<<std::endl;
+// declaration of the static member fLastIndex
+Int_t Particle::fLastIndex;
+void InitialStateHydjet::Initialize(List_t &source, ParticleAllocator & allocator) {
+ // Generate initial particles from the soft and hard components
- if(fParams.fNhsel != 0) {
- //generate non-equilibrated part event
+ // Initialize the static "last index variable"
+ Particle::InitIndexing();
- //std::cout<<"in InitialStateHydjet before hyevnt"<<std::endl;
+ //----- high-pt part------------------------------
+ TLorentzVector partJMom, partJPos, zeroVec;
- hyevnt_();
-
- //std::cout<<"in InitialStateHydjet after hyevnt"<<std::endl;
-
-
+ // run a HYDJET event
+ hyevnt_();
+
+ if(fParams.fNhsel != 0) {
//get number of particles in jets
Int_t numbJetPart = HYPART.njp;
- Double_t Bgen = HYFPAR.bgen;
- Int_t Njet = HYJPAR.njet;
- Int_t Nbcol = HYFPAR.nbcol;
-
- // std::cout<<"in InitialStateHydjet::Initialize bgen "<<Bgen<<" njet "<<Njet<<" "<<" Nbcol "<<Nbcol<<std::endl;
- // std::cout<<"in InitialStateHydjet::Initialize numb jet part"<<numbJetPart<<std::endl;
-
- for(Int_t i = 0; i <numbJetPart; ++i) {
- Int_t pdg = int(HYPART.ppart[i][1]);
- 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
+ for(Int_t i = 0; i <numbJetPart; i++) {
+ Int_t pdg = Int_t(HYPART.ppart[i][1]);
Double_t px = HYPART.ppart[i][2];
Double_t py = HYPART.ppart[i][3];
Double_t pz = HYPART.ppart[i][4];
Double_t vy = HYPART.ppart[i][7];
Double_t vz = HYPART.ppart[i][8];
Double_t vt = HYPART.ppart[i][9];
- partJMom.SetXYZT(px, py, pz, e);
- partJPos.SetXYZT(vx, vy, vz, vt);
- // std::cout<<" --InitialStateHydjet pdg "<<pdg<<" px "<<px<<" py "<<py<<" pz "<<pz<<" e"<<e<<std::endl;
- // std::cout<<" vx in fm "<<vx<<" vy "<<vy<<" vz "<<vz<<"vt"<<vt<<std::endl;
ParticlePDG *partDef = fDatabase->GetPDGParticle(pdg);
Int_t type =1; //from jet
- if(partDef)allocator.AddParticle(Particle(partDef, partJPos, partJMom, 0, 0,type,-1, zeroVec, zeroVec), source);
+ if(partDef) {
+ partJMom.SetXYZT(px, py, pz, e);
+ partJPos.SetXYZT(vx, vy, vz, vt);
+ Particle *particle=new Particle(partDef, partJPos, partJMom, 0, 0, type, -1, zeroVec, zeroVec);
+ particle->SetIndex();
+ allocator.AddParticle(*particle, source);
+ delete particle;
+ }
}
- } //nhsel !=0 not only hydro!
-
- //std::cout<<"in InitialStateHydjet::Initialize 2"<<std::endl;
+ } //nhsel !=0 not only hydro!
//----------HYDRO part------------------------------------------------
TLorentzVector partPos, partMom, n1, p0;
TVector3 vec3;
- const TLorentzVector zeroVec;
//set maximal hadron energy
const Double_t eMax = 5.;
//-------------------------------------
// get impact parameter
- // Double_t impactParameter = HYFPAR.bgen;
//effective volume for central
double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]
if(fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ; //Gaussian distr.
- 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);
+ 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);
//effective volume for non-central Simpson2
- Double_t VolEffnoncent = fParams.fTau * dYl * SimpsonIntegrator2(0., 2.*TMath::Pi());
+ Double_t volEffnoncent = fParams.fTau * dYl * SimpsonIntegrator2(0., 2.*TMath::Pi());
+ fVolEff = volEffcent * HYFPAR.npart/HYFPAR.npart0;
- fVolEff = VolEffcent * HYFPAR.npart/HYFPAR.npart0;
+ Double_t coeffRB = TMath::Sqrt(volEffcent * HYFPAR.npart/HYFPAR.npart0/volEffnoncent);
+ Double_t coeffR1 = HYFPAR.npart/HYFPAR.npart0;
+ coeffR1 = TMath::Power(coeffR1, 0.333333);
- Double_t coeff_RB = TMath::Sqrt(VolEffcent * HYFPAR.npart/HYFPAR.npart0/VolEffnoncent);
- Double_t coeff_R1 = HYFPAR.npart/HYFPAR.npart0;
- coeff_R1 = TMath::Power(coeff_R1, 0.333333);
+ double veff=fVolEff;
- //std::cout<<"HYFPAR.npart"<<HYFPAR.npart<<std::endl;
-
- double Veff=fVolEff;
-
- std::cout<<"Veff "<<Veff<<std::endl;
-
//------------------------------------
//cycle on particles types
for(Int_t i = 0; i < fParams.fNPartTypes; ++i) {
- double Mparam = fParams.fPartMult[2 * i] * Veff;
- Int_t multiplicity = gRandom->Poisson(Mparam);
+ Double_t mparam = fParams.fPartMult[2 * i] * veff;
+ Int_t multiplicity = gRandom->Poisson(mparam);
const Int_t encoding = fParams.fPartEnc[i];
if(multiplicity > 0) {
}
//no charm now !
if(partDef->GetCharmQNumber()!=0 || partDef->GetCharmAQNumber()!=0){
- cout<<"charmed particle, don't use now ! "<<encoding<<endl;
continue;
}
//prepare histogram to sample hadron energy:
Double_t h = (eMax - mass) / nBins;
Double_t x = mass + 0.5 * h;
- Int_t i;
- for(i = 0; i < nBins; ++i) {
- if(x>=mu && fParams.fThFO>0)probList[i] = x * TMath::Sqrt(x * x - mass * mass) /
+ Int_t ii;
+ for(ii = 0; ii < nBins; ++ii) {
+ if(x>=mu && fParams.fThFO>0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) /
(TMath::Exp((x - mu) / (fParams.fThFO)) + d);
- if(x>=mu && fParams.fThFO<=0)probList[i] = x * TMath::Sqrt(x * x - mass * mass) /
+ if(x>=mu && fParams.fThFO<=0)probList[ii] = x * TMath::Sqrt(x * x - mass * mass) /
(TMath::Exp((x - mu) / (fParams.fT)) + d);
- if(x<mu)probList[i] = 0.;
+ if(x<mu)probList[ii] = 0.;
x += h;
}
arrayFunctDistE.PrepareTable(probList);
h = (fParams.fR) / nBins;
x = 0.5 * h;
Double_t param = (fParams.fUmax) / (fParams.fR);
- for(i = 0; i < nBins; ++i) {
- probList[i] = x * TMath::CosH(param*x);
+ for(ii = 0; ii < nBins; ++ii) {
+ probList[ii] = x * TMath::CosH(param*x);
x += h;
}
arrayFunctDistR.PrepareTable(probList);
//loop over hadrons, assign hadron coordinates and momenta
Double_t weight = 0., yy = 0., px0 = 0., py0 = 0., pz0 = 0.;
Double_t e = 0., x0 = 0., y0 = 0., z0 = 0., t0 = 0., etaF = 0.;
- Double_t r, RB, phiF;
+ Double_t r, rB, phiF;
for(Int_t j = 0; j < multiplicity; ++j) {
do {
n1.SetXYZT(0.,0.,TMath::SinH(etaF),TMath::CosH(etaF));
if(TMath::Abs(etaF)>5.)continue;
- //old
- // double RBold = fParams.fR * TMath::Sqrt(1-fParams.fEpsilon);
-
- RB = fParams.fR * coeff_RB * coeff_R1;
+ rB = fParams.fR * coeffRB * coeffR1;
Double_t rho = TMath::Sqrt(gRandom->Rndm());
Double_t phi = TMath::TwoPi() * gRandom->Rndm();
- Double_t Rx = TMath::Sqrt(1-fParams.fEpsilon)*RB;
- Double_t Ry = TMath::Sqrt(1+fParams.fEpsilon)*RB;
+ Double_t rx = TMath::Sqrt(1-fParams.fEpsilon)*rB;
+ Double_t ry = TMath::Sqrt(1+fParams.fEpsilon)*rB;
- x0 = Rx * rho * TMath::Cos(phi);
- y0 = Ry * rho * TMath::Sin(phi);
+ x0 = rx * rho * TMath::Cos(phi);
+ y0 = ry * rho * TMath::Sin(phi);
r = TMath::Sqrt(x0*x0+y0*y0);
phiF = TMath::Abs(TMath::ATan(y0/x0));
if(x0>0&&y0<0)phiF = 2.*TMath::Pi()-phiF;
//proper time with emission duration
- Double_t tau = coeff_R1 * fParams.fTau + sqrt(2.) * fParams.fSigmaTau * coeff_R1 * (gRandom->Gaus());
+ Double_t tau = coeffR1 * fParams.fTau + sqrt(2.) * fParams.fSigmaTau * coeffR1 * (gRandom->Gaus());
z0 = tau * TMath::SinH(etaF);
t0 = tau * TMath::CosH(etaF);
- Double_t rhou = fParams.fUmax * r / RB;
+ Double_t rhou = fParams.fUmax * r / rB;
- //double_t rold= r/coeff_RB;
- //Double_t rhou_old = fParams.fUmax * rold / RBold;
- //std::cout<<"rhou"<<rhou<<"rhou_old"<<rhou_old<<std::endl;
Double_t uxf = TMath::SinH(rhou)*TMath::Sqrt(1+fParams.fDelta)*TMath::Cos(phiF);
Double_t uyf = TMath::SinH(rhou)*TMath::Sqrt(1-fParams.fDelta)*TMath::Sin(phiF);
Double_t php0 = TMath::TwoPi() * gRandom->Rndm();
Double_t ctp0 = 2. * gRandom->Rndm() - 1.;
- Double_t stp0 = TMath::Sqrt(1. - ctp0 * ctp0);
+ Double_t stp0 = TMath::Sqrt((1.-ctp0)*(1.+ctp0));
e = mass + (eMax - mass) * arrayFunctDistE();
- Double_t pp0 = TMath::Sqrt(e * e - mass * mass);
+ Double_t pp0 = TMath::Sqrt((e-mass)*(e+mass));
px0 = pp0 * stp0 * TMath::Sin(php0);
py0 = pp0 * stp0 * TMath::Cos(php0);
pz0 = pp0 * ctp0;
weight = (n1 * p0) /e; // weight for rdr gammar: weight = (n1 * p0) / n1[3] / e;
} while(yy >= weight);
- // if(abs(z0)>1000 || abs(x0)>1000)std::cout<<"====== etaF==== "<<etaF<<std::endl;
partMom.SetXYZT(px0, py0, pz0, e);
partPos.SetXYZT(x0, y0, z0, t0);
partMom.Boost(vec3);
Int_t type =0; //hydro
- allocator.AddParticle(Particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec), source);
-
+
+ Particle *particle=new Particle(partDef, partPos, partMom, 0., 0, type, -1, zeroVec, zeroVec);
+ particle->SetIndex();
+ allocator.AddParticle(*particle, source);
+ delete particle;
} //nhsel==4 , no hydro part
}
}
}
- std::cout<<"in InitialStateHydjet::Initialize OUT"<<std::endl;
-
}
Bool_t InitialStateHydjet::ReadParams() {
- std::cout<<"\nWelcome to Hydjet++ hadron freezeout generator!\nInput: \n\n";
+ // Read parameters from an input file in ascii
+
Float_t par[200] = {0.};
Int_t i = 0;
std::string s(40,' ');
}
Bool_t InitialStateHydjet::MultIni() {
+ // Calculate average multiplicities, chemical potentials (if necessary),
+ // initialize pyquen
+
//check and redefine input parameters
if(fParams.fTMuType>0 && fParams.fSqrtS > 2.24) {
if(fParams.fSqrtS < 2.24){
//compute strangeness potential
if(fParams.fMuB > 0.01)
fParams.fMuS = psp->CalculateStrangePotential();
- cout << "fMuS = " << fParams.fMuS << endl;
//if user choose fYlmax larger then allowed by kinematics at the specified beam energy sqrt(s)
if(fParams.fYlmax > TMath::Log(fParams.fSqrtS/0.94)){
if(fParams.fCorrS <= 0.) {
//see F. Becattini, J. Mannien, M. Gazdzicki, Phys Rev. C73 044905 (2006)
fParams.fCorrS = 1. - 0.386* TMath::Exp(-1.23*fParams.fT/fParams.fMuB);
- std::cout<<"The phenomenological f-la F. Becattini et al. PRC73 044905 (2006) for CorrS was used." << std::endl;
- std::cout<<"Strangeness suppression parameter = "<<fParams.fCorrS << std::endl;
}
std::cout<<"The phenomenological f-la J. Cleymans et al. PRC73 034905 (2006) for Tch mu_B was used." << std::endl;
}
- std::cout<<"Used eta_max = "<<fParams.fYlmax<< std::endl;
- std::cout<<"maximal allowed eta_max TMath::Log(fParams.fSqrtS/0.94)= "<<TMath::Log(fParams.fSqrtS/0.94)<<std::endl;
-
-
-
//initialisation of high-pt part
HYJPAR.nhsel = fParams.fNhsel;
PYQPAR.ianglu = fParams.fIanglu;
- myini_();
+ myini_(); //
// calculation of multiplicities of different particle species
// according to the grand canonical approach
GrandCanonical gc(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
- GrandCanonical gc_ch(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
- GrandCanonical gc_pi_th(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
- GrandCanonical gc_th_0(15, fParams.fThFO, 0., 0., 0.);
-
- // std::ofstream outMult("densities.txt");
- // outMult<<"encoding particle density chemical potential "<<std::endl;
-
+ GrandCanonical gcCh(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
+ GrandCanonical gcPiTh(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
+ GrandCanonical gcTh0(15, fParams.fThFO, 0., 0., 0.);
//effective volume for central
double dYl= 2 * fParams.fYlmax; //uniform distr. [-Ylmax; Ylmax]
if (fParams.fEtaType >0) dYl = TMath::Sqrt(2 * TMath::Pi()) * fParams.fYlmax ; //Gaussian distr.
fVolEff = 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);
- std::cout <<"central Effective volume = " << fVolEff << " [fm^3]" << std::endl;
- Double_t particleDensity_pi_ch=0;
- Double_t particleDensity_pi_th=0;
- // Double_t particleDensity_th_0=0;
+ Double_t particleDensityPiCh=0;
+ Double_t particleDensityPiTh=0;
if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
- GrandCanonical gc_ch(15, fParams.fT, fParams.fMuB, fParams.fMuS, fParams.fMuI3);
- GrandCanonical gc_pi_th(15, fParams.fThFO, 0., 0., fParams.fMu_th_pip);
- GrandCanonical gc_th_0(15, fParams.fThFO, 0., 0., 0.);
- particleDensity_pi_ch = gc_ch.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
- particleDensity_pi_th = gc_pi_th.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
+ particleDensityPiCh = gcCh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
+ particleDensityPiTh = gcPiTh.ParticleNumberDensity(fDatabase->GetPDGParticle(211));
}
for(Int_t particleIndex = 0; particleIndex < fDatabase->GetNParticles(); particleIndex++) {
Int_t encoding = currParticle->GetPDG();
//strangeness supression
Double_t gammaS = 1;
- Int_t S = Int_t(currParticle->GetStrangeness());
- if(encoding == 333)S = 2;
- if(fParams.fCorrS < 1. && S != 0)gammaS = TMath::Power(fParams.fCorrS,-TMath::Abs(S));
+ Int_t s = Int_t(currParticle->GetStrangeness());
+ if(encoding == 333)
+ s = 2;
+ if(fParams.fCorrS < 1. && s != 0)
+ gammaS = TMath::Power(fParams.fCorrS,-TMath::Abs(s));
//average densities
- Double_t particleDensity = gammaS*gc.ParticleNumberDensity(currParticle);
+ Double_t particleDensity = gc.ParticleNumberDensity(currParticle)/gammaS;
//compute chemical potential for single f.o. mu==mu_ch
Double_t mu = fParams.fMuB * Int_t(currParticle->GetBaryonNumber()) +
//thermal f.o.
if(fParams.fThFO != fParams.fT && fParams.fThFO > 0){
- Double_t particleDensity_ch = gc_ch.ParticleNumberDensity(currParticle);
- Double_t particleDensity_th_0 = gc_th_0.ParticleNumberDensity(currParticle);
- Double_t numb_dens_bolt = particleDensity_pi_th*particleDensity_ch/particleDensity_pi_ch;
- // Double_t coeff = ((currParticle->GetSpin() + 1.) * (currParticle->GetMass()) *
- // (currParticle->GetMass()) * fParams.fThFO / hbarc / hbarc / hbarc) /
- // (2. * TMath::Pi() * TMath::Pi()) * TMath::BesselK(2,(currParticle->GetMass())/fParams.fThFO);
- //Double_t mumy = fParams.fThFO*TMath::Log(numb_dens_bolt/coeff);
- mu = fParams.fThFO*TMath::Log(numb_dens_bolt/particleDensity_th_0);
+ Double_t particleDensityCh = gcCh.ParticleNumberDensity(currParticle);
+ Double_t particleDensityTh0 = gcTh0.ParticleNumberDensity(currParticle);
+ Double_t numbDensBolt = particleDensityPiTh*particleDensityCh/particleDensityPiCh;
+ mu = fParams.fThFO*TMath::Log(numbDensBolt/particleDensityTh0);
if(abs(encoding)==211 || encoding==111)mu= fParams.fMu_th_pip;
- particleDensity = numb_dens_bolt;
+ particleDensity = numbDensBolt;
}
- if(abs(encoding)==22)particleDensity=0;
+ // set particle number density to zero for some species
+ // photons
+ if(abs(encoding)==22)
+ particleDensity=0;
+ // K0L and K0S
+ if(abs(encoding)==130 || abs(encoding)==310) {
+ particleDensity=0;
+ }
if(particleDensity > 0.) {
// outMult<<encoding<< " " <<particleDensity<< " "<<mu<<std::endl;
}
Double_t InitialStateHydjet::SimpsonIntegrator2(Double_t a, Double_t b) {
+ // Simpson integration
Int_t nsubIntervals=10000;
Double_t h = (b - a)/nsubIntervals; //0-pi, phi
Double_t s=0;
- // Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal RB ?
-
Double_t x = 0; //phi
for(Int_t j = 1; j < nsubIntervals; j++) {
x += h; // phi
Double_t e = fParams.fEpsilon;
- Double_t RsB = fParams.fR; //test: podstavit' *coefff_RB
- Double_t RB = RsB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 RB
- Double_t sr = SimpsonIntegrator(0,RB,x);
+ Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
+ Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
+ Double_t sr = SimpsonIntegrator(0,rB,x);
s += sr;
}
return s*h;
}
Double_t InitialStateHydjet::SimpsonIntegrator(Double_t a, Double_t b, Double_t phi) {
+ // Simpson integration
Int_t nsubIntervals=100;
Double_t h = (b - a)/nsubIntervals;
Double_t s = f2(phi,a + 0.5*h);
//f2=f(phi,r)
Double_t InitialStateHydjet::f2(Double_t x, Double_t y) {
- Double_t RsB = fParams.fR; //test: podstavit' *coefff_RB
- Double_t rhou = fParams.fUmax * y / RsB;
+ // formula
+ Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
+ Double_t rhou = fParams.fUmax * y / rSB;
Double_t ff = y*TMath::CosH(rhou)*
TMath::Sqrt(1+fParams.fDelta*TMath::Cos(2*x)*TMath::TanH(rhou)*TMath::TanH(rhou));
//n_mu u^mu f-la 20
Double_t InitialStateHydjet::MidpointIntegrator2(Double_t a, Double_t b) {
+ // Perform integration through the mid-point method
Int_t nsubIntervals=2000;
Int_t nsubIntervals2=1;
Double_t h = (b - a)/nsubIntervals; //0-pi , phi
- Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal RB ?
-
+ Double_t h2 = (fParams.fR)/nsubIntervals; //0-R maximal rB ?
Double_t x = a + 0.5*h;
Double_t y = 0;
-
Double_t t = f2(x,y);
-
Double_t e = fParams.fEpsilon;
-
for(Int_t j = 1; j < nsubIntervals; j++) {
x += h; // integr phi
-
- Double_t RsB = fParams.fR; //test: podstavit' *coefff_RB
- Double_t RB = RsB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 RB
-
- nsubIntervals2 = Int_t(RB / h2)+1;
+ Double_t rSB = fParams.fR; //test: podstavit' *coefff_RB
+ Double_t rB = rSB *(TMath::Sqrt(1-e*e)/TMath::Sqrt(1+e*TMath::Cos(2*x))); //f-la7 rB
+ nsubIntervals2 = Int_t(rB / h2)+1;
// integr R
y=0;
for(Int_t i = 1; i < nsubIntervals2; i++)
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff