50G4bool G4eDPWAElasticDCS::gIsGridLoaded =
false;
51G4String G4eDPWAElasticDCS::gDataDirectory =
"";
53std::size_t G4eDPWAElasticDCS::gNumEnergies = 106;
54std::size_t G4eDPWAElasticDCS::gIndxEnergyLim = 35;
55std::size_t G4eDPWAElasticDCS::gNumThetas1 = 247;
56std::size_t G4eDPWAElasticDCS::gNumThetas2 = 128;
57G4double G4eDPWAElasticDCS::gLogMinEkin = 1.0;
58G4double G4eDPWAElasticDCS::gInvDelLogEkin = 1.0;
60std::vector<G4double> G4eDPWAElasticDCS::gTheEnergies(G4eDPWAElasticDCS::gNumEnergies);
61std::vector<G4double> G4eDPWAElasticDCS::gTheMus1(G4eDPWAElasticDCS::gNumThetas1);
62std::vector<G4double> G4eDPWAElasticDCS::gTheMus2(G4eDPWAElasticDCS::gNumThetas2);
63std::vector<G4double> G4eDPWAElasticDCS::gTheU1(G4eDPWAElasticDCS::gNumThetas1);
64std::vector<G4double> G4eDPWAElasticDCS::gTheU2(G4eDPWAElasticDCS::gNumThetas2);
67const G4double G4eDPWAElasticDCS::gXGL[] = {
68 1.98550718E-02, 1.01666761E-01, 2.37233795E-01, 4.08282679E-01,
69 5.91717321E-01, 7.62766205E-01, 8.98333239E-01, 9.80144928E-01
71const G4double G4eDPWAElasticDCS::gWGL[] = {
72 5.06142681E-02, 1.11190517E-01, 1.56853323E-01, 1.81341892E-01,
73 1.81341892E-01, 1.56853323E-01, 1.11190517E-01, 5.06142681E-02
81: fIsRestrictedSamplingRequired(isrestricted), fIsElectron(iselectron) {
82 fDCS.resize(gMaxZ+1,
nullptr);
83 fDCSLow.resize(gMaxZ+1,
nullptr);
84 fSamplingTables.resize(gMaxZ+1,
nullptr);
90 for (std::size_t i=0; i<fDCS.size(); ++i) {
91 if (fDCS[i])
delete fDCS[i];
93 for (std::size_t i=0; i<fDCSLow.size(); ++i) {
94 if (fDCSLow[i])
delete fDCSLow[i];
96 for (std::size_t i=0; i<fSamplingTables.size(); ++i) {
97 if (fSamplingTables[i])
delete fSamplingTables[i];
100 for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
101 if (fSCPCPerMatCuts[imc]) {
102 fSCPCPerMatCuts[imc]->fVSCPC.clear();
103 delete fSCPCPerMatCuts[imc];
106 fSCPCPerMatCuts.clear();
113 if (!gIsGridLoaded) {
117 BuildSmplingTableForZ(iz);
123void G4eDPWAElasticDCS::LoadGrid() {
124 G4String fname = FindDirectoryPath() +
"grid.dat";
125 std::ifstream infile(fname.c_str());
126 if (!infile.is_open()) {
128 " Problem while trying to read " + fname +
" file.\n"+
129 " G4LEDATA version should be G4EMLOW7.12 or later.\n";
130 G4Exception(
"G4eDPWAElasticDCS::ReadCompressedFile",
"em0006",
135 infile >> gNumEnergies;
136 infile >> gNumThetas1;
137 infile >> gNumThetas2;
141 gTheEnergies.resize(gNumEnergies);
142 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
144 gTheEnergies[ie] =
G4Log(dum*CLHEP::MeV);
145 if (gTheEnergies[ie]<
G4Log(2.0*CLHEP::keV)) gIndxEnergyLim = ie;
149 gLogMinEkin = gTheEnergies[0];
150 gInvDelLogEkin = (gNumEnergies-1)/(gTheEnergies[gNumEnergies-1]-gTheEnergies[0]);
152 gTheMus1.resize(gNumThetas1);
153 gTheU1.resize(gNumThetas1);
154 const double theA = 0.01;
155 for (std::size_t it=0; it<gNumThetas1; ++it) {
157 gTheMus1[it] = 0.5*(1.0-std::cos(dum*CLHEP::degree));
158 gTheU1[it] = (theA+1.0)*gTheMus1[it]/(theA+gTheMus1[it]);
161 gTheMus2.resize(gNumThetas2);
162 gTheU2.resize(gNumThetas2);
163 for (std::size_t it=0; it<gNumThetas2; ++it) {
165 gTheMus2[it] = 0.5*(1.0-std::cos(dum*CLHEP::degree));
166 gTheU2[it] = (theA+1.0)*gTheMus2[it]/(theA+gTheMus2[it]);
170 gIsGridLoaded =
true;
175void G4eDPWAElasticDCS::LoadDCSForZ(
G4int iz) {
177 if (fDCS[iz])
return;
183 const std::size_t hNumEnergries = gNumEnergies-gIndxEnergyLim;
186 for (std::size_t it=0; it<gNumThetas2; ++it) {
187 v2DHigh->
PutX(it, gTheMus2[it]);
189 for (std::size_t ie=0; ie<hNumEnergries; ++ie) {
190 v2DHigh->
PutY(ie, gTheEnergies[gIndxEnergyLim+ie]);
192 std::ostringstream ossh;
193 ossh << FindDirectoryPath() <<
"dcss/el/dcs_"<< iz<<
"_h";
194 std::istringstream finh(std::ios::in);
195 ReadCompressedFile(ossh.str(), finh);
197 for (std::size_t it=0; it<gNumThetas2; ++it) {
199 for (std::size_t ie=0; ie<hNumEnergries; ++ie) {
210 for (std::size_t it=0; it<gNumThetas1; ++it) {
211 v2DLow->
PutX(it, gTheMus1[it]);
213 for (std::size_t ie=0; ie<gIndxEnergyLim+1; ++ie) {
214 v2DLow->
PutY(ie, gTheEnergies[ie]);
216 std::ostringstream ossl;
217 ossl << FindDirectoryPath() <<
"dcss/el/dcs_"<< iz<<
"_l";
218 std::istringstream finl(std::ios::in);
219 ReadCompressedFile(ossl.str(), finl);
220 for (std::size_t it=0; it<gNumThetas1; ++it) {
222 for (std::size_t ie=0; ie<gIndxEnergyLim; ++ie) {
230 for (std::size_t it=0; it<gNumThetas1; ++it) {
231 const G4double val = v2DHigh->
Value(gTheMus1[it], gTheEnergies[gIndxEnergyLim], ix, iy);
232 v2DLow->
PutValue(it, gIndxEnergyLim, val);
235 fDCSLow[iz] = v2DLow;
241 for (std::size_t it=0; it<gNumThetas2; ++it) {
242 v2D->
PutX(it, gTheMus2[it]);
244 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
245 v2D->
PutY(ie, gTheEnergies[ie]);
247 std::ostringstream oss;
248 oss << FindDirectoryPath() <<
"dcss/pos/dcs_"<< iz;
249 std::istringstream fin(std::ios::in);
250 ReadCompressedFile(oss.str(), fin);
252 for (std::size_t it=0; it<gNumThetas2; ++it) {
254 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
276 mumin = std::max(0.0, std::min(1.0, mumin));
277 mumax = std::max(0.0, std::min(1.0, mumax));
278 if (mumin>=mumax)
return;
280 const G4double lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1],
G4Log(ekin)));
282 const G4bool isLowerGrid = (fIsElectron && lekin<gTheEnergies[gIndxEnergyLim]);
283 const std::vector<G4double>& theMuVector = (isLowerGrid) ? gTheMus1 : gTheMus2;
287 const std::size_t iMuStart = (mumin == 0.0) ? 0 : std::distance( theMuVector.begin(), std::upper_bound(theMuVector.begin(), theMuVector.end(), mumin) )-1 ;
289 const std::size_t iMuEnd = (mumax == 1.0) ? theMuVector.size()-2 : std::distance( theMuVector.begin(), std::upper_bound(theMuVector.begin(), theMuVector.end(), mumax) )-1 ;
294 for (std::size_t imu=iMuStart; imu<=iMuEnd; ++imu) {
298 const G4double low = (imu==iMuStart) ? mumin : theMuVector[imu];
299 const G4double del = (imu==iMuEnd) ? mumax-low : theMuVector[imu+1]-low;
301 for (std::size_t igl=0; igl<8; ++igl) {
302 const double mu = low + del*gXGL[igl];
303 const double dcs =
G4Exp(the2DDCS->
Value(mu, lekin, ix, iy));
304 elcsPar += gWGL[igl]*dcs;
305 tr1csPar += gWGL[igl]*dcs*mu;
306 tr2csPar += gWGL[igl]*dcs*mu*(1.0-mu);
309 tr1cs += del*tr1csPar;
310 tr2cs += del*tr2csPar;
312 elcs *= 2.0*CLHEP::twopi;
313 tr1cs *= 4.0*CLHEP::twopi;
314 tr2cs *= 12.0*CLHEP::twopi;
342 std::vector<G4double>
fW;
344 std::vector<G4double>
fA;
345 std::vector<G4double>
fB;
346 std::vector<G4int>
fI;
351void G4eDPWAElasticDCS::BuildSmplingTableForZ(
G4int iz) {
353 if (fSamplingTables[iz])
return;
356 std::vector<OneSamplingTable>* sTables =
new std::vector<OneSamplingTable>(gNumEnergies);
358 std::ostringstream oss;
359 const G4String fname = fIsElectron ?
"stables/el/" :
"stables/pos/";
360 oss << FindDirectoryPath() << fname <<
"stable_" << iz;
361 std::istringstream fin(std::ios::in);
362 ReadCompressedFile(oss.str(), fin);
363 std::size_t ndata = 0;
364 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
368 aTable.
SetSize(ndata, !fIsRestrictedSamplingRequired);
372 if (!fIsRestrictedSamplingRequired) {
373 for (std::size_t
id=0;
id<ndata; ++id) {
374 fin >> aTable.
fW[id];
376 for (std::size_t
id=0;
id<ndata; ++id) {
377 fin >> aTable.
fI[id];
380 for (std::size_t
id=0;
id<ndata; ++id) {
381 fin >> aTable.
fCum[id];
383 for (std::size_t
id=0;
id<ndata; ++id) {
384 fin >> aTable.
fA[id];
386 for (std::size_t
id=0;
id<ndata; ++id) {
387 fin >> aTable.
fB[id];
390 fSamplingTables[iz] = sTables;
405 lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], lekin));
408 const G4double rem = (lekin-gLogMinEkin)*gInvDelLogEkin;
409 const std::size_t k = (std::size_t)rem;
410 const std::size_t iekin = (r1 < rem-k) ? k+1 : k;
412 const double mu = SampleMu(iz, iekin, r2, r3);
413 return std::max(-1.0, std::min(1.0, 1.0-2.0*mu));
430 lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], lekin));
433 const G4double rem = (lekin-gLogMinEkin)*gInvDelLogEkin;
434 const std::size_t k = (size_t)rem;
435 const std::size_t iekin = (r1 < rem-k) ? k : k+1;
437 const G4double mu = SampleMu(iz, iekin, r2, 0.5*(1.0-costMax), 0.5*(1.0-costMin));
438 return std::max(-1.0, std::min(1.0, 1.0-2.0*mu));
444G4eDPWAElasticDCS::SampleMu(std::size_t izet, std::size_t ie,
G4double r1,
G4double r2) {
448 std::size_t indxl = (std::size_t)rest;
450 if (rtn.
fW[indxl] < dum0) indxl = rtn.
fI[indxl];
455 const G4double dum1 = (1.0 + rtn.
fA[indxl] + rtn.
fB[indxl]) * delta * aval;
456 const G4double dum2 = delta * delta + rtn.
fA[indxl] * delta * aval + rtn.
fB[indxl] * aval * aval;
457 const std::vector<G4double>& theUVect = (fIsElectron && ie<gIndxEnergyLim) ? gTheU1 : gTheU2;
458 const G4double u = theUVect[indxl] + dum1 / dum2 * (theUVect[indxl + 1] - theUVect[indxl]);
467 const std::vector<G4double>& uvect) {
468 const std::size_t iLow = std::distance( uvect.begin(), std::upper_bound(uvect.begin(), uvect.end(), u) )-1;
469 const G4double tau = (u-uvect[iLow])/(uvect[iLow+1]-uvect[iLow]);
470 const G4double dum0 = (1.0+stable.
fA[iLow]*(1.0-tau)+stable.
fB[iLow]);
471 const G4double dum1 = 2.0*stable.
fB[iLow]*tau;
472 const G4double dum2 = 1.0 - std::sqrt(std::max(0.0, 1.0-2.0*dum1*tau/(dum0*dum0)));
473 return std::min(stable.
fCum[iLow+1], std::max(stable.
fCum[iLow], stable.
fCum[iLow]+dum0*dum2*(stable.
fCum[iLow+1]-stable.
fCum[iLow])/dum1 ));
477G4double G4eDPWAElasticDCS::SampleMu(std::size_t izet, std::size_t ie,
G4double r1,
482 const std::vector<G4double>& theUVect = (fIsElectron && ie<gIndxEnergyLim) ? gTheU1 : gTheU2;
483 const G4double xiMin = (muMin > 0.0) ? FindCumValue((theA+1.0)*muMin/(theA+muMin), rtn, theUVect) : 0.0;
484 const G4double xiMax = (muMax < 1.0) ? FindCumValue((theA+1.0)*muMax/(theA+muMax), rtn, theUVect) : 1.0;
486 const G4double xi = xiMin+r1*(xiMax-xiMin);
487 const std::size_t iLow = std::distance( rtn.
fCum.begin(), std::upper_bound(rtn.
fCum.begin(), rtn.
fCum.end(), xi) )-1;
491 const G4double dum1 = (1.0 + rtn.
fA[iLow] + rtn.
fB[iLow]) * delta * aval;
492 const G4double dum2 = delta * delta + rtn.
fA[iLow] * delta * aval + rtn.
fB[iLow] * aval * aval;
493 const G4double u = theUVect[iLow] + dum1 / dum2 * (theUVect[iLow + 1] - theUVect[iLow]);
494 return theA*u/(theA+1.0-u);
500const G4String& G4eDPWAElasticDCS::FindDirectoryPath() {
502 if (gDataDirectory.empty()) {
503 const char* path = std::getenv(
"G4LEDATA");
505 std::ostringstream ost;
506 ost << path <<
"/dpwa/";
507 gDataDirectory = ost.str();
509 G4Exception(
"G4eDPWAElasticDCS::FindDirectoryPath()",
"em0006",
511 "Environment variable G4LEDATA not defined");
514 return gDataDirectory;
521G4eDPWAElasticDCS::ReadCompressedFile(
G4String fname, std::istringstream &iss) {
525 std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
528 G4int fileSize = in.tellg();
530 in.seekg(0,std::ios::beg);
532 Bytef *compdata =
new Bytef[fileSize];
534 in.read((
char*)compdata, fileSize);
537 uLongf complen = (uLongf)(fileSize*4);
538 Bytef *uncompdata =
new Bytef[complen];
539 while (
Z_OK!=
uncompress(uncompdata, &complen, compdata, fileSize)) {
543 uncompdata =
new Bytef[complen];
548 dataString =
new G4String((
char*)uncompdata, (
long)complen);
550 delete [] uncompdata;
553 " Problem while trying to read " + fname +
" data file.\n"+
554 " G4LEDATA version should be G4EMLOW7.12 or later.\n";
555 G4Exception(
"G4eDPWAElasticDCS::ReadCompressedFile",
"em0006",
561 iss.str(*dataString);
575 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<=fSCPCPerMatCuts[imc]->fPrCut) {
580 G4double remaining = (lekin-fSCPCPerMatCuts[imc]->fLEmin)*fSCPCPerMatCuts[imc]->fILDel;
583 G4int imax = fSCPCPerMatCuts[imc]->fVSCPC.size()-1;
585 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[imax];
587 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[lindx] + remaining*(fSCPCPerMatCuts[imc]->fVSCPC[lindx+1]-fSCPCPerMatCuts[imc]->fVSCPC[lindx]);
599 for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
600 if (fSCPCPerMatCuts[imc]) {
601 fSCPCPerMatCuts[imc]->fVSCPC.clear();
602 delete fSCPCPerMatCuts[imc];
603 fSCPCPerMatCuts[imc] =
nullptr;
608 fSCPCPerMatCuts.resize(numMatCuts,
nullptr);
610 for (std::size_t imc=0; imc<numMatCuts; ++imc) {
615 const G4double limit = fIsElectron ? 2.0*ecut : ecut;
616 const G4double min = std::max(limit,lowEnergyLimit);
617 const G4double max = highEnergyLimit;
619 fSCPCPerMatCuts[imc] =
new SCPCorrection();
620 fSCPCPerMatCuts[imc]->fIsUse =
false;
621 fSCPCPerMatCuts[imc]->fPrCut = min;
624 G4int numEbins = fNumSPCEbinPerDec*
G4lrint(std::log10(max/min));
625 numEbins = std::max(numEbins,3);
628 fSCPCPerMatCuts[imc] =
new SCPCorrection();
629 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbins,1.0);
630 fSCPCPerMatCuts[imc]->fIsUse =
true;
631 fSCPCPerMatCuts[imc]->fPrCut = min;
632 fSCPCPerMatCuts[imc]->fLEmin = lmin;
633 fSCPCPerMatCuts[imc]->fILDel = 1./ldel;
637 ComputeMParams(mat, moliereBc, moliereXc2);
639 for (
G4int ie=0; ie<numEbins; ++ie) {
644 const G4double tau = ekin/CLHEP::electron_mass_c2;
645 const G4double tauCut = ecut/CLHEP::electron_mass_c2;
647 const G4double A = moliereXc2/(4.0*tau*(tau+2.)*moliereBc);
649 const G4double dum0 = (tau+2.)/(tau+1.);
651 G4double gm =
G4Log(0.5*tau/tauCut) + (1.+dum0*dum0)*
G4Log(2.*(tau-tauCut+2.)/(tau+4.))
652 - 0.25*(tau+2.)*( tau+2.+2.*(2.*tau+1.)/(dum1*dum1))*
653 G4Log((tau+4.)*(tau-tauCut)/tau/(tau-tauCut+2.))
654 + 0.5*(tau-2*tauCut)*(tau+2.)*(1./(tau-tauCut)-1./(dum1*dum1));
661 scpCorr = 1.-gm*z0/(z0*(z0+1.));
663 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCorr;
674 const G4double finstrc2 = 5.325135453E-5;
687 for(
G4int ielem = 0; ielem < numelems; ielem++) {
688 const G4double zet = (*theElemVect)[ielem]->GetZ();
689 const G4double iwa = (*theElemVect)[ielem]->GetN();
690 const G4double ipz = theNbAtomsPerVolVect[ielem]/theTotNbAtomsPerVol;
691 const G4double dum = ipz*zet*(zet+1.0);
693 ze += dum*(-2.0/3.0)*
G4Log(zet);
694 zx += dum*
G4Log(1.0+3.34*finstrc2*zet*zet);
699 theBc = const1*density*zs/sa*
G4Exp(ze/zs)/
G4Exp(zx/zs);
700 theXc2 = const2*density*zs/sa;
702 theBc *= 1.0/CLHEP::cm;
703 theXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
double A(double temperature)
std::vector< G4Element * > G4ElementVector
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
G4double G4Log(G4double x)
G4double GetZeffective() const
const G4Material * GetMaterial() const
G4double GetDensity() const
const G4ElementVector * GetElementVector() const
G4double GetTotNbOfAtomsPerVolume() const
G4IonisParamMat * GetIonisation() const
size_t GetNumberOfElements() const
const G4double * GetVecNbOfAtomsPerVolume() const
void PutY(std::size_t idy, G4double value)
G4double Value(G4double x, G4double y, std::size_t &lastidx, std::size_t &lastidy) const
void SetBicubicInterpolation(G4bool)
void PutValue(std::size_t idx, std::size_t idy, G4double value)
void PutX(std::size_t idx, G4double value)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
std::size_t GetTableSize() const
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
static G4ProductionCutsTable * GetProductionCutsTable()
void ComputeCSPerAtom(G4int iz, G4double ekin, G4double &elcs, G4double &tr1cs, G4double &tr2cs, G4double mumin=0.0, G4double mumax=1.0)
void InitialiseForZ(std::size_t iz)
void InitSCPCorrection(G4double lowEnergyLimit, G4double highEnergyLimit)
G4double ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut, G4double ekin)
G4double SampleCosineTheta(std::size_t iz, G4double lekin, G4double r1, G4double r2, G4double r3)
G4double SampleCosineThetaRestricted(std::size_t iz, G4double lekin, G4double r1, G4double r2, G4double costMax, G4double costMin)
G4eDPWAElasticDCS(G4bool iselectron=true, G4bool isrestricted=false)
std::vector< G4double > fB
std::vector< G4double > fA
void SetSize(std::size_t nx, G4bool useAlias)
std::vector< G4double > fW
std::vector< G4double > fCum
int ZEXPORT uncompress(Bytef *dest, uLongf *destLen, const Bytef *source, uLong sourceLen)