56 logAtomicCrossSection(0),atomicFormFactor(0),logFormFactorTable(0),
57 pMaxTable(0),samplingTable(0),fLocalTable(
false)
81 G4double logfactor1 = std::log(10.)/250.;
85 if (logenergy < logtransitionenergy)
86 logenergy += logfactor1;
88 logenergy += logfactor2;
90 }
while (logenergy < logmaxenergy);
102 if (item.second)
delete item.second;
109 if (item.second)
delete item.second;
130 if (item.second)
delete item.second;
138 if (item.second)
delete item.second;
146 if (item.second)
delete item.second;
160 G4cout <<
"Calling G4PenelopeRayleighModel::Initialise()" <<
G4endl;
171 G4cout <<
"Calling G4PenelopeRayleighModel::Initialise() [master]" <<
G4endl;
185 pMaxTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
187 samplingTable =
new std::map<const G4Material*,G4PenelopeSamplingData*>;
201 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
222 G4cout <<
"Penelope Rayleigh model v2008 is initialized " << G4endl
241 G4cout <<
"Calling G4PenelopeRayleighModel::InitialiseLocal()" <<
G4endl;
285 G4cout <<
"Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" <<
G4endl;
307 ed <<
"Unable to retrieve the cross section table for Z=" << iZ <<
G4endl;
308 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
309 G4Exception(
"G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
324 ed <<
"Unable to find Z=" << iZ <<
" in the atomic cross section table" <<
G4endl;
325 G4Exception(
"G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
331 cross =
G4Exp(logXS);
335 G4cout <<
"Rayleigh cross section at " << energy/
keV <<
" keV for Z=" << Z <<
355 std::vector<G4double> *StechiometricFactors =
new std::vector<G4double>;
356 for (
G4int i=0;i<nElements;i++)
358 G4double fraction = fractionVector[i];
359 G4double atomicWeigth = (*elementVector)[i]->GetA()/(
g/
mole);
360 StechiometricFactors->push_back(fraction/atomicWeigth);
363 G4double MaxStechiometricFactor = 0.;
364 for (
G4int i=0;i<nElements;i++)
366 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
367 MaxStechiometricFactor = (*StechiometricFactors)[i];
369 if (MaxStechiometricFactor<1
e-16)
372 ed <<
"Inconsistent data of atomic composition for " <<
374 G4Exception(
"G4PenelopeRayleighModel::BuildFormFactorTable()",
378 for (
G4int i=0;i<nElements;i++)
379 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
383 for (
G4int i=0;i<nElements;i++)
384 atomsPerMolecule += (*StechiometricFactors)[i];
395 for (
G4int i=0;i<nElements;i++)
397 G4int iZ = (
G4int) (*elementVector)[i]->GetZ();
400 ff2 += f*f*(*StechiometricFactors)[i];
407 delete StechiometricFactors;
435 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeRayleighModel" <<
G4endl;
468 pMaxTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
470 samplingTable =
new std::map<const G4Material*,G4PenelopeSamplingData*>;
481 ed <<
"Unable to find the samplingTable data for " <<
483 ed <<
"This can happen only in Unit Tests" <<
G4endl;
484 G4Exception(
"G4PenelopeRayleighModel::SampleSecondaries()",
492 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
540 G4double LastQ2inTheTable = theDataTable->
GetX(nData-1);
558 cosTheta = 1.0-2.0*xx/q2max;
565 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
569 G4double dirX = sinTheta*std::cos(phi);
570 G4double dirY = sinTheta*std::sin(phi);
576 photonDirection1.
rotateUz(photonDirection0);
591 G4cout <<
"G4PenelopeRayleighModel::ReadDataFile()" <<
G4endl;
592 G4cout <<
"Going to read Rayleigh data files for Z=" << Z <<
G4endl;
595 char* path = std::getenv(
"G4LEDATA");
598 G4String excep =
"G4LEDATA environment variable not set!";
599 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
607 std::ostringstream ost;
609 ost << path <<
"/penelope/rayleigh/pdgra" << Z <<
".p08";
611 ost << path <<
"/penelope/rayleigh/pdgra0" << Z <<
".p08";
612 std::ifstream
file(ost.str().c_str());
616 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
621 file >> readZ >> nPoints;
623 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
626 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
627 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
633 for (
size_t i=0;i<nPoints;i++)
639 theVec->
PutValue(i,std::log(ene),std::log(xs));
640 if (
file.eof() && i != (nPoints-1))
643 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
644 ed <<
"Found less than " << nPoints <<
"entries " <<
G4endl;
645 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
651 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
652 "em2044",
FatalException,
"Unable to allocate the atomic cross section table");
662 std::ostringstream ost2;
664 ost2 << path <<
"/penelope/rayleigh/pdaff" << Z <<
".p08";
666 ost2 << path <<
"/penelope/rayleigh/pdaff0" << Z <<
".p08";
667 file.open(ost2.str().c_str());
671 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
674 file >> readZ >> nPoints;
676 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
679 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
680 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
690 for (
size_t i=0;i<nPoints;i++)
699 if (
file.eof() && i != (nPoints-1))
702 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
703 ed <<
"Found less than " << nPoints <<
"entries " <<
G4endl;
704 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
710 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
712 "Unable to allocate the atomicFormFactor data table");
729 G4double logQSquared = (QSquared>1
e-10) ? std::log(QSquared) : -23.;
740 ed <<
"Unable to retrieve F squared table for " << mat->
GetName() <<
G4endl;
741 G4Exception(
"G4PenelopeRayleighModel::GetFSquared()",
745 if (logQSquared < -20)
750 else if (logQSquared > maxlogQ2)
762 G4cout <<
"Q^2 = " << QSquared <<
" (units of 1/(m_e*c); F^2 = " << f2 <<
G4endl;
774 const size_t np = 150;
786 size_t nReducedPoints = np/4;
791 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
793 "Too few points to initialize the sampling algorithm");
795 if (q2min > (q2max-1
e-10))
797 G4cout <<
"q2min= " << q2min <<
" q2max= " << q2max <<
G4endl;
798 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
800 "Too narrow grid to initialize the sampling algorithm");
813 const G4int nip = 51;
817 for (
size_t i=1;i<NUNIF-1;i++)
825 G4cout <<
"Vector x has " << x->size() <<
" points, while NUNIF = " << NUNIF <<
G4endl;
834 for (
size_t i=0;i<NUNIF-1;i++)
847 pdfi->push_back(pdfk);
853 pdfih->push_back(pdfIK);
864 G4double next = previous + cons*((*pdfi)[
k-1]+4.0*(*pdfih)[
k-1]+(*pdfi)[
k]);
865 sumi->push_back(next);
868 G4double lastIntegral = (*sumi)[sumi->size()-1];
869 area->push_back(lastIntegral);
872 for (
size_t k=0;
k<sumi->size();
k++)
873 (*sumi)[
k] *= factor;
876 if ((*pdfi)[0] < 1
e-35)
877 (*pdfi)[0] = 1
e-5*pdfmax;
878 if ((*pdfi)[pdfi->size()-1] < 1
e-35)
879 (*pdfi)[pdfi->size()-1] = 1
e-5*pdfmax;
882 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
885 G4double C_temp = 1.0+A_temp+B_temp;
894 a->push_back(A_temp);
895 b->push_back(B_temp);
896 c->push_back(C_temp);
911 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
912 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
913 if (
k == 0 ||
k == nip-1)
914 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[
k]);
916 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
921 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
937 (*x)[x->size()-1] = q2max;
944 if (x->size() != NUNIF || a->size() != NUNIF ||
945 err->size() != NUNIF || area->size() != NUNIF)
948 ed <<
"Problem in building the Table for Sampling: array dimensions do not match" <<
G4endl;
949 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
961 for (
size_t i=0;i<err->size()-2;i++)
964 if ((*err)[i] > maxError)
966 maxError = (*err)[i];
972 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
974 x->insert(x->begin()+iErrMax+1,newx);
976 area->insert(area->begin()+iErrMax+1,0.);
977 a->insert(a->begin()+iErrMax+1,0.);
978 b->insert(b->begin()+iErrMax+1,0.);
979 c->insert(c->begin()+iErrMax+1,0.);
980 err->insert(err->begin()+iErrMax+1,0.);
983 for (
size_t i=iErrMax;i<=iErrMax+1;i++)
990 G4double dxLocal = (*x)[i+1]-(*x)[i];
997 pdfi->push_back(pdfk);
1003 pdfih->push_back(pdfIK);
1010 sumi->push_back(0.);
1014 G4double next = previous + cons*((*pdfi)[
k-1]+4.0*(*pdfih)[
k-1]+(*pdfi)[
k]);
1015 sumi->push_back(next);
1017 G4double lastIntegral = (*sumi)[sumi->size()-1];
1018 (*area)[i] = lastIntegral;
1021 G4double factor = 1.0/lastIntegral;
1022 for (
size_t k=0;
k<sumi->size();
k++)
1023 (*sumi)[
k] *= factor;
1026 if ((*pdfi)[0] < 1
e-35)
1027 (*pdfi)[0] = 1
e-5*pdfmax;
1028 if ((*pdfi)[pdfi->size()-1] < 1
e-35)
1029 (*pdfi)[pdfi->size()-1] = 1
e-5*pdfmax;
1032 G4double pui = (*pdfi)[pdfi->size()-1]*factor;
1033 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
1034 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
1035 G4double C_temp = 1.0+A_temp+B_temp;
1059 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1060 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1061 if (
k == 0 ||
k == nip-1)
1062 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[
k]);
1064 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1069 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1082 }
while(x->size() < np);
1084 if (x->size() != np || a->size() != np ||
1085 err->size() != np || area->size() != np)
1087 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1089 "Problem in building the extended Table for Sampling: array dimensions do not match ");
1096 for (
size_t i=0;i<np-1;i++)
1100 for (
size_t i=0;i<np-1;i++)
1104 errMax =
std::max(errMax,(*err)[i]);
1110 for (
size_t i=0;i<np-1;i++)
1113 PAC->push_back(previous+(*area)[i]);
1115 (*PAC)[PAC->size()-1] = 1.;
1122 std::vector<size_t> *ITTL =
new std::vector<size_t>;
1123 std::vector<size_t> *ITTU =
new std::vector<size_t>;
1126 for (
size_t i=0;i<np;i++)
1134 for (
size_t i=1;i<(np-1);i++)
1138 for (
size_t j=(*ITTL)[i-1];j<np && !found;j++)
1140 if ((*PAC)[j] > ptst)
1148 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1149 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1150 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1152 if (ITTU->size() != np || ITTU->size() != np)
1154 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1156 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1164 for (
size_t i=0;i<np;i++)
1166 theTable->
AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1171 G4cout <<
"*************************************************************************" <<
1173 G4cout <<
"Sampling table for Penelope Rayleigh scattering in " << mat->
GetName() <<
G4endl;
1202 G4cout <<
"G4PenelopeRayleighModel::BuildPMaxTable" <<
G4endl;
1203 G4cout <<
"Going to instanziate the pMaxTable !" <<
G4endl;
1205 pMaxTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
1214 G4Exception(
"G4PenelopeRayleighModel::GetPMaxTable()",
1216 "SamplingTable is not properly instantiated");
1224 ed <<
"Sampling table for material " << mat->
GetName() <<
" not found";
1225 G4Exception(
"G4PenelopeRayleighModel::GetPMaxTable()",
1237 const size_t nip = 51;
1253 size_t lowerBound = 0;
1254 size_t upperBound = tablePoints-1;
1255 while (lowerBound <= upperBound)
1257 size_t midBin = (lowerBound + upperBound)/2;
1258 if( Qm2 < theTable->GetX(midBin))
1259 { upperBound = midBin-1; }
1261 { lowerBound = midBin+1; }
1271 for (
size_t k=0;
k<nip;
k++)
1275 (theTable->
GetX(upperBound+1)-Q1);
1280 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1281 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1285 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1286 ((1.0-theB*etap*etap)*ci*(theTable->
GetX(upperBound+1)-Q1));
1287 fun->push_back(theFun);
1295 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1296 sum->push_back(secondPoint);
1297 for (
size_t hh=2;
hh<nip-1;
hh++)
1300 G4double next = previous+(13.0*((*fun)[
hh-1]+(*fun)[
hh])-
1301 (*fun)[
hh+1]-(*fun)[
hh-2])*HCONS;
1302 sum->push_back(next);
1304 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1305 (*fun)[nip-3])*CONS;
1306 sum->push_back(last);
1307 thePMax = thePAC + (*sum)[sum->size()-1];
1318 thePMax = theTable->
GetPAC(0);
1322 theVec->
PutValue(ie,energy,thePMax);
1325 pMaxTable->insert(std::make_pair(mat,theVec));
1334 G4cout <<
"*****************************************************************" <<
G4endl;
1338 G4cout <<
"*****************************************************************" <<
G4endl;