| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 //
>> 27 // $Id: G4StrawTubeXTRadiator.cc,v 1.4 2006/06/29 19:56:13 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01-patch-01 $
>> 29 //
>> 30
>> 31 #include <complex>
26 32
27 #include "G4StrawTubeXTRadiator.hh" 33 #include "G4StrawTubeXTRadiator.hh"
>> 34 #include "Randomize.hh"
28 35
29 #include "G4Gamma.hh" 36 #include "G4Gamma.hh"
30 #include "G4PhysicalConstants.hh" << 37
31 #include "G4SystemOfUnits.hh" << 38 using namespace std;
32 39
33 ////////////////////////////////////////////// 40 ////////////////////////////////////////////////////////////////////////////
>> 41 //
34 // Constructor, destructor 42 // Constructor, destructor
35 G4StrawTubeXTRadiator::G4StrawTubeXTRadiator(G << 43
36 G << 44 G4StrawTubeXTRadiator::G4StrawTubeXTRadiator(G4LogicalVolume *anEnvelope,
37 G << 45 G4Material* foilMat,G4Material* gasMat,
38 G << 46 G4double a, G4double b, G4Material* mediumMat,
39 G << 47 G4bool unishut,
40 c << 48 const G4String& processName) :
41 : G4VXTRenergyLoss(anEnvelope, foilMat, gasM << 49 G4VXTRenergyLoss(anEnvelope,foilMat,gasMat,a,b,1,processName)
42 { 50 {
43 if(verboseLevel > 0) << 51 G4cout<<"Straw tube X-ray TR radiator EM process is called"<<G4endl;
44 G4cout << "Straw tube X-ray TR radiator E <<
45 52
46 if(unishut) << 53 if( unishut )
47 { 54 {
48 fAlphaPlate = 1. / 3.; << 55 fAlphaPlate = 1./3.;
49 fAlphaGas = 12.4; 56 fAlphaGas = 12.4;
50 if(verboseLevel > 0) << 57 G4cout<<"straw uniform shooting: "<<"fAlphaPlate = "
51 G4cout << "straw uniform shooting: " << 58 <<fAlphaPlate<<" ; fAlphaGas = "<<fAlphaGas<<G4endl;
52 << "fAlphaPlate = " << fAlphaPlat << 59
53 << " ; fAlphaGas = " << fAlphaGas <<
54 } 60 }
55 else 61 else
56 { 62 {
57 fAlphaPlate = 0.5; 63 fAlphaPlate = 0.5;
58 fAlphaGas = 5.; 64 fAlphaGas = 5.;
59 if(verboseLevel > 0) << 65 G4cout<<"straw isotropical shooting: "<<"fAlphaPlate = "
60 G4cout << "straw isotropical shooting: " << 66 <<fAlphaPlate<<" ; fAlphaGas = "<<fAlphaGas<<G4endl;
61 << "fAlphaPlate = " << fAlphaPlat << 67
62 << " ; fAlphaGas = " << fAlphaGas <<
63 } <<
64 68
>> 69 }
65 // index of medium material 70 // index of medium material
66 fMatIndex3 = (G4int)mediumMat->GetIndex(); <<
67 if(verboseLevel > 0) <<
68 G4cout << "medium material = " << mediumMa <<
69 <<
70 // plasma energy squared for plate material <<
71 fSigma3 = fPlasmaCof * mediumMat->GetElectro <<
72 if(verboseLevel > 0) <<
73 G4cout << "medium plasma energy = " << std <<
74 << G4endl; <<
75 71
76 // Compute cofs for preparation of linear ph << 72 fMatIndex3 = mediumMat->GetIndex();
>> 73 G4cout<<"medium material = "<<mediumMat->GetName()<<G4endl;
>> 74
>> 75 // plasma energy squared for plate material
>> 76
>> 77 fSigma3 = fPlasmaCof*mediumMat->GetElectronDensity();
>> 78 G4cout<<"medium plasma energy = "<<sqrt(fSigma3)/eV<<" eV"<<G4endl;
>> 79
>> 80 // Compute cofs for preparation of linear photo absorption in external medium
>> 81
77 ComputeMediumPhotoAbsCof(); 82 ComputeMediumPhotoAbsCof();
>> 83
>> 84
>> 85
>> 86
>> 87 // Build energy and angular integral spectra of X-ray TR photons from
>> 88 // a radiator
>> 89
>> 90 // BuildTable();
78 } 91 }
79 92
80 ////////////////////////////////////////////// 93 ///////////////////////////////////////////////////////////////////////////
81 G4StrawTubeXTRadiator::~G4StrawTubeXTRadiator( <<
82 94
83 void G4StrawTubeXTRadiator::ProcessDescription << 95 G4StrawTubeXTRadiator::~G4StrawTubeXTRadiator()
84 { 96 {
85 out << "Simulation of forward X-ray transiti << 97 ;
86 "a straw tube radiator.\n"; <<
87 } 98 }
88 99
>> 100
>> 101
89 ////////////////////////////////////////////// 102 ///////////////////////////////////////////////////////////////////////////
>> 103 //
90 // Approximation for radiator interference fac 104 // Approximation for radiator interference factor for the case of
91 // straw tube radiator. The plate (window, str << 105 // straw tube radiator. The plate (window, straw wall) and gas (inside straw)
92 // gap thicknesses are gamma distributed. << 106 // gap thicknesses are gamma distributed.
93 // The mean values of the plate and gas gap th << 107 // The mean values of the plate and gas gap thicknesses
94 // are supposed to be about XTR formation zone 108 // are supposed to be about XTR formation zone.
95 G4double G4StrawTubeXTRadiator::GetStackFactor << 109
96 << 110 G4double
>> 111 G4StrawTubeXTRadiator::GetStackFactor( G4double energy,
>> 112 G4double gamma, G4double varAngle )
97 { 113 {
98 G4double result, L2, L3, M2, M3; <<
99 114
100 L2 = GetPlateFormationZone(energy, gamma, va << 115
101 L3 = GetGasFormationZone(energy, gamma, varA << 116 G4double result, L2, L3, M2, M3;
>> 117
>> 118 L2 = GetPlateFormationZone(energy,gamma,varAngle);
>> 119 L3 = GetGasFormationZone(energy,gamma,varAngle);
102 120
103 M2 = GetPlateLinearPhotoAbs(energy); 121 M2 = GetPlateLinearPhotoAbs(energy);
104 M3 = GetGasLinearPhotoAbs(energy); 122 M3 = GetGasLinearPhotoAbs(energy);
105 123
106 G4complex C2(1.0 + 0.5 * fPlateThick * M2 / <<
107 fPlateThick / L2 / fAlphaPlate) <<
108 G4complex C3(1.0 + 0.5 * fGasThick * M3 / fA <<
109 fGasThick / L3 / fAlphaGas); <<
110 <<
111 G4complex H2 = std::pow(C2, -fAlphaPlate); <<
112 G4complex H3 = std::pow(C3, -fAlphaGas); <<
113 G4complex H = H2 * H3; <<
114 <<
115 G4complex Z1 = GetMediumComplexFZ(energy, ga <<
116 G4complex Z2 = GetPlateComplexFZ(energy, gam <<
117 G4complex Z3 = GetGasComplexFZ(energy, gamma <<
118 <<
119 G4complex R = (Z1 - Z2) * (Z1 - Z2) * (1. - <<
120 (Z2 - Z3) * (Z2 - Z3) * (1. - <<
121 2. * (Z1 - Z2) * (Z2 - Z3) * H <<
122 124
123 result = 2.0 * std::real(R) * (varAngle * en << 125 G4complex C2(1.0 + 0.5*fPlateThick*M2/fAlphaPlate, fPlateThick/L2/fAlphaPlate);
>> 126 G4complex C3(1.0 + 0.5*fGasThick*M3/fAlphaGas, fGasThick/L3/fAlphaGas);
>> 127
>> 128 G4complex H2 = pow(C2,-fAlphaPlate);
>> 129 G4complex H3 = pow(C3,-fAlphaGas);
>> 130 G4complex H = H2*H3;
>> 131
>> 132 G4complex Z1 = GetMediumComplexFZ(energy,gamma,varAngle);
>> 133 G4complex Z2 = GetPlateComplexFZ(energy,gamma,varAngle);
>> 134 G4complex Z3 = GetGasComplexFZ(energy,gamma,varAngle);
>> 135
>> 136
>> 137 G4complex R = ( Z1 - Z2 )*( Z1 - Z2 )*( 1. - H2*H ) +
>> 138 ( Z2 - Z3 )*( Z2 - Z3 )*( 1. - H3 ) +
>> 139 2.*( Z1 - Z2 )*( Z2 - Z3 )*H2*( 1. - H3 ) ;
>> 140
>> 141 result = 2.0*real(R)*(varAngle*energy/hbarc/hbarc);
>> 142
>> 143 return result;
124 144
125 return result; <<
126 } 145 }
127 146
128 ////////////////////////////////////////////// << 147
>> 148 //////////////////////////////////////////////////////////////////////
>> 149 //////////////////////////////////////////////////////////////////////
>> 150 //////////////////////////////////////////////////////////////////////
>> 151 //
129 // Calculates formation zone for external medi 152 // Calculates formation zone for external medium. Omega is energy !!!
130 G4double G4StrawTubeXTRadiator::GetMediumForma << 153
131 << 154 G4double G4StrawTubeXTRadiator::GetMediumFormationZone( G4double omega ,
132 << 155 G4double gamma ,
>> 156 G4double varAngle )
133 { 157 {
134 G4double cof, lambda; 158 G4double cof, lambda;
135 lambda = 1.0 / gamma / gamma + varAngle + fS << 159 lambda = 1.0/gamma/gamma + varAngle + fSigma3/omega/omega;
136 cof = 2.0 * hbarc / omega / lambda; << 160 cof = 2.0*hbarc/omega/lambda ;
137 return cof; << 161 return cof ;
138 } 162 }
139 163
140 ////////////////////////////////////////////// << 164 //////////////////////////////////////////////////////////////////////
>> 165 //
141 // Calculates complex formation zone for exter 166 // Calculates complex formation zone for external medium. Omega is energy !!!
142 G4complex G4StrawTubeXTRadiator::GetMediumComp << 167
143 << 168 G4complex G4StrawTubeXTRadiator::GetMediumComplexFZ( G4double omega ,
144 << 169 G4double gamma ,
>> 170 G4double varAngle )
145 { 171 {
146 G4double cof, length, delta, real_v, image_v << 172 G4double cof, length,delta, real, image;
147 173
148 length = 0.5 * GetMediumFormationZone(omega, << 174 length = 0.5*GetMediumFormationZone(omega,gamma,varAngle);
149 delta = length * GetMediumLinearPhotoAbs(om << 175 delta = length*GetMediumLinearPhotoAbs(omega);
150 cof = 1.0 / (1.0 + delta * delta); << 176 cof = 1.0/(1.0 + delta*delta);
151 177
152 real_v = length * cof; << 178 real = length*cof;
153 image_v = real_v * delta; << 179 image = real*delta;
154 180
155 G4complex zone(real_v, image_v); << 181 G4complex zone(real,image);
156 return zone; 182 return zone;
157 } 183 }
158 184
159 ////////////////////////////////////////////// 185 ////////////////////////////////////////////////////////////////////////
>> 186 //
160 // Computes matrix of Sandia photo absorption 187 // Computes matrix of Sandia photo absorption cross section coefficients for
161 // medium material 188 // medium material
162 void G4StrawTubeXTRadiator::ComputeMediumPhoto << 189
>> 190 void G4StrawTubeXTRadiator::ComputeMediumPhotoAbsCof()
163 { 191 {
164 const G4MaterialTable* theMaterialTable = G4 << 192 G4int i, j, numberOfElements;
165 const G4Material* mat = (* << 193 static const G4MaterialTable*
166 fMediumPhotoAbsCof = ma << 194 theMaterialTable = G4Material::GetMaterialTable();
>> 195
>> 196 G4SandiaTable thisMaterialSandiaTable(fMatIndex3);
>> 197 numberOfElements = (*theMaterialTable)[fMatIndex3]->GetNumberOfElements();
>> 198 G4int* thisMaterialZ = new G4int[numberOfElements];
>> 199
>> 200 for(i=0;i<numberOfElements;i++)
>> 201 {
>> 202 thisMaterialZ[i] = (G4int)(*theMaterialTable)[fMatIndex3]->
>> 203 GetElement(i)->GetZ() ;
>> 204 }
>> 205 fMediumIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
>> 206 (thisMaterialZ,numberOfElements) ;
>> 207
>> 208 fMediumIntervalNumber = thisMaterialSandiaTable.SandiaMixing
>> 209 ( thisMaterialZ ,
>> 210 (*theMaterialTable)[fMatIndex3]->GetFractionVector() ,
>> 211 numberOfElements,fMediumIntervalNumber);
>> 212
>> 213 fMediumPhotoAbsCof = new G4double*[fMediumIntervalNumber];
>> 214
>> 215 for(i=0;i<fMediumIntervalNumber;i++)
>> 216 {
>> 217 fMediumPhotoAbsCof[i] = new G4double[5];
>> 218 }
>> 219 for(i=0;i<fMediumIntervalNumber;i++)
>> 220 {
>> 221 fMediumPhotoAbsCof[i][0] = thisMaterialSandiaTable.
>> 222 GetPhotoAbsorpCof(i+1,0);
>> 223
>> 224 for(j=1;j<5;j++)
>> 225 {
>> 226 fMediumPhotoAbsCof[i][j] = thisMaterialSandiaTable.
>> 227 GetPhotoAbsorpCof(i+1,j)*
>> 228 (*theMaterialTable)[fMatIndex3]->GetDensity();
>> 229 }
>> 230 }
>> 231 delete[] thisMaterialZ;
>> 232 return;
167 } 233 }
168 234
169 ////////////////////////////////////////////// 235 //////////////////////////////////////////////////////////////////////
170 // Returns the value of linear photo absorptio << 236 //
>> 237 // Returns the value of linear photo absorption coefficient (in reciprocal
171 // length) for medium for given energy of X-ra 238 // length) for medium for given energy of X-ray photon omega
172 G4double G4StrawTubeXTRadiator::GetMediumLinea <<
173 { <<
174 G4double omega2, omega3, omega4; <<
175 239
176 omega2 = omega * omega; << 240 G4double G4StrawTubeXTRadiator::GetMediumLinearPhotoAbs(G4double omega)
177 omega3 = omega2 * omega; << 241 {
178 omega4 = omega2 * omega2; << 242 G4int i ;
>> 243 G4double omega2, omega3, omega4;
179 244
180 const G4double* SandiaCof = << 245 omega2 = omega*omega;
181 fMediumPhotoAbsCof->GetSandiaCofForMateria << 246 omega3 = omega2*omega;
>> 247 omega4 = omega2*omega2;
182 248
183 G4double cross = SandiaCof[0] / omega + Sand << 249 for(i=0;i<fMediumIntervalNumber;i++)
184 SandiaCof[2] / omega3 + San << 250 {
185 return cross; << 251 if( omega < fMediumPhotoAbsCof[i][0] ) break;
>> 252 }
>> 253 if( i == 0 )
>> 254 {
>> 255 G4Exception("Invalid (<I1) energy in G4VXTRenergyLoss::GetMediumLinearPhotoAbs");
>> 256 }
>> 257 else i-- ;
>> 258
>> 259 return fMediumPhotoAbsCof[i][1]/omega + fMediumPhotoAbsCof[i][2]/omega2 +
>> 260 fMediumPhotoAbsCof[i][3]/omega3 + fMediumPhotoAbsCof[i][4]/omega4 ;
186 } 261 }
>> 262
>> 263
>> 264
>> 265
>> 266 //
>> 267 //
>> 268 ////////////////////////////////////////////////////////////////////////////
>> 269
>> 270
>> 271
>> 272
>> 273
>> 274
>> 275
>> 276
187 277