| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // 19.04.24 V. Grichine, first version
27 //
28
29 #include "G4XTRGaussRadModel.hh"
30 #include "G4PhysicalConstants.hh"
31
32 #include "G4Material.hh"
33 #include "G4Element.hh"
34 #include "G4NistManager.hh"
35
36 using namespace std;
37 using namespace CLHEP;
38
39 ////////////////////////////////////////////////////////////////////////////
40 // Constructor, destructor
41
42 G4XTRGaussRadModel::G4XTRGaussRadModel(
43 G4LogicalVolume* anEnvelope, G4double alphaPlate, G4double alphaGas, G4Material* foilMat, G4Material* gasMat,
44 G4double aa, G4double b, G4int n, const G4String& processName)
45 : G4VXTRenergyLoss(anEnvelope, foilMat, gasMat, aa, b, n, processName)
46 {
47 if( verboseLevel > 0 )
48 G4cout << "G4XTRGaussRadModel EM process is called"
49 << G4endl;
50
51 fAlphaPlate = alphaPlate; // ~40
52 fAlphaGas = alphaGas; // ~10
53 fExitFlux = true; // XTR photons are moved to the end of radiator
54 }
55
56 ///////////////////////////////////////////////////////////////////////////
57
58 G4XTRGaussRadModel::~G4XTRGaussRadModel() = default;
59
60 ///////////////////////////////////////////////////////////////////////////
61
62 void G4XTRGaussRadModel::ProcessDescription(std::ostream& out) const
63 {
64 out << "Simulation of forward X-ray transition radiation generated by\n"
65 "relativistic charged particles crossing the interface between\n"
66 "two materials.\n";
67 }
68
69 ///////////////////////////////////////////////////////////////////////////
70
71 G4double G4XTRGaussRadModel::SpectralXTRdEdx(G4double energy)
72 {
73 static constexpr G4double cofPHC = 4. * pi * hbarc;
74 G4double result, sum = 0., tmp, cof1, cof2, cofMin, theta2, theta2k;
75 G4double aMa, bMb, sigma, dump;
76 G4int k, kMax, kMin;
77
78 aMa = fPlateThick * GetPlateLinearPhotoAbs(energy);
79 bMb = fGasThick * GetGasLinearPhotoAbs(energy);
80 sigma = 0.5 * (aMa + bMb);
81 dump = std::exp(-fPlateNumber * sigma);
82 if(verboseLevel > 2)
83 G4cout << " dump = " << dump << G4endl;
84 tmp = (fSigma1 - fSigma2) / cofPHC / energy;
85 cof1 = fPlateThick * tmp;
86 cof2 = fGasThick * tmp;
87
88 cofMin = energy * (fPlateThick + fGasThick) / fGamma / fGamma;
89 cofMin += (fPlateThick * fSigma1 + fGasThick * fSigma2) / energy;
90 cofMin /= cofPHC;
91
92 theta2 = cofPHC / (energy * (fPlateThick + fGasThick));
93
94 kMin = G4int(cofMin);
95 if(cofMin > kMin)
96 kMin++;
97
98 kMax = kMin + 200; // 99; // 49; //
99
100 if(verboseLevel > 2)
101 {
102 G4cout << cof1 << " " << cof2 << " " << cofMin << G4endl;
103 G4cout << "kMin = " << kMin << "; kMax = " << kMax << G4endl;
104 }
105 for(k = kMin; k <= kMax; ++k)
106 {
107 tmp = pi * fPlateThick * (k + cof2) / (fPlateThick + fGasThick);
108 result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2);
109 if(k == kMin && kMin == G4int(cofMin))
110 {
111 sum +=
112 0.5 * std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result;
113 }
114 else
115 {
116 sum += std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result;
117 }
118 theta2k = std::sqrt(theta2 * std::abs(k - cofMin));
119
120 if(verboseLevel > 2)
121 {
122 G4cout << k << " " << theta2k << " "
123 << std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result
124 << " " << sum << G4endl;
125 }
126 }
127 result = 2 * (cof1 + cof2) * (cof1 + cof2) * sum / energy;
128 result *= dump * (-1 + dump + 2 * fPlateNumber);
129 return result;
130 }
131
132 ///////////////////////////////////////////////////////////////////////////
133 //
134 // Approximation for radiator interference factor for the case of
135 // Gauss-distributed regular radiator. The plate and gas gap thicknesses
136 // are Gauss distributed with RMS
137 // sa and sb for plate and gas, respectively.
138 // The mean values of the plate and gas gap thicknesses
139 // are supposed to be about XTR formation zones.
140 // The XTR photons are moved to the end of radiator
141
142
143 G4double G4XTRGaussRadModel::GetStackFactor(G4double energy,
144 G4double gamma,
145 G4double varAngle)
146 {
147 G4double result(0.);
148
149 G4double Ma, Mb, aMa, bMb, sigma;
150
151 G4double sa = fPlateThick/fAlphaPlate;
152 G4double sb = fGasThick/fAlphaGas;
153
154 Ma = GetPlateLinearPhotoAbs(energy);
155 aMa = fPlateThick * Ma;
156 Mb = GetGasLinearPhotoAbs(energy);
157 bMb = fGasThick * Mb;
158 sigma = aMa + bMb; // m1*t1+m2*t2 dimensionless
159 G4double nn = G4double( fPlateNumber );
160
161 // Gauss fluctuation of foil and gas gaps according to
162 // RMS = sa = a/fAlphaPlate and RMS = sb = b/fAlphaGas
163
164 G4complex med(0.,1.);
165 G4complex Z1 = GetPlateComplexFZ( energy, gamma, varAngle);
166 G4complex por1 = -0.5*med*fPlateThick/Z1 - 0.125*sa*sa/Z1/Z1;
167
168 G4complex Z2 = GetGasComplexFZ( energy, gamma, varAngle);
169 G4complex por2 = -0.5*med*fGasThick/Z2 - 0.125*sb*sb/Z2/Z2;
170
171 G4complex npor = (por1+por2)*nn;
172
173 G4complex Ha = exp(por1);
174 G4complex Hb = exp(por2);
175
176 G4complex H = Ha*Hb;
177 G4complex Hn = exp(npor);
178
179 G4complex A = ( 1.0 - Ha ) * ( 1.0 - Hb ) / ( 1. - H );
180 G4complex B1 = ( 1.0 - Ha ) * ( 1.0 - Ha ) * Hb / ( 1. - H );
181
182 G4complex R = B1*( exp(-sigma*nn) - Hn )/( exp(-sigma) - H );
183
184 R += A * ( 1 - exp(-sigma*nn) ) / ( 1. - exp(-sigma) ); // -> A*nn
185
186 R *= OneInterfaceXTRdEdx(energy, gamma, varAngle);
187
188 result = 2.0 * std::real(R);
189
190 return result;
191 }
192