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