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 "G4RegularXTRadiator.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "Randomize.hh"
00035
00036 #include "G4Gamma.hh"
00037
00039
00040
00041
00042 G4RegularXTRadiator::G4RegularXTRadiator(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 G4cout<<"Regular X-ray TR radiator EM process is called"<<G4endl ;
00049
00050
00051
00052
00053 fAlphaPlate = 10000;
00054 fAlphaGas = 1000;
00055 G4cout<<"fAlphaPlate = "<<fAlphaPlate<<" ; fAlphaGas = "<<fAlphaGas<<G4endl ;
00056
00057
00058 }
00059
00061
00062 G4RegularXTRadiator::~G4RegularXTRadiator()
00063 {
00064 ;
00065 }
00066
00068
00069
00070
00071 G4double G4RegularXTRadiator::SpectralXTRdEdx(G4double energy)
00072 {
00073 G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC, theta2, theta2k;
00074 G4double aMa, bMb ,sigma, dump;
00075 G4int k, kMax, kMin;
00076
00077 aMa = fPlateThick*GetPlateLinearPhotoAbs(energy);
00078 bMb = fGasThick*GetGasLinearPhotoAbs(energy);
00079 sigma = 0.5*(aMa + bMb);
00080 dump = std::exp(-fPlateNumber*sigma);
00081 if(verboseLevel > 2) G4cout<<" dump = "<<dump<<G4endl;
00082 cofPHC = 4*pi*hbarc;
00083 tmp = (fSigma1 - fSigma2)/cofPHC/energy;
00084 cof1 = fPlateThick*tmp;
00085 cof2 = fGasThick*tmp;
00086
00087 cofMin = energy*(fPlateThick + fGasThick)/fGamma/fGamma;
00088 cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
00089 cofMin /= cofPHC;
00090
00091 theta2 = cofPHC/(energy*(fPlateThick + fGasThick));
00092
00093
00094
00095
00096
00097 kMin = G4int(cofMin);
00098 if (cofMin > kMin) kMin++;
00099
00100
00101
00102
00103
00104
00105
00106
00107 kMax = kMin + 49;
00108
00109
00110
00111
00112 if(verboseLevel > 2)
00113 {
00114 G4cout<<cof1<<" "<<cof2<<" "<<cofMin<<G4endl;
00115 G4cout<<"kMin = "<<kMin<<"; kMax = "<<kMax<<G4endl;
00116 }
00117 for( k = kMin; k <= kMax; k++ )
00118 {
00119 tmp = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
00120 result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
00121
00122 if( k == kMin && kMin == G4int(cofMin) )
00123 {
00124 sum += 0.5*std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
00125 }
00126 else
00127 {
00128 sum += std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
00129 }
00130 theta2k = std::sqrt(theta2*std::abs(k-cofMin));
00131
00132 if(verboseLevel > 2)
00133 {
00134
00135
00136 G4cout<<k<<" "<<theta2k<<" "<<std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result
00137 <<" "<<sum<<G4endl;
00138 }
00139 }
00140 result = 2*( cof1 + cof2 )*( cof1 + cof2 )*sum/energy;
00141
00142
00143 result *= ( 1 - dump + 2*dump*fPlateNumber );
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157 return result;
00158 }
00159
00160
00161
00163
00164
00165
00166
00167
00168
00169
00170 G4double
00171 G4RegularXTRadiator::GetStackFactor( G4double energy,
00172 G4double gamma, G4double varAngle )
00173 {
00174
00175
00176
00177 G4double result, Za, Zb, Ma, Mb;
00178
00179 Za = GetPlateFormationZone(energy,gamma,varAngle);
00180 Zb = GetGasFormationZone(energy,gamma,varAngle);
00181
00182 Ma = GetPlateLinearPhotoAbs(energy);
00183 Mb = GetGasLinearPhotoAbs(energy);
00184
00185
00186 G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPlate,fPlateThick/Za/fAlphaPlate);
00187 G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas,fGasThick/Zb/fAlphaGas);
00188
00189 G4complex Ha = std::pow(Ca,-fAlphaPlate);
00190 G4complex Hb = std::pow(Cb,-fAlphaGas);
00191 G4complex H = Ha*Hb;
00192
00193 G4complex F1 = (1.0 - Ha)*(1.0 - Hb )/(1.0 - H)
00194 * G4double(fPlateNumber);
00195
00196 G4complex F2 = (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H)/(1.0-H)
00197 * (1.0 - std::pow(H,fPlateNumber));
00198
00199 G4complex R = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
00200
00201 result = 2.0*std::real(R);
00202
00203 return result;
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 }
00246
00247
00248
00249
00251
00252
00253
00254
00255
00256
00257
00258