65 , fGammaCutInKineticEnergy(nullptr)
66 , fAngleDistrTable(nullptr)
67 , fEnergyDistrTable(nullptr)
68 , fAngleForEnergyTable(nullptr)
69 , fPlatePhotoAbsCof(nullptr)
70 , fGasPhotoAbsCof(nullptr)
109 G4cout <<
"### G4VXTRenergyLoss: the number of TR radiator plates = "
113 G4Exception(
"G4VXTRenergyLoss::G4VXTRenergyLoss()",
"VXTRELoss01",
186 out <<
"Base class for 'fast' parameterisation model describing X-ray "
188 "radiation. Angular distribution is very rough.\n";
205 G4double charge, chargeSq, massRatio, TkinScaled;
217 gamma = 1.0 + kinEnergy / mass;
223 if(std::fabs(gamma -
fGamma) < 0.05 * gamma)
228 chargeSq = charge * charge;
230 TkinScaled = kinEnergy * massRatio;
232 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
234 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
245 sigma = (*(*fEnergyDistrTable)(iPlace))(0) * chargeSq;
252 W1 = (E2 - TkinScaled) * W;
253 W2 = (TkinScaled - E1) * W;
254 sigma = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 +
280 G4Exception(
"G4VXTRenergyLoss::BuildPhysicsTable",
"Notification",
281 JustWarning,
"XTR initialisation for neutral particle ?!");
290 <<
"Build angle for energy distribution according the current radiator"
301 G4int iTkin, iTR, iPlace;
332 G4cout <<
"Lorentz Factor"
334 <<
"XTR photon number" <<
G4endl;
337 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
356 for(iTR =
fBinTR - 2; iTR >= 0; --iTR)
360 energySum += radiatorCof *
fCofTR *
384 G4cout <<
"total time for build X-ray TR energy loss tables = "
431 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
443 for(iTR = 0; iTR <
fBinTR; ++iTR)
455 for(i =
fBinTR - 2; i >= 0; --i)
475 G4cout <<
"total time for build X-ray TR angle for energy loss tables = "
511 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
528 for(iTR = 0; iTR <
fBinTR; ++iTR)
543 G4cout <<
"total time for build XTR angle for given energy tables = "
555 G4double theta = 0., result, tmp = 0., cof1, cof2, cofMin, cofPHC,
578 <<
"; tmp = " << 0. <<
"; angleSum = " << angleSum <<
G4endl;
581 for(iTheta =
n - 1; iTheta >= 1; --iTheta)
583 k = iTheta - 1 +
kMin;
585 result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2);
586 tmp = std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result;
593 else if(iTheta ==
n - 1)
603 G4cout <<
"iTheta = " << iTheta <<
"; k = " << k
604 <<
"; theta = " << std::sqrt(theta) *
fGamma <<
"; tmp = " << tmp
605 <<
"; angleSum = " << angleSum <<
G4endl;
607 angleVector->
PutValue(iTheta, theta, angleSum);
617 G4cout <<
"iTheta = " << iTheta <<
"; theta = " << std::sqrt(theta) *
fGamma
618 <<
"; tmp = " << tmp <<
"; angleSum = " << angleSum <<
G4endl;
620 angleVector->
PutValue(iTheta, theta, angleSum);
630 G4int iTkin, iTR, iPlace;
654 G4cout <<
"Lorentz Factor"
656 <<
"XTR photon number" <<
G4endl;
659 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
686 for(iTR =
fBinTR - 2; iTR >= 0; --iTR)
688 angleSum += radiatorCof *
fCofTR *
693 angleVector->
PutValue(iTR, angleSum);
707 G4cout <<
"total time for build X-ray TR angle tables = "
722 G4double energyTR, theta, theta2, phi, dirX, dirY, dirZ;
728 G4cout <<
"Start of G4VXTRenergyLoss::PostStepDoIt " <<
G4endl;
729 G4cout <<
"name of current material = "
737 G4cout <<
"Go out from G4VXTRenergyLoss::PostStepDoIt: wrong volume "
750 G4double gamma = 1.0 + kinEnergy / mass;
757 G4double TkinScaled = kinEnergy * massRatio;
762 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
764 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
772 G4cout <<
"Go out from G4VXTRenergyLoss::PostStepDoIt:iTkin = " << iTkin
791 theta = std::sqrt(theta2);
803 dirX = std::sin(theta) * std::cos(phi);
804 dirY = std::sin(theta) * std::sin(phi);
805 dirZ = std::cos(theta);
832 G4cout <<
"distance to exit = " << distance /
mm <<
" mm" <<
G4endl;
835 startTime += distance /
c_light;
881 static constexpr G4int iMax = 8;
884 G4double lim[iMax] = { 0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0 };
886 for(i = 0; i < iMax; ++i)
894 for(i = 0; i < iMax - 1; ++i)
896 angleSum += integral.Legendre96(
918 G4double sum = 0., tmp1, tmp2, tmp = 0., cof1, cof2, cofMin, cofPHC, energy1,
927 cofMin = std::sqrt(cof1 * cof2);
939 tmp2 = std::sqrt(tmp1 * tmp1 - cof1 * cof2);
940 energy1 = (tmp1 + tmp2) / cof1;
941 energy2 = (tmp1 - tmp2) / cof1;
943 for(i = 0; i < 2; ++i)
953 tmp2 = std::sin(tmp1);
954 tmp = energy1 * tmp2 * tmp2;
959 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2);
960 tmp1 = cof1 / (4. *
hbarc) - cof2 / (4. *
hbarc * energy1 * energy1);
961 tmp2 = std::abs(tmp1);
976 tmp2 = std::sin(tmp1);
977 tmp = energy2 * tmp2 * tmp2;
982 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2);
983 tmp1 = cof1 / (4. *
hbarc) - cof2 / (4. *
hbarc * energy2 * energy2);
984 tmp2 = std::abs(tmp1);
1006 lambda = 1.0 / gamma / gamma + varAngle +
fSigma1 / omega / omega;
1016 G4double cof, length, delta, real_v, image_v;
1020 cof = 1.0 / (1.0 + delta * delta);
1022 real_v = length * cof;
1023 image_v = real_v * delta;
1048 omega2 = omega * omega;
1049 omega3 = omega2 * omega;
1050 omega4 = omega2 * omega2;
1053 G4double cross = SandiaCof[0] / omega + SandiaCof[1] / omega2 +
1054 SandiaCof[2] / omega3 + SandiaCof[3] / omega4;
1064 lambda = 1.0 / gamma / gamma + varAngle +
fSigma2 / omega / omega;
1074 G4double cof, length, delta, real_v, image_v;
1078 cof = 1.0 / (1.0 + delta * delta);
1080 real_v = length * cof;
1081 image_v = real_v * delta;
1105 omega2 = omega * omega;
1106 omega3 = omega2 * omega;
1107 omega4 = omega2 * omega2;
1110 G4double cross = SandiaCof[0] / omega + SandiaCof[1] / omega2 +
1111 SandiaCof[2] / omega3 + SandiaCof[3] / omega4;
1129 std::ofstream outPlate(
"plateZmu.dat", std::ios::out);
1130 outPlate.setf(std::ios::scientific, std::ios::floatfield);
1135 varAngle = 1 / gamma / gamma;
1137 G4cout <<
"energy, keV" <<
"\t" <<
"Zmu for plate" <<
G4endl;
1138 for(i = 0; i < 100; ++i)
1140 omega = (1.0 + i) *
keV;
1145 outPlate << omega /
keV <<
"\t\t"
1166 std::ofstream outGas(
"gasZmu.dat", std::ios::out);
1167 outGas.setf(std::ios::scientific, std::ios::floatfield);
1171 varAngle = 1 / gamma / gamma;
1173 G4cout <<
"energy, keV" <<
"\t" <<
"Zmu for gas" <<
G4endl;
1174 for(i = 0; i < 100; ++i)
1176 omega = (1.0 + i) *
keV;
1181 outGas << omega /
keV <<
"\t\t"
1191 G4int i, numberOfElements;
1192 G4double xSection = 0., nowZ, sumZ = 0.;
1195 numberOfElements = (*theMaterialTable)[
fMatIndex1]->GetNumberOfElements();
1197 for(i = 0; i < numberOfElements; ++i)
1199 nowZ = (*theMaterialTable)[
fMatIndex1]->GetElement(i)->GetZ();
1204 xSection *= (*theMaterialTable)[
fMatIndex1]->GetElectronDensity();
1212 G4int i, numberOfElements;
1213 G4double xSection = 0., nowZ, sumZ = 0.;
1216 numberOfElements = (*theMaterialTable)[
fMatIndex2]->GetNumberOfElements();
1218 for(i = 0; i < numberOfElements; ++i)
1220 nowZ = (*theMaterialTable)[
fMatIndex2]->GetElement(i)->GetZ();
1225 xSection *= (*theMaterialTable)[
fMatIndex2]->GetElectronDensity();
1236 return CrossSection;
1237 if(GammaEnergy < 0.1 *
keV)
1238 return CrossSection;
1239 if(GammaEnergy > (100. *
GeV /
Z))
1240 return CrossSection;
1242 static constexpr G4double a = 20.0;
1243 static constexpr G4double b = 230.0;
1244 static constexpr G4double c = 440.0;
1247 d3 = 6.7527 *
barn, d4 = -1.9798e+1 *
barn,
1250 f1 = -3.9178e-7 *
barn, f2 = 6.8241e-5 *
barn,
1251 f3 = 6.0480e-5 *
barn, f4 = 3.0274e-4 *
barn;
1264 p1Z * std::log(1. + 2. * X) / X +
1265 (p2Z + p3Z * X + p4Z * X * X) / (1. + a * X + b * X * X + c * X * X * X);
1268 if(GammaEnergy <
T0)
1273 p1Z * std::log(1. + 2. * X) / X +
1274 (p2Z + p3Z * X + p4Z * X * X) / (1. + a * X + b * X * X + c * X * X * X);
1275 G4double c1 = -
T0 * (sigma - CrossSection) / (CrossSection * dT0);
1278 c2 = 0.375 - 0.0556 * std::log(
Z);
1280 CrossSection *= std::exp(-y * (c1 + c2 * y));
1282 return CrossSection;
1296 G4double formationLength1, formationLength2;
1301 return (varAngle /
energy) * (formationLength1 - formationLength2) *
1302 (formationLength1 - formationLength2);
1359 std::ofstream outEn(
"numberE.dat", std::ios::out);
1360 outEn.setf(std::ios::scientific, std::ios::floatfield);
1362 std::ofstream outAng(
"numberAng.dat", std::ios::out);
1363 outAng.setf(std::ios::scientific, std::ios::floatfield);
1365 for(iTkin = 0; iTkin <
fTotBin; ++iTkin)
1369 numberE = (*(*fEnergyDistrTable)(iTkin))(0);
1371 G4cout << gamma <<
"\t\t" << numberE <<
"\t" <<
G4endl;
1373 outEn << gamma <<
"\t\t" << numberE <<
G4endl;
1383 G4int iTransfer, iPlace;
1392 for(iTransfer = 0;; ++iTransfer)
1403 W = 1.0 / (E2 - E1);
1404 W1 = (E2 - scaledTkin) * W;
1405 W2 = (scaledTkin - E1) * W;
1407 position = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 +
1411 for(iTransfer = 0;; ++iTransfer)
1433 result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1437 y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer - 1);
1438 y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer);
1440 x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1);
1441 x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1471 for(iTR = 0; iTR <
fBinTR; ++iTR)
1473 if(energyXTR < fXTREnergyVector->GetLowEdgeEnergy(iTR))
1482 for(iAngle = 0;; ++iAngle)
1501 if( iTransfer == 0 )
1504 result = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1508 y1 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer - 1);
1509 y2 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer);
1511 x1 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1);
1512 x2 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1514 if(x1 == x2) result = x2;
1520 result = x1 + (
position - y1) * (x2 - x1) / (y2 - y1);
static const G4double e4[47]
static const G4double e1[44]
static const G4double e2[44]
static const G4double e3[45]
G4double condition(const G4ErrorSymMatrix &m)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::vector< G4Material * > G4MaterialTable
static const G4double T0[78]
static constexpr double twopi
static constexpr double barn
static constexpr double mm
static constexpr double keV
static constexpr double eV
static constexpr double GeV
static constexpr double pi
static constexpr double cm
static const G4double angle[DIMMOTT]
std::complex< G4double > G4complex
G4GLOB_DLL std::ostream G4cout
Hep3Vector & rotateUz(const Hep3Vector &)
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4VSolid * GetSolid() const
G4Material * GetMaterial() const
G4SandiaTable * GetSandiaTable() const
G4double GetElectronDensity() const
static G4MaterialTable * GetMaterialTable()
const G4String & GetName() const
void AddSecondary(G4Track *aSecondary)
void ProposeEnergy(G4double finalEnergy)
virtual void Initialize(const G4Track &)
G4double GetPDGMass() const
G4double GetPDGCharge() const
void PutValue(const std::size_t index, const G4double e, const G4double value)
static G4int GetModelID(const G4int modelIndex)
void insertAt(std::size_t, G4PhysicsVector *)
G4double GetLowEdgeEnergy(const std::size_t index) const
void PutValue(const std::size_t index, const G4double value)
G4double GetSandiaCofForMaterial(G4int, G4int) const
const G4VTouchable * GetTouchable() const
G4double GetGlobalTime() const
const G4ThreeVector & GetPosition() const
const G4TouchableHandle & GetTouchableHandle() const
G4StepPoint * GetPostStepPoint() const
G4double GetUserElapsed() const
G4VPhysicalVolume * GetVolume() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
const G4DynamicParticle * GetDynamicParticle() const
void SetParentID(const G4int aValue)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
void SetNumberOfSecondaries(G4int totSecondaries)
G4LogicalVolume * GetLogicalVolume() const
void SetProcessSubType(G4int)
G4VParticleChange * pParticleChange
const G4String & GetProcessName() const
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const =0
virtual const G4ThreeVector & GetTranslation(G4int depth=0) const =0
virtual const G4RotationMatrix * GetRotation(G4int depth=0) const =0
G4double GetPlateLinearPhotoAbs(G4double)
G4double AngleSpectralXTRdEdx(G4double energy)
G4double XTRNAngleDensity(G4double varAngle)
G4PhysicsTable * fEnergyDistrTable
void ComputeGasPhotoAbsCof()
G4double GetAngleXTR(G4int iTR, G4double position, G4int iAngle)
virtual G4bool IsApplicable(const G4ParticleDefinition &) override
G4PhysicsTable * fAngleForEnergyTable
static constexpr G4double fMaxProtonTkin
G4PhysicsLogVector * fProtonEnergyVector
G4PhysicsLogVector * fXTREnergyVector
G4double GetGasFormationZone(G4double, G4double, G4double)
G4SandiaTable * fPlatePhotoAbsCof
G4double OneBoundaryXTRNdensity(G4double energy, G4double gamma, G4double varAngle) const
G4PhysicsFreeVector * GetAngleVector(G4double energy, G4int n)
G4PhysicsTable * fAngleDistrTable
G4double GetRandomAngle(G4double energyXTR, G4int iTkin)
G4complex GetPlateComplexFZ(G4double, G4double, G4double)
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
G4complex OneInterfaceXTRdEdx(G4double energy, G4double gamma, G4double varAngle)
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
G4double GetXTRenergy(G4int iPlace, G4double position, G4int iTransfer)
G4VXTRenergyLoss(G4LogicalVolume *anEnvelope, G4Material *, G4Material *, G4double, G4double, G4int, const G4String &processName="XTRenergyLoss", G4ProcessType type=fElectromagnetic)
void BuildGlobalAngleTable()
G4double XTRNAngleSpectralDensity(G4double energy)
G4double SpectralAngleXTRdEdx(G4double varAngle)
G4double GetXTRrandomEnergy(G4double scaledTkin, G4int iTkin)
G4SandiaTable * fGasPhotoAbsCof
static constexpr G4double fCofTR
void GetPlateZmuProduct()
void BuildAngleForEnergyBank()
G4LogicalVolume * fEnvelope
std::vector< G4PhysicsTable * > fAngleBank
static constexpr G4double fPlasmaCof
virtual ~G4VXTRenergyLoss()
virtual G4double SpectralXTRdEdx(G4double energy)
G4double GetComptonPerAtom(G4double, G4double)
void GetNumberOfPhotons()
static constexpr G4double fMinProtonTkin
G4double GetGasCompton(G4double)
G4double XTRNSpectralAngleDensity(G4double varAngle)
virtual void ProcessDescription(std::ostream &) const override
void ComputePlatePhotoAbsCof()
G4double GetPlateFormationZone(G4double, G4double, G4double)
G4double GetGasLinearPhotoAbs(G4double)
G4double XTRNSpectralDensity(G4double energy)
virtual G4double GetStackFactor(G4double energy, G4double gamma, G4double varAngle)
G4ParticleDefinition * fPtrGamma
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
G4ParticleChange fParticleChange
G4complex GetGasComplexFZ(G4double, G4double, G4double)
G4double GetPlateCompton(G4double)
G4double AngleXTRdEdx(G4double varAngle)
static constexpr double keV
ThreeVector shoot(const G4int Ap, const G4int Af)
G4double energy(const ThreeVector &p, const G4double m)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static const G4double Z1[5]