70 fGammaCutInKineticEnergy(nullptr),
72 fAngleDistrTable(nullptr),
73 fEnergyDistrTable(nullptr),
74 fPlatePhotoAbsCof(nullptr),
75 fGasPhotoAbsCof(nullptr),
76 fAngleForEnergyTable(nullptr)
116 G4cout<<
"### G4VXTRenergyLoss: the number of TR radiator plates = "
120 G4Exception(
"G4VXTRenergyLoss::G4VXTRenergyLoss()",
"VXTRELoss01",
210 G4double charge, chargeSq, massRatio, TkinScaled;
221 gamma = 1.0 + kinEnergy/mass;
227 if ( std::fabs( gamma -
fGamma ) < 0.05*gamma ) lambda =
fLambda;
231 chargeSq = charge*charge;
233 TkinScaled = kinEnergy*massRatio;
235 for(iTkin = 0; iTkin <
fTotBin; iTkin++)
237 if( TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
break;
241 if(iTkin == 0) lambda =
DBL_MAX;
246 sigma = (*(*fEnergyDistrTable)(iPlace))(0)*chargeSq;
253 W1 = (E2 - TkinScaled)*W;
254 W2 = (TkinScaled - E1)*W;
255 sigma = ( (*(*fEnergyDistrTable)(iPlace ))(0)*W1 +
260 else lambda = 1./sigma;
282 "XTR initialisation for neutral particle ?!" );
290 G4cout<<
"Build angle for energy distribution according the current radiator"
304 G4int iTkin, iTR, iPlace;
330 G4cout<<
"Lorentz Factor"<<
"\t"<<
"XTR photon number"<<
G4endl;
333 for( iTkin = 0; iTkin <
fTotBin; iTkin++ )
353 for( iTR =
fBinTR - 2; iTR >= 0; iTR-- )
383 G4cout<<
"total time for build X-ray TR energy loss tables = "
428 for( iTkin = 0; iTkin <
fTotBin; iTkin++ )
443 for( iTR = 0; iTR <
fBinTR; iTR++ )
451 angleVector ->
PutValue(fBinTR - 1, angleSum);
453 for( i = fBinTR - 2; i >= 0; i-- )
462 angleVector ->
PutValue(i, angleSum);
473 G4cout<<
"total time for build X-ray TR angle for energy loss tables = "
506 G4cout<<
"Lorentz Factor"<<
"\t"<<
"XTR photon number"<<
G4endl;
509 for( iTkin = 0; iTkin <
fTotBin; iTkin++ )
527 for( iTR = 0; iTR <
fBinTR; iTR++ )
546 G4cout<<
"total time for build XTR angle for given energy tables = "
574 kMin =
G4int(cofMin);
575 if (cofMin > kMin) kMin++;
581 G4cout<<
"n-1 = "<<n-1<<
"; theta = "
584 <<
"; angleSum = "<<angleSum<<
G4endl;
588 for( iTheta = n - 1; iTheta >= 1; iTheta-- )
591 k = iTheta - 1 +
kMin;
595 result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
597 tmp = std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/
result;
599 if( k == kMin && kMin ==
G4int(cofMin) )
603 else if(iTheta == n-1);
608 theta = std::abs(k-cofMin)*cofPHC/energy/(fPlateThick +
fGasThick);
612 G4cout<<
"iTheta = "<<iTheta<<
"; k = "<<k<<
"; theta = "
613 <<std::sqrt(theta)*
fGamma<<
"; tmp = "
615 <<
"; angleSum = "<<angleSum<<
G4endl;
617 angleVector->
PutValue( iTheta, theta, angleSum );
626 G4cout<<
"iTheta = "<<iTheta<<
"; theta = "
627 <<std::sqrt(theta)*
fGamma<<
"; tmp = "
629 <<
"; angleSum = "<<angleSum<<
G4endl;
631 angleVector->
PutValue( iTheta, theta, angleSum );
642 G4int iTkin, iTR, iPlace;
662 G4cout<<
"Lorentz Factor"<<
"\t"<<
"XTR photon number"<<
G4endl;
665 for( iTkin = 0; iTkin <
fTotBin; iTkin++ )
691 for( iTR =
fBinTR - 2; iTR >= 0; iTR-- )
699 angleVector ->
PutValue(iTR,angleSum);
717 G4cout<<
"total time for build X-ray TR angle tables = "
735 G4double energyTR, theta,theta2, phi, dirX, dirY, dirZ;
742 G4cout<<
"Start of G4VXTRenergyLoss::PostStepDoIt "<<
G4endl;
743 G4cout<<
"name of current material = "
750 G4cout<<
"Go out from G4VXTRenergyLoss::PostStepDoIt: wrong volume "<<
G4endl;
763 G4double gamma = 1.0 + kinEnergy/mass;
770 G4double TkinScaled = kinEnergy*massRatio;
775 for( iTkin = 0; iTkin <
fTotBin; iTkin++ )
777 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
break;
785 G4cout<<
"Go out from G4VXTRenergyLoss::PostStepDoIt:iTkin = "<<iTkin<<
G4endl;
805 if( theta2 > 0.) theta = std::sqrt(theta2);
812 if( theta >= 0.1 ) theta = 0.1;
818 dirX = std::sin(theta)*std::cos(phi);
819 dirY = std::sin(theta)*std::sin(phi);
820 dirZ = std::cos(theta);
827 directionTR, energyTR);
847 position += distance*directionTR;
851 startTime, position );
895 if(result < 0.0) result = 0.0;
908 G4double lim[8] = { 0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0 };
934 for( i = 0; i < iMax-1; i++ )
953 if(result < 0) result = 0.0;
966 G4double sum = 0., tmp1, tmp2,
tmp=0., cof1, cof2, cofMin, cofPHC, energy1, energy2;
976 cofMin = std::sqrt(cof1*cof2);
979 kMin =
G4int(cofMin);
980 if (cofMin > kMin) kMin++;
986 for( k = kMin; k <=
kMax; k++ )
989 tmp2 = std::sqrt(tmp1*tmp1-cof1*cof2);
990 energy1 = (tmp1+tmp2)/cof1;
991 energy2 = (tmp1-tmp2)/cof1;
993 for(i = 0; i < 2; i++)
1001 tmp2 = std::sin(tmp1);
1002 tmp = energy1*tmp2*tmp2;
1004 tmp1 = hbarc*energy1/( energy1*energy1*(1./
fGamma/
fGamma + varAngle) + fSigma2 );
1005 tmp *= (tmp1-tmp2)*(tmp1-tmp2);
1006 tmp1 = cof1/(4.*
hbarc) - cof2/(4.*hbarc*energy1*energy1);
1007 tmp2 = std::abs(tmp1);
1009 if(tmp2 > 0.) tmp /= tmp2;
1018 tmp2 = std::sin(tmp1);
1019 tmp = energy2*tmp2*tmp2;
1021 tmp1 = hbarc*energy2/( energy2*energy2*(1./
fGamma/
fGamma + varAngle) + fSigma2 );
1022 tmp *= (tmp1-tmp2)*(tmp1-tmp2);
1023 tmp1 = cof1/(4.*
hbarc) - cof2/(4.*hbarc*energy2*energy2);
1024 tmp2 = std::abs(tmp1);
1026 if(tmp2 > 0.) tmp /= tmp2;
1035 result /= hbarc*
hbarc;
1057 lambda = 1.0/gamma/gamma + varAngle +
fSigma1/omega/omega;
1070 G4double cof, length,delta, real_v, image_v;
1074 cof = 1.0/(1.0 + delta*delta);
1076 real_v = length*cof;
1077 image_v = real_v*delta;
1109 omega2 = omega*omega;
1110 omega3 = omega2*omega;
1111 omega4 = omega2*omega2;
1114 G4double cross = SandiaCof[0]/omega + SandiaCof[1]/omega2 +
1115 SandiaCof[2]/omega3 + SandiaCof[3]/omega4;
1129 lambda = 1.0/gamma/gamma + varAngle +
fSigma2/omega/omega;
1143 G4double cof, length,delta, real_v, image_v;
1147 cof = 1.0/(1.0 + delta*delta);
1149 real_v = length*cof;
1150 image_v = real_v*delta;
1180 omega2 = omega*omega;
1181 omega3 = omega2*omega;
1182 omega4 = omega2*omega2;
1185 G4double cross = SandiaCof[0]/omega + SandiaCof[1]/omega2 +
1186 SandiaCof[2]/omega3 + SandiaCof[3]/omega4;
1210 std::ofstream outPlate(
"plateZmu.dat", std::ios::out );
1211 outPlate.setf( std::ios::scientific, std::ios::floatfield );
1216 varAngle = 1/gamma/gamma;
1221 omega = (1.0 + i)*
keV;
1248 std::ofstream outGas(
"gasZmu.dat", std::ios::out );
1249 outGas.setf( std::ios::scientific, std::ios::floatfield );
1253 varAngle = 1/gamma/gamma;
1258 omega = (1.0 + i)*
keV;
1273 G4int i, numberOfElements;
1274 G4double xSection = 0., nowZ, sumZ = 0.;
1277 numberOfElements = (*theMaterialTable)[
fMatIndex1]->GetNumberOfElements();
1279 for( i = 0; i < numberOfElements; i++ )
1281 nowZ = (*theMaterialTable)[
fMatIndex1]->GetElement(i)->GetZ();
1286 xSection *= (*theMaterialTable)[
fMatIndex1]->GetElectronDensity();
1297 G4int i, numberOfElements;
1298 G4double xSection = 0., nowZ, sumZ = 0.;
1301 numberOfElements = (*theMaterialTable)[
fMatIndex2]->GetNumberOfElements();
1303 for( i = 0; i < numberOfElements; i++ )
1305 nowZ = (*theMaterialTable)[
fMatIndex2]->GetElement(i)->GetZ();
1310 xSection *= (*theMaterialTable)[
fMatIndex2]->GetElectronDensity();
1322 if ( Z < 0.9999 )
return CrossSection;
1323 if ( GammaEnergy < 0.1*
keV )
return CrossSection;
1324 if ( GammaEnergy > (100.*
GeV/Z) )
return CrossSection;
1326 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
1330 e1= 1.9756e-5*
barn, e2=-1.0205e-2*
barn, e3=-7.3913e-2*
barn, e4= 2.7079e-2*
barn,
1333 G4double p1Z = Z*(d1 + e1*Z + f1*Z*
Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
1334 p3Z = Z*(d3 + e3*Z + f3*Z*
Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
1337 if (Z < 1.5) T0 = 40.0*
keV;
1340 CrossSection = p1Z*std::log(1.+2.*X)/X
1341 + (p2Z + p3Z*X + p4Z*X*
X)/(1. + a*X + b*X*X + c*X*X*X);
1345 if (GammaEnergy < T0)
1349 G4double sigma = p1Z*std::log(1.+2.*X)/X
1350 + (p2Z + p3Z*X + p4Z*X*
X)/(1. + a*X + b*X*X + c*X*X*X);
1351 G4double c1 = -T0*(sigma-CrossSection)/(CrossSection*dT0);
1353 if (Z > 1.5) c2 = 0.375-0.0556*std::log(Z);
1355 CrossSection *= std::exp(-y*(c1+c2*y));
1358 return CrossSection;
1376 G4double formationLength1, formationLength2;
1377 formationLength1 = 1.0/
1381 formationLength2 = 1.0/
1385 return (varAngle/energy)*(formationLength1 - formationLength2)
1386 *(formationLength1 - formationLength2);
1456 std::ofstream outEn(
"numberE.dat", std::ios::out );
1457 outEn.setf( std::ios::scientific, std::ios::floatfield );
1459 std::ofstream outAng(
"numberAng.dat", std::ios::out );
1460 outAng.setf( std::ios::scientific, std::ios::floatfield );
1462 for(iTkin=0;iTkin<
fTotBin;iTkin++)
1466 numberE = (*(*fEnergyDistrTable)(iTkin))(0);
1469 G4cout<<gamma<<
"\t\t"<<numberE<<
"\t"
1472 outEn<<gamma<<
"\t\t"<<numberE<<
G4endl;
1484 G4int iTransfer, iPlace;
1495 for(iTransfer=0;;iTransfer++)
1506 W1 = (E2 - scaledTkin)*W;
1507 W2 = (scaledTkin - E1)*W;
1509 position =( (*(*fEnergyDistrTable)(iPlace))(0)*W1 +
1514 for(iTransfer=0;;iTransfer++)
1524 if(transfer < 0.0 ) transfer = 0.0;
1541 result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1545 y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer-1);
1546 y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer);
1548 x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1549 x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1551 if ( x1 == x2 ) result =
x2;
1574 if (iTkin ==
fTotBin) iTkin--;
1578 for( iTR = 0; iTR <
fBinTR; iTR++ )
1580 if( energyXTR < fXTREnergyVector->GetLowEdgeEnergy(iTR) )
break;
1582 if (iTR == fBinTR) iTR--;
1584 position = (*(*fAngleForEnergyTable)(iTR))(0)*
G4UniformRand();
1586 for( iAngle = 0;; iAngle++)
1605 if( iTransfer == 0 )
1607 result = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1611 y1 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer-1);
1612 y2 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer);
1614 x1 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1);
1615 x2 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1617 if ( x1 == x2 ) result =
x2;
1623 result = x1 + (position -
y1)*(x2 - x1)/(y2 -
y1);
G4ParticleChange fParticleChange
G4double GetGasCompton(G4double)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double OneBoundaryXTRNdensity(G4double energy, G4double gamma, G4double varAngle) const
const G4VTouchable * GetTouchable() const
std::vector< ExP01TrackerHit * > a
virtual G4double SpectralXTRdEdx(G4double energy)
G4LogicalVolume * GetLogicalVolume() const
G4double GetLowEdgeEnergy(size_t binNumber) const
void BuildGlobalAngleTable()
Float_t y1[n_points_granero]
static G4MaterialTable * GetMaterialTable()
G4complex GetGasComplexFZ(G4double, G4double, G4double)
static constexpr double keV
G4ParticleDefinition * fPtrGamma
static constexpr double mm
Float_t x1[n_points_granero]
const G4ThreeVector & GetMomentumDirection() const
static constexpr double hbarc
G4double GetXTRenergy(G4int iPlace, G4double position, G4int iTransfer)
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
G4double GetPDGCharge() const
virtual const G4ThreeVector & GetTranslation(G4int depth=0) const =0
void GetPlateZmuProduct()
Hep3Vector & rotateUz(const Hep3Vector &)
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
virtual const G4RotationMatrix * GetRotation(G4int depth=0) const =0
void insertAt(size_t, G4PhysicsVector *)
G4double condition(const G4ErrorSymMatrix &m)
G4double GetPDGMass() const
void BuildAngleForEnergyBank()
const G4TouchableHandle & GetTouchableHandle() const
void AddSecondary(G4Track *aSecondary)
G4Material * GetMaterial() const
G4PhysicsLogVector * fXTREnergyVector
Float_t y2[n_points_geant4]
G4double AngleXTRdEdx(G4double varAngle)
static constexpr double proton_mass_c2
const G4String & GetName() const
static constexpr double TeV
G4double SpectralAngleXTRdEdx(G4double varAngle)
std::complex< G4double > G4complex
G4PhysicsTable * fEnergyDistrTable
G4ParticleDefinition * GetDefinition() const
G4VParticleChange * pParticleChange
void ComputeGasPhotoAbsCof()
virtual void Initialize(const G4Track &)
G4double GetGasFormationZone(G4double, G4double, G4double)
static const G4double T0[78]
static constexpr double electron_mass_c2
const G4ThreeVector & GetPosition() const
virtual G4bool IsApplicable(const G4ParticleDefinition &) override
static constexpr double twopi
G4SandiaTable * fGasPhotoAbsCof
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4double GetAngleXTR(G4int iTR, G4double position, G4int iAngle)
const G4String & GetProcessName() const
G4double GetUserElapsed() const
G4VXTRenergyLoss(G4LogicalVolume *anEnvelope, G4Material *, G4Material *, G4double, G4double, G4int, const G4String &processName="XTRenergyLoss", G4ProcessType type=fElectromagnetic)
static constexpr double eV
G4double GetGlobalTime() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep) override
G4StepPoint * GetPostStepPoint() const
void ProposeEnergy(G4double finalEnergy)
G4double GetPlateLinearPhotoAbs(G4double)
G4double GetPlateFormationZone(G4double, G4double, G4double)
G4double GetSandiaCofForMaterial(G4int, G4int) const
void ComputePlatePhotoAbsCof()
G4PhysicsLogVector * fProtonEnergyVector
ThreeVector shoot(const G4int Ap, const G4int Af)
std::vector< G4PhysicsTable * > fAngleBank
G4double GetXTRrandomEnergy(G4double scaledTkin, G4int iTkin)
G4double G4ParticleHPJENDLHEData::G4double result
G4complex GetPlateComplexFZ(G4double, G4double, G4double)
std::vector< G4Material * > G4MaterialTable
G4double XTRNAngleSpectralDensity(G4double energy)
G4double GetComptonPerAtom(G4double, G4double)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double GetRandomAngle(G4double energyXTR, G4int iTkin)
void GetNumberOfPhotons()
G4PhysicsTable * fAngleDistrTable
static constexpr double c_light
void SetParentID(const G4int aValue)
static constexpr double barn
virtual G4double GetStackFactor(G4double energy, G4double gamma, G4double varAngle)
G4PhysicsFreeVector * GetAngleVector(G4double energy, G4int n)
G4double GetPlateCompton(G4double)
virtual G4double GetMeanFreePath(const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition) override
void SetProcessSubType(G4int)
G4double GetKineticEnergy() const
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4VSolid * GetSolid() const
G4SandiaTable * GetSandiaTable() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &) override
void PutValue(size_t index, G4double energy, G4double dataValue)
static constexpr double pi
G4double XTRNSpectralDensity(G4double energy)
static constexpr double fine_structure_const
Float_t x2[n_points_geant4]
G4double GetGasLinearPhotoAbs(G4double)
G4PhysicsTable * fAngleForEnergyTable
G4double AngleSpectralXTRdEdx(G4double energy)
G4double XTRNSpectralAngleDensity(G4double varAngle)
G4double Legendre10(T &typeT, F f, G4double a, G4double b)
static constexpr double GeV
const G4DynamicParticle * GetDynamicParticle() const
G4double GetElectronDensity() const
G4LogicalVolume * fEnvelope
void PutValue(size_t index, G4double theValue)
virtual ~G4VXTRenergyLoss()
G4VPhysicalVolume * GetVolume() const
G4complex OneInterfaceXTRdEdx(G4double energy, G4double gamma, G4double varAngle)
G4double XTRNAngleDensity(G4double varAngle)
void SetNumberOfSecondaries(G4int totSecondaries)
G4SandiaTable * fPlatePhotoAbsCof