00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include <complex>
00031
00032 #include "G4TransparentRegXTRadiator.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "Randomize.hh"
00035 #include "G4Integrator.hh"
00036 #include "G4Gamma.hh"
00037
00039
00040
00041
00042 G4TransparentRegXTRadiator::G4TransparentRegXTRadiator(G4LogicalVolume *anEnvelope,
00043 G4Material* foilMat,G4Material* gasMat,
00044 G4double a, G4double b, G4int n,
00045 const G4String& processName) :
00046 G4VXTRenergyLoss(anEnvelope,foilMat,gasMat,a,b,n,processName)
00047 {
00048 if(verboseLevel > 0)
00049 G4cout<<"Regular transparent X-ray TR radiator EM process is called"<<G4endl;
00050
00051
00052
00053
00054 fAlphaPlate = 10000;
00055 fAlphaGas = 1000;
00056
00057
00058 }
00059
00061
00062 G4TransparentRegXTRadiator::~G4TransparentRegXTRadiator()
00063 {
00064 ;
00065 }
00066
00068
00069
00070
00071 G4double G4TransparentRegXTRadiator::SpectralXTRdEdx(G4double energy)
00072 {
00073 G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC, theta2, theta2k ;
00074 G4int k, kMax, kMin;
00075
00076
00077
00078
00079
00080 cofPHC = 4*pi*hbarc;
00081 tmp = (fSigma1 - fSigma2)/cofPHC/energy;
00082 cof1 = fPlateThick*tmp;
00083 cof2 = fGasThick*tmp;
00084
00085 cofMin = energy*(fPlateThick + fGasThick)/fGamma/fGamma;
00086 cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
00087 cofMin /= cofPHC;
00088
00089 theta2 = cofPHC/(energy*(fPlateThick + fGasThick));
00090
00091
00092
00093
00094
00095 kMin = G4int(cofMin);
00096 if (cofMin > kMin) kMin++;
00097
00098
00099
00100
00101
00102
00103
00104
00105 kMax = kMin + 49;
00106
00107
00108
00109
00110 if(verboseLevel > 2)
00111 {
00112 G4cout<<cof1<<" "<<cof2<<" "<<cofMin<<G4endl;
00113 G4cout<<"kMin = "<<kMin<<"; kMax = "<<kMax<<G4endl;
00114 }
00115 for( k = kMin; k <= kMax; k++ )
00116 {
00117 tmp = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
00118 result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
00119
00120 if( k == kMin && kMin == G4int(cofMin) )
00121 {
00122 sum += 0.5*std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
00123 }
00124 else
00125 {
00126 sum += std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
00127 }
00128 theta2k = std::sqrt(theta2*std::abs(k-cofMin));
00129
00130 if(verboseLevel > 2)
00131 {
00132
00133
00134 G4cout<<k<<" "<<theta2k<<" "<<std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result
00135 <<" "<<sum<<G4endl;
00136 }
00137 }
00138 result = 4*( cof1 + cof2 )*( cof1 + cof2 )*sum/energy;
00139
00140
00141 result *= fPlateNumber;
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156 return result;
00157 }
00158
00159
00161
00162
00163
00164
00165
00166
00167
00168 G4double
00169 G4TransparentRegXTRadiator::GetStackFactor( G4double energy,
00170 G4double gamma, G4double varAngle )
00171 {
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, D, sigma;
00200
00201 aZa = fPlateThick/GetPlateFormationZone(energy,gamma,varAngle);
00202 bZb = fGasThick/GetGasFormationZone(energy,gamma,varAngle);
00203 aMa = fPlateThick*GetPlateLinearPhotoAbs(energy);
00204 bMb = fGasThick*GetGasLinearPhotoAbs(energy);
00205 sigma = aMa*fPlateThick + bMb*fGasThick;
00206 Qa = std::exp(-0.5*aMa);
00207 Qb = std::exp(-0.5*bMb);
00208 Q = Qa*Qb;
00209
00210 G4complex Ha( Qa*std::cos(aZa), -Qa*std::sin(aZa) );
00211 G4complex Hb( Qb*std::cos(bZb), -Qb*std::sin(bZb) );
00212 G4complex H = Ha*Hb;
00213 G4complex Hs = conj(H);
00214 D = 1.0 /( (1 - Q)*(1 - Q) +
00215 4*Q*std::sin(0.5*(aZa + bZb))*std::sin(0.5*(aZa + bZb)) );
00216 G4complex F1 = (1.0 - Ha)*(1.0 - Hb)*(1.0 - Hs)
00217 * G4double(fPlateNumber)*D;
00218 G4complex F2 = (1.0 - Ha)*(1.0 - Ha)*Hb*(1.0 - Hs)*(1.0 - Hs)
00219
00220 * (1.0 - std::exp(-0.5*fPlateNumber*sigma)) * D*D;
00221 G4complex R = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
00222 result = 2.0*std::real(R);
00223 return result;
00224
00225 }
00226
00227
00228
00229
00231
00232
00233
00234
00235
00236
00237
00238