56 :
G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
57 logAtomicCrossSection(0),atomicFormFactor(0),logFormFactorTable(0),
58 pMaxTable(0),samplingTable(0),fLocalTable(false)
82 G4double logfactor1 = std::log(10.)/250.;
86 if (logenergy < logtransitionenergy)
87 logenergy += logfactor1;
89 logenergy += logfactor2;
91 }
while (logenergy < logmaxenergy);
100 std::map <G4int,G4PhysicsFreeVector*>::iterator i;
104 if (i->second)
delete i->second;
112 if (i->second)
delete i->second;
131 std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i;
136 if (i->second)
delete i->second;
144 if (i->second)
delete i->second;
149 std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii;
153 if (ii->second)
delete ii->second;
167 G4cout <<
"Calling G4PenelopeRayleighModel::Initialise()" <<
G4endl;
178 G4cout <<
"Calling G4PenelopeRayleighModel::Initialise() [master]" <<
G4endl;
192 pMaxTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
194 samplingTable =
new std::map<const G4Material*,G4PenelopeSamplingData*>;
208 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
229 G4cout <<
"Penelope Rayleigh model v2008 is initialized " << G4endl
248 G4cout <<
"Calling G4PenelopeRayleighModel::InitialiseLocal()" <<
G4endl;
292 G4cout <<
"Calling CrossSectionPerAtom() of G4PenelopeRayleighModel" <<
G4endl;
314 ed <<
"Unable to retrieve the cross section table for Z=" << iZ <<
G4endl;
315 ed <<
"This can happen only in Unit Tests or via G4EmCalculator" <<
G4endl;
316 G4Exception(
"G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
331 ed <<
"Unable to find Z=" << iZ <<
" in the atomic cross section table" <<
G4endl;
332 G4Exception(
"G4PenelopeRayleighModel::ComputeCrossSectionPerAtom()",
338 cross = std::exp(logXS);
341 G4cout <<
"Rayleigh cross section at " << energy/
keV <<
" keV for Z=" << Z <<
359 std::vector<G4double> *StechiometricFactors =
new std::vector<G4double>;
360 for (
G4int i=0;i<nElements;i++)
362 G4double fraction = fractionVector[i];
363 G4double atomicWeigth = (*elementVector)[i]->GetA()/(
g/
mole);
364 StechiometricFactors->push_back(fraction/atomicWeigth);
367 G4double MaxStechiometricFactor = 0.;
368 for (
G4int i=0;i<nElements;i++)
370 if ((*StechiometricFactors)[i] > MaxStechiometricFactor)
371 MaxStechiometricFactor = (*StechiometricFactors)[i];
373 if (MaxStechiometricFactor<1e-16)
376 ed <<
"Inconsistent data of atomic composition for " <<
378 G4Exception(
"G4PenelopeRayleighModel::BuildFormFactorTable()",
382 for (
G4int i=0;i<nElements;i++)
383 (*StechiometricFactors)[i] /= MaxStechiometricFactor;
387 for (
G4int i=0;i<nElements;i++)
388 atomsPerMolecule += (*StechiometricFactors)[i];
399 for (
G4int i=0;i<nElements;i++)
401 G4int iZ = (
G4int) (*elementVector)[i]->GetZ();
404 ff2 += f*f*(*StechiometricFactors)[i];
411 delete StechiometricFactors;
439 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeRayleighModel" <<
G4endl;
472 pMaxTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
474 samplingTable =
new std::map<const G4Material*,G4PenelopeSamplingData*>;
485 ed <<
"Unable to find the samplingTable data for " <<
487 ed <<
"This can happen only in Unit Tests" <<
G4endl;
488 G4Exception(
"G4PenelopeRayleighModel::SampleSecondaries()",
496 G4int iZ = (
G4int) theElementVector->at(j)->GetZ();
527 G4double qmax = 2.0*photonEnergy0/electron_mass_c2;
536 G4double G = 0.5*(1+cosTheta*cosTheta);
544 G4double LastQ2inTheTable = theDataTable->
GetX(nData-1);
562 cosTheta = 1.0-2.0*xx/q2max;
563 G4double G = 0.5*(1+cosTheta*cosTheta);
569 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
573 G4double dirX = sinTheta*std::cos(phi);
574 G4double dirY = sinTheta*std::sin(phi);
580 photonDirection1.rotateUz(photonDirection0);
595 G4cout <<
"G4PenelopeRayleighModel::ReadDataFile()" <<
G4endl;
596 G4cout <<
"Going to read Rayleigh data files for Z=" << Z <<
G4endl;
599 char* path = getenv(
"G4LEDATA");
602 G4String excep =
"G4LEDATA environment variable not set!";
603 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
611 std::ostringstream ost;
613 ost << path <<
"/penelope/rayleigh/pdgra" << Z <<
".p08";
615 ost << path <<
"/penelope/rayleigh/pdgra0" << Z <<
".p08";
616 std::ifstream file(ost.str().c_str());
620 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
625 file >> readZ >> nPoints;
627 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
630 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
631 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
637 for (
size_t i=0;i<nPoints;i++)
639 file >> ene >>
f1 >>
f2 >> xs;
643 theVec->
PutValue(i,std::log(ene),std::log(xs));
644 if (file.eof() && i != (nPoints-1))
647 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
648 ed <<
"Found less than " << nPoints <<
"entries " <<
G4endl;
649 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
655 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
656 "em2044",
FatalException,
"Unable to allocate the atomic cross section table");
666 std::ostringstream ost2;
668 ost2 << path <<
"/penelope/rayleigh/pdaff" << Z <<
".p08";
670 ost2 << path <<
"/penelope/rayleigh/pdaff0" << Z <<
".p08";
671 file.open(ost2.str().c_str());
675 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
678 file >> readZ >> nPoints;
680 if (readZ != Z || nPoints <= 0 || nPoints >= 5000)
683 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
684 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
694 for (
size_t i=0;i<nPoints;i++)
696 file >> q >> ff >> incoh;
703 if (file.eof() && i != (nPoints-1))
706 ed <<
"Corrupted data file for Z=" << Z <<
G4endl;
707 ed <<
"Found less than " << nPoints <<
"entries " <<
G4endl;
708 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
714 G4Exception(
"G4PenelopeRayleighModel::ReadDataFile()",
716 "Unable to allocate the atomicFormFactor data table");
733 G4double logQSquared = (QSquared>1e-10) ? std::log(QSquared) : -23.;
744 ed <<
"Unable to retrieve F squared table for " << mat->
GetName() <<
G4endl;
745 G4Exception(
"G4PenelopeRayleighModel::GetFSquared()",
749 if (logQSquared < -20)
752 f2 = std::exp(logf2);
754 else if (logQSquared > maxlogQ2)
760 f2 = std::exp(logf2);
766 G4cout <<
"Q^2 = " << QSquared <<
" (units of 1/(m_e*c); F^2 = " << f2 <<
G4endl;
778 const size_t np = 150;
790 size_t nReducedPoints = np/4;
795 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
797 "Too few points to initialize the sampling algorithm");
799 if (q2min > (q2max-1e-10))
801 G4cout <<
"q2min= " << q2min <<
" q2max= " << q2max <<
G4endl;
802 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
804 "Too narrow grid to initialize the sampling algorithm");
817 const G4int nip = 51;
821 for (
size_t i=1;i<NUNIF-1;i++)
829 G4cout <<
"Vector x has " << x->size() <<
" points, while NUNIF = " << NUNIF <<
G4endl;
838 for (
size_t i=0;i<NUNIF-1;i++)
847 for (
G4int k=0;k<nip;k++)
851 pdfi->push_back(pdfk);
857 pdfih->push_back(pdfIK);
865 for (
G4int k=1;k<nip;k++)
868 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
869 sumi->push_back(next);
872 G4double lastIntegral = (*sumi)[sumi->size()-1];
873 area->push_back(lastIntegral);
876 for (
size_t k=0;k<sumi->size();k++)
880 if ((*pdfi)[0] < 1e-35)
881 (*pdfi)[0] = 1e-5*pdfmax;
882 if ((*pdfi)[pdfi->size()-1] < 1e-35)
883 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
887 G4double B_temp = 1.0-1.0/(pli*pui*dx*dx);
888 G4double A_temp = (1.0/(pli*dx))-1.0-B_temp;
889 G4double C_temp = 1.0+A_temp+B_temp;
898 a->push_back(A_temp);
899 b->push_back(B_temp);
900 c->push_back(C_temp);
912 for (
G4int k=0;k<nip;k++)
915 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
916 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
917 if (k == 0 || k == nip-1)
918 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
920 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
925 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
941 (*x)[x->size()-1] = q2max;
948 if (x->size() != NUNIF || a->size() != NUNIF ||
949 err->size() != NUNIF || area->size() != NUNIF)
952 ed <<
"Problem in building the Table for Sampling: array dimensions do not match" <<
G4endl;
953 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
965 for (
size_t i=0;i<err->size()-2;i++)
968 if ((*err)[i] > maxError)
970 maxError = (*err)[i];
976 G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]);
978 x->insert(x->begin()+iErrMax+1,newx);
980 area->insert(area->begin()+iErrMax+1,0.);
981 a->insert(a->begin()+iErrMax+1,0.);
982 b->insert(b->begin()+iErrMax+1,0.);
983 c->insert(c->begin()+iErrMax+1,0.);
984 err->insert(err->begin()+iErrMax+1,0.);
987 for (
size_t i=iErrMax;i<=iErrMax+1;i++)
994 G4double dxLocal = (*x)[i+1]-(*x)[i];
997 for (
G4int k=0;k<nip;k++)
1001 pdfi->push_back(pdfk);
1007 pdfih->push_back(pdfIK);
1014 sumi->push_back(0.);
1015 for (
G4int k=1;k<nip;k++)
1018 G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]);
1019 sumi->push_back(next);
1021 G4double lastIntegral = (*sumi)[sumi->size()-1];
1022 (*area)[i] = lastIntegral;
1026 for (
size_t k=0;k<sumi->size();k++)
1030 if ((*pdfi)[0] < 1e-35)
1031 (*pdfi)[0] = 1e-5*pdfmax;
1032 if ((*pdfi)[pdfi->size()-1] < 1e-35)
1033 (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax;
1037 G4double B_temp = 1.0-1.0/(pli*pui*dxLocal*dxLocal);
1038 G4double A_temp = (1.0/(pli*dxLocal))-1.0-B_temp;
1039 G4double C_temp = 1.0+A_temp+B_temp;
1060 for (
G4int k=0;k<nip;k++)
1063 G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/
1064 ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i]));
1065 if (k == 0 || k == nip-1)
1066 (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]);
1068 (*err)[i] += std::fabs(pap-(*pdfi)[k]);
1073 if ((*err)[i] > 0.1*(*area)[i] && icase == 1)
1086 }
while(x->size() < np);
1088 if (x->size() != np || a->size() != np ||
1089 err->size() != np || area->size() != np)
1091 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1093 "Problem in building the extended Table for Sampling: array dimensions do not match ");
1100 for (
size_t i=0;i<np-1;i++)
1104 for (
size_t i=0;i<np-1;i++)
1108 errMax =
std::max(errMax,(*err)[i]);
1114 for (
size_t i=0;i<np-1;i++)
1117 PAC->push_back(previous+(*area)[i]);
1119 (*PAC)[PAC->size()-1] = 1.;
1126 std::vector<size_t> *ITTL =
new std::vector<size_t>;
1127 std::vector<size_t> *ITTU =
new std::vector<size_t>;
1130 for (
size_t i=0;i<np;i++)
1138 for (
size_t i=1;i<(np-1);i++)
1142 for (
size_t j=(*ITTL)[i-1];j<np && !found;j++)
1144 if ((*PAC)[j] > ptst)
1152 (*ITTU)[ITTU->size()-2] = ITTU->size()-1;
1153 (*ITTU)[ITTU->size()-1] = ITTU->size()-1;
1154 (*ITTL)[ITTL->size()-1] = ITTU->size()-2;
1156 if (ITTU->size() != np || ITTU->size() != np)
1158 G4Exception(
"G4PenelopeRayleighModel::InitializeSamplingAlgorithm()",
1160 "Problem in building the Limit Tables for Sampling: array dimensions do not match");
1168 for (
size_t i=0;i<np;i++)
1170 theTable->
AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]);
1175 G4cout <<
"*************************************************************************" <<
1177 G4cout <<
"Sampling table for Penelope Rayleigh scattering in " << mat->
GetName() <<
G4endl;
1206 G4cout <<
"G4PenelopeRayleighModel::BuildPMaxTable" <<
G4endl;
1207 G4cout <<
"Going to instanziate the pMaxTable !" <<
G4endl;
1209 pMaxTable =
new std::map<const G4Material*,G4PhysicsFreeVector*>;
1218 G4Exception(
"G4PenelopeRayleighModel::GetPMaxTable()",
1220 "SamplingTable is not properly instantiated");
1228 ed <<
"Sampling table for material " << mat->
GetName() <<
" not found";
1229 G4Exception(
"G4PenelopeRayleighModel::GetPMaxTable()",
1241 const size_t nip = 51;
1246 G4double Qm = 2.0*energy/electron_mass_c2;
1257 size_t lowerBound = 0;
1258 size_t upperBound = tablePoints-1;
1259 while (lowerBound <= upperBound)
1261 size_t midBin = (lowerBound + upperBound)/2;
1262 if( Qm2 < theTable->GetX(midBin))
1263 { upperBound = midBin-1; }
1265 { lowerBound = midBin+1; }
1275 for (
size_t k=0;k<nip;k++)
1279 (theTable->
GetX(upperBound+1)-Q1);
1284 if (std::fabs(con1) > 1.0e-16*std::fabs(con2))
1285 etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1;
1289 (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/
1290 ((1.0-theB*etap*etap)*ci*(theTable->
GetX(upperBound+1)-Q1));
1291 fun->push_back(theFun);
1299 (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS;
1300 sum->push_back(secondPoint);
1301 for (
size_t hh=2;hh<nip-1;hh++)
1304 G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])-
1305 (*fun)[hh+1]-(*fun)[hh-2])*HCONS;
1306 sum->push_back(next);
1308 G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]-
1309 (*fun)[nip-3])*CONS;
1310 sum->push_back(last);
1311 thePMax = thePAC + (*sum)[sum->size()-1];
1322 thePMax = theTable->
GetPAC(0);
1326 theVec->
PutValue(ie,energy,thePMax);
1329 pMaxTable->insert(std::make_pair(mat,theVec));
1338 G4cout <<
"*****************************************************************" <<
G4endl;
1342 G4cout <<
"*****************************************************************" <<
G4endl;
G4double LowEnergyLimit() const
std::vector< G4Element * > G4ElementVector
void PutValue(size_t binNumber, G4double binValue, G4double dataValue)
std::ostringstream G4ExceptionDescription
virtual ~G4PenelopeRayleighModel()
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double GetB(size_t index)
std::map< const G4Material *, G4PhysicsFreeVector * > * logFormFactorTable
const G4String & GetName() const
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
G4double fIntrinsicHighEnergyLimit
G4double GetPAC(size_t index)
size_t GetVectorLength() const
void BuildFormFactorTable(const G4Material *)
G4double GetLowEdgeEnergy(size_t binNumber) const
const G4ElementVector * GetElementVector() const
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4MUTEX_INITIALIZER
G4double GetFSquared(const G4Material *, const G4double)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
size_t GetNumberOfStoredPoints()
G4GLOB_DLL std::ostream G4cout
size_t GetTableSize() const
void GetPMaxTable(const G4Material *)
void InitializeSamplingAlgorithm(const G4Material *)
const G4ThreeVector & GetMomentumDirection() const
std::map< G4int, G4PhysicsFreeVector * > * atomicFormFactor
static const double twopi
void DumpFormFactorTable(const G4Material *)
G4double Value(G4double theEnergy, size_t &lastidx) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
static G4ProductionCutsTable * GetProductionCutsTable()
G4double SampleValue(G4double rndm)
static const G4double factor
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
G4double fIntrinsicLowEnergyLimit
const G4double x[NPOINTSGL]
G4double GetA(size_t index)
G4DataVector logEnergyGridPMax
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4double GetX(size_t index)
std::map< G4int, G4PhysicsFreeVector * > * logAtomicCrossSection
G4ParticleChangeForGamma * fParticleChange
void SetProposedKineticEnergy(G4double proposedKinEnergy)
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
std::map< const G4Material *, G4PhysicsFreeVector * > * pMaxTable
size_t GetNumberOfElements() const
void ProposeTrackStatus(G4TrackStatus status)
G4PenelopeRayleighModel(const G4ParticleDefinition *p=0, const G4String &processName="PenRayleigh")
void AddPoint(G4double x0, G4double pac0, G4double a0, G4double b0, size_t ITTL0, size_t ITTU0)
G4ThreeVector G4ParticleMomentum
const G4double * GetFractionVector() const
G4DataVector logQSquareGrid
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
std::map< const G4Material *, G4PenelopeSamplingData * > * samplingTable
const G4ParticleDefinition * fParticle
G4ParticleChangeForGamma * GetParticleChangeForGamma()
const G4Material * GetMaterial() const
void SetParticle(const G4ParticleDefinition *)