| 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 // History: 26 // History:
27 // 1st version 11.09.97 V. Grichine (Vladimir. 27 // 1st version 11.09.97 V. Grichine (Vladimir.Grichine@cern.ch )
28 // 2nd version 17.12.97 V. Grichine 28 // 2nd version 17.12.97 V. Grichine
29 // 17-09-01, migration of Materials to pure ST 29 // 17-09-01, migration of Materials to pure STL (mma)
30 // 10-03-03, migration to "cut per region" (V. 30 // 10-03-03, migration to "cut per region" (V.Ivanchenko)
31 // 03.06.03, V.Ivanchenko fix compilation warn 31 // 03.06.03, V.Ivanchenko fix compilation warnings
32 32
33 #include "G4ForwardXrayTR.hh" 33 #include "G4ForwardXrayTR.hh"
34 34
35 #include "globals.hh" 35 #include "globals.hh"
36 #include "G4Gamma.hh" 36 #include "G4Gamma.hh"
37 #include "G4GeometryTolerance.hh" 37 #include "G4GeometryTolerance.hh"
38 #include "G4Material.hh" 38 #include "G4Material.hh"
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4PhysicsLinearVector.hh" 40 #include "G4PhysicsLinearVector.hh"
41 #include "G4PhysicsLogVector.hh" 41 #include "G4PhysicsLogVector.hh"
42 #include "G4PhysicsTable.hh" 42 #include "G4PhysicsTable.hh"
43 #include "G4PhysicsVector.hh" 43 #include "G4PhysicsVector.hh"
44 #include "G4Poisson.hh" 44 #include "G4Poisson.hh"
45 #include "G4ProductionCutsTable.hh" 45 #include "G4ProductionCutsTable.hh"
46 #include "G4SystemOfUnits.hh" 46 #include "G4SystemOfUnits.hh"
47 #include "G4PhysicsModelCatalog.hh" 47 #include "G4PhysicsModelCatalog.hh"
48 48
49 ////////////////////////////////////////////// 49 //////////////////////////////////////////////////////////////////////
50 // 50 //
51 // Constructor for creation of physics tables 51 // Constructor for creation of physics tables (angle and energy TR
52 // distributions) for a couple of selected mat 52 // distributions) for a couple of selected materials.
53 // 53 //
54 // Recommended for use in applications with ma 54 // Recommended for use in applications with many materials involved,
55 // when only few (usually couple) materials ar 55 // when only few (usually couple) materials are interested for generation
56 // of TR on the interface between them 56 // of TR on the interface between them
57 G4ForwardXrayTR::G4ForwardXrayTR(const G4Strin 57 G4ForwardXrayTR::G4ForwardXrayTR(const G4String& matName1,
58 const G4Strin 58 const G4String& matName2,
59 const G4Strin 59 const G4String& processName)
60 : G4TransitionRadiation(processName) 60 : G4TransitionRadiation(processName)
61 { 61 {
62 secID = G4PhysicsModelCatalog::GetModelID("m 62 secID = G4PhysicsModelCatalog::GetModelID("model_XrayTR");
63 fPtrGamma = nullptr; 63 fPtrGamma = nullptr;
64 fGammaCutInKineticEnergy = nullptr; 64 fGammaCutInKineticEnergy = nullptr;
65 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR 65 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = 0.0;
66 fGamma = fSigma1 = fSigma2 = 0.0; 66 fGamma = fSigma1 = fSigma2 = 0.0;
67 fAngleDistrTable = nullptr; 67 fAngleDistrTable = nullptr;
68 fEnergyDistrTable = nullptr; 68 fEnergyDistrTable = nullptr;
69 fMatIndex1 = fMatIndex2 = 0; 69 fMatIndex1 = fMatIndex2 = 0;
70 70
71 // Proton energy vector initialization 71 // Proton energy vector initialization
72 fProtonEnergyVector = 72 fProtonEnergyVector =
73 new G4PhysicsLogVector(fMinProtonTkin, fMa 73 new G4PhysicsLogVector(fMinProtonTkin, fMaxProtonTkin, fTotBin);
74 G4int iMat; 74 G4int iMat;
75 const G4ProductionCutsTable* theCoupleTable 75 const G4ProductionCutsTable* theCoupleTable =
76 G4ProductionCutsTable::GetProductionCutsTa 76 G4ProductionCutsTable::GetProductionCutsTable();
77 G4int numOfCouples = (G4int)theCoupleTable-> 77 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
78 78
79 G4bool build = true; 79 G4bool build = true;
80 80
81 for(iMat = 0; iMat < numOfCouples; ++iMat) 81 for(iMat = 0; iMat < numOfCouples; ++iMat) // check first material name
82 { 82 {
83 const G4MaterialCutsCouple* couple = 83 const G4MaterialCutsCouple* couple =
84 theCoupleTable->GetMaterialCutsCouple(iM 84 theCoupleTable->GetMaterialCutsCouple(iMat);
85 if(matName1 == couple->GetMaterial()->GetN 85 if(matName1 == couple->GetMaterial()->GetName())
86 { 86 {
87 fMatIndex1 = couple->GetIndex(); 87 fMatIndex1 = couple->GetIndex();
88 break; 88 break;
89 } 89 }
90 } 90 }
91 if(iMat == numOfCouples) 91 if(iMat == numOfCouples)
92 { 92 {
93 G4Exception("G4ForwardXrayTR::G4ForwardXra 93 G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR01",
94 JustWarning, 94 JustWarning,
95 "Invalid first material name i 95 "Invalid first material name in G4ForwardXrayTR constructor!");
96 build = false; 96 build = false;
97 } 97 }
98 98
99 if(build) 99 if(build)
100 { 100 {
101 for(iMat = 0; iMat < numOfCouples; ++iMat) 101 for(iMat = 0; iMat < numOfCouples; ++iMat) // check second material name
102 { 102 {
103 const G4MaterialCutsCouple* couple = 103 const G4MaterialCutsCouple* couple =
104 theCoupleTable->GetMaterialCutsCouple( 104 theCoupleTable->GetMaterialCutsCouple(iMat);
105 if(matName2 == couple->GetMaterial()->Ge 105 if(matName2 == couple->GetMaterial()->GetName())
106 { 106 {
107 fMatIndex2 = couple->GetIndex(); 107 fMatIndex2 = couple->GetIndex();
108 break; 108 break;
109 } 109 }
110 } 110 }
111 if(iMat == numOfCouples) 111 if(iMat == numOfCouples)
112 { 112 {
113 G4Exception( 113 G4Exception(
114 "G4ForwardXrayTR::G4ForwardXrayTR", "F 114 "G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR02", JustWarning,
115 "Invalid second material name in G4For 115 "Invalid second material name in G4ForwardXrayTR constructor!");
116 build = false; 116 build = false;
117 } 117 }
118 } 118 }
119 if(build) 119 if(build)
120 { 120 {
121 BuildXrayTRtables(); 121 BuildXrayTRtables();
122 } 122 }
123 } 123 }
124 124
125 ////////////////////////////////////////////// 125 /////////////////////////////////////////////////////////////////////////
126 // Constructor used by X-ray transition radiat 126 // Constructor used by X-ray transition radiation parametrisation models
127 G4ForwardXrayTR::G4ForwardXrayTR(const G4Strin 127 G4ForwardXrayTR::G4ForwardXrayTR(const G4String& processName)
128 : G4TransitionRadiation(processName) 128 : G4TransitionRadiation(processName)
129 { 129 {
130 fPtrGamma = nullptr; 130 fPtrGamma = nullptr;
131 fGammaCutInKineticEnergy = nullptr; 131 fGammaCutInKineticEnergy = nullptr;
132 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR 132 fGammaTkinCut = fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = 0.0;
133 fGamma = fSigma1 = fSigma2 = 0.0; 133 fGamma = fSigma1 = fSigma2 = 0.0;
134 fAngleDistrTable = nullptr; 134 fAngleDistrTable = nullptr;
135 fEnergyDistrTable = nullptr; 135 fEnergyDistrTable = nullptr;
136 fMatIndex1 = fMatIndex2 = 0; 136 fMatIndex1 = fMatIndex2 = 0;
137 137
138 // Proton energy vector initialization 138 // Proton energy vector initialization
139 fProtonEnergyVector = 139 fProtonEnergyVector =
140 new G4PhysicsLogVector(fMinProtonTkin, fMa 140 new G4PhysicsLogVector(fMinProtonTkin, fMaxProtonTkin, fTotBin);
141 } 141 }
142 142
143 ////////////////////////////////////////////// 143 //////////////////////////////////////////////////////////////////////
144 // Destructor 144 // Destructor
145 G4ForwardXrayTR::~G4ForwardXrayTR() 145 G4ForwardXrayTR::~G4ForwardXrayTR()
146 { 146 {
147 delete fAngleDistrTable; 147 delete fAngleDistrTable;
148 delete fEnergyDistrTable; 148 delete fEnergyDistrTable;
149 delete fProtonEnergyVector; 149 delete fProtonEnergyVector;
150 } 150 }
151 151
152 void G4ForwardXrayTR::ProcessDescription(std:: 152 void G4ForwardXrayTR::ProcessDescription(std::ostream& out) const
153 { 153 {
154 out << "Simulation of forward X-ray transiti 154 out << "Simulation of forward X-ray transition radiation generated by\n"
155 "relativistic charged particles cross 155 "relativistic charged particles crossing the interface between\n"
156 "two materials.\n"; 156 "two materials.\n";
157 } 157 }
158 158
159 G4double G4ForwardXrayTR::GetMeanFreePath(cons 159 G4double G4ForwardXrayTR::GetMeanFreePath(const G4Track&, G4double,
160 G4Fo 160 G4ForceCondition* condition)
161 { 161 {
162 *condition = Forced; 162 *condition = Forced;
163 return DBL_MAX; // so TR doesn't limit mean 163 return DBL_MAX; // so TR doesn't limit mean free path
164 } 164 }
165 165
166 ////////////////////////////////////////////// 166 //////////////////////////////////////////////////////////////////////////////
167 // Build physics tables for energy and angular 167 // Build physics tables for energy and angular distributions of X-ray TR photon
168 void G4ForwardXrayTR::BuildXrayTRtables() 168 void G4ForwardXrayTR::BuildXrayTRtables()
169 { 169 {
170 G4int iMat, jMat, iTkin, iTR, iPlace; 170 G4int iMat, jMat, iTkin, iTR, iPlace;
171 const G4ProductionCutsTable* theCoupleTable 171 const G4ProductionCutsTable* theCoupleTable =
172 G4ProductionCutsTable::GetProductionCutsTa 172 G4ProductionCutsTable::GetProductionCutsTable();
173 G4int numOfCouples = (G4int)theCoupleTable-> 173 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
174 174
175 fGammaCutInKineticEnergy = theCoupleTable->G 175 fGammaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
176 176
177 fAngleDistrTable = new G4PhysicsTable(2 * f 177 fAngleDistrTable = new G4PhysicsTable(2 * fTotBin);
178 fEnergyDistrTable = new G4PhysicsTable(2 * f 178 fEnergyDistrTable = new G4PhysicsTable(2 * fTotBin);
179 179
180 for(iMat = 0; iMat < numOfCouples; 180 for(iMat = 0; iMat < numOfCouples;
181 ++iMat) // loop over pairs of different 181 ++iMat) // loop over pairs of different materials
182 { 182 {
183 if(iMat != fMatIndex1 && iMat != fMatIndex 183 if(iMat != fMatIndex1 && iMat != fMatIndex2)
184 continue; 184 continue;
185 185
186 for(jMat = 0; jMat < numOfCouples; ++jMat) 186 for(jMat = 0; jMat < numOfCouples; ++jMat) // transition iMat -> jMat !!!
187 { 187 {
188 if(iMat == jMat || (jMat != fMatIndex1 & 188 if(iMat == jMat || (jMat != fMatIndex1 && jMat != fMatIndex2))
189 { 189 {
190 continue; 190 continue;
191 } 191 }
192 else 192 else
193 { 193 {
194 const G4MaterialCutsCouple* iCouple = 194 const G4MaterialCutsCouple* iCouple =
195 theCoupleTable->GetMaterialCutsCoupl 195 theCoupleTable->GetMaterialCutsCouple(iMat);
196 const G4MaterialCutsCouple* jCouple = 196 const G4MaterialCutsCouple* jCouple =
197 theCoupleTable->GetMaterialCutsCoupl 197 theCoupleTable->GetMaterialCutsCouple(jMat);
198 const G4Material* mat1 = iCouple->GetM 198 const G4Material* mat1 = iCouple->GetMaterial();
199 const G4Material* mat2 = jCouple->GetM 199 const G4Material* mat2 = jCouple->GetMaterial();
200 200
201 fSigma1 = fPlasmaCof * (mat1->GetElect 201 fSigma1 = fPlasmaCof * (mat1->GetElectronDensity());
202 fSigma2 = fPlasmaCof * (mat2->GetElect 202 fSigma2 = fPlasmaCof * (mat2->GetElectronDensity());
203 203
204 fGammaTkinCut = 0.0; 204 fGammaTkinCut = 0.0;
205 205
206 if(fGammaTkinCut > fTheMinEnergyTR) / 206 if(fGammaTkinCut > fTheMinEnergyTR) // setting of min/max TR energies
207 { 207 {
208 fMinEnergyTR = fGammaTkinCut; 208 fMinEnergyTR = fGammaTkinCut;
209 } 209 }
210 else 210 else
211 { 211 {
212 fMinEnergyTR = fTheMinEnergyTR; 212 fMinEnergyTR = fTheMinEnergyTR;
213 } 213 }
214 if(fGammaTkinCut > fTheMaxEnergyTR) 214 if(fGammaTkinCut > fTheMaxEnergyTR)
215 { 215 {
216 fMaxEnergyTR = 2.0 * fGammaTkinCut; 216 fMaxEnergyTR = 2.0 * fGammaTkinCut; // usually very low TR rate
217 } 217 }
218 else 218 else
219 { 219 {
220 fMaxEnergyTR = fTheMaxEnergyTR; 220 fMaxEnergyTR = fTheMaxEnergyTR;
221 } 221 }
222 for(iTkin = 0; iTkin < fTotBin; ++iTki 222 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // Lorentz factor loop
223 { 223 {
224 auto energyVector = 224 auto energyVector =
225 new G4PhysicsLogVector(fMinEnergyT 225 new G4PhysicsLogVector(fMinEnergyTR, fMaxEnergyTR, fBinTR);
226 226
227 fGamma = 1.0 + (fProtonEnergyVector- 227 fGamma = 1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) /
228 proton_mass_c2); 228 proton_mass_c2);
229 229
230 fMaxThetaTR = 10000.0 / (fGamma * fG 230 fMaxThetaTR = 10000.0 / (fGamma * fGamma);
231 231
232 if(fMaxThetaTR > fTheMaxAngle) 232 if(fMaxThetaTR > fTheMaxAngle)
233 { 233 {
234 fMaxThetaTR = fTheMaxAngle; 234 fMaxThetaTR = fTheMaxAngle;
235 } 235 }
236 else 236 else
237 { 237 {
238 if(fMaxThetaTR < fTheMinAngle) 238 if(fMaxThetaTR < fTheMinAngle)
239 { 239 {
240 fMaxThetaTR = fTheMinAngle; 240 fMaxThetaTR = fTheMinAngle;
241 } 241 }
242 } 242 }
243 auto angleVector = 243 auto angleVector =
244 new G4PhysicsLinearVector(0.0, fMa 244 new G4PhysicsLinearVector(0.0, fMaxThetaTR, fBinTR);
245 G4double energySum = 0.0; 245 G4double energySum = 0.0;
246 G4double angleSum = 0.0; 246 G4double angleSum = 0.0;
247 247
248 energyVector->PutValue(fBinTR - 1, e 248 energyVector->PutValue(fBinTR - 1, energySum);
249 angleVector->PutValue(fBinTR - 1, an 249 angleVector->PutValue(fBinTR - 1, angleSum);
250 250
251 for(iTR = fBinTR - 2; iTR >= 0; --iT 251 for(iTR = fBinTR - 2; iTR >= 0; --iTR)
252 { 252 {
253 energySum += 253 energySum +=
254 fCofTR * EnergySum(energyVector- 254 fCofTR * EnergySum(energyVector->GetLowEdgeEnergy(iTR),
255 energyVector- 255 energyVector->GetLowEdgeEnergy(iTR + 1));
256 256
257 angleSum += 257 angleSum +=
258 fCofTR * AngleSum(angleVector->G 258 fCofTR * AngleSum(angleVector->GetLowEdgeEnergy(iTR),
259 angleVector->G 259 angleVector->GetLowEdgeEnergy(iTR + 1));
260 260
261 energyVector->PutValue(iTR, energy 261 energyVector->PutValue(iTR, energySum);
262 angleVector->PutValue(iTR, angleSu 262 angleVector->PutValue(iTR, angleSum);
263 } 263 }
264 264
265 if(jMat < iMat) 265 if(jMat < iMat)
266 { 266 {
267 iPlace = fTotBin + iTkin; 267 iPlace = fTotBin + iTkin;
268 } 268 }
269 else // jMat > iMat right part of m 269 else // jMat > iMat right part of matrices (jMat-1) !
270 { 270 {
271 iPlace = iTkin; 271 iPlace = iTkin;
272 } 272 }
273 fEnergyDistrTable->insertAt(iPlace, 273 fEnergyDistrTable->insertAt(iPlace, energyVector);
274 fAngleDistrTable->insertAt(iPlace, a 274 fAngleDistrTable->insertAt(iPlace, angleVector);
275 } // iTkin 275 } // iTkin
276 } // jMat != iMat 276 } // jMat != iMat
277 } // jMat 277 } // jMat
278 } // iMat 278 } // iMat
279 } 279 }
280 280
281 ////////////////////////////////////////////// 281 ///////////////////////////////////////////////////////////////////////
282 // 282 //
283 // This function returns the spectral and angl 283 // This function returns the spectral and angle density of TR quanta
284 // in X-ray energy region generated forward wh 284 // in X-ray energy region generated forward when a relativistic
285 // charged particle crosses interface between 285 // charged particle crosses interface between two materials.
286 // The high energy small theta approximation i 286 // The high energy small theta approximation is applied.
287 // (matter1 -> matter2) 287 // (matter1 -> matter2)
288 // varAngle =2* (1 - std::cos(Theta)) or appro 288 // varAngle =2* (1 - std::cos(Theta)) or approximately = Theta*Theta
289 // 289 //
290 G4double G4ForwardXrayTR::SpectralAngleTRdensi 290 G4double G4ForwardXrayTR::SpectralAngleTRdensity(G4double energy,
291 291 G4double varAngle) const
292 { 292 {
293 G4double formationLength1, formationLength2; 293 G4double formationLength1, formationLength2;
294 formationLength1 = 294 formationLength1 =
295 1.0 / (1.0 / (fGamma * fGamma) + fSigma1 / 295 1.0 / (1.0 / (fGamma * fGamma) + fSigma1 / (energy * energy) + varAngle);
296 formationLength2 = 296 formationLength2 =
297 1.0 / (1.0 / (fGamma * fGamma) + fSigma2 / 297 1.0 / (1.0 / (fGamma * fGamma) + fSigma2 / (energy * energy) + varAngle);
298 return (varAngle / energy) * (formationLengt 298 return (varAngle / energy) * (formationLength1 - formationLength2) *
299 (formationLength1 - formationLength2) 299 (formationLength1 - formationLength2);
300 } 300 }
301 301
302 ////////////////////////////////////////////// 302 //////////////////////////////////////////////////////////////////
303 // Analytical formula for angular density of X 303 // Analytical formula for angular density of X-ray TR photons
304 G4double G4ForwardXrayTR::AngleDensity(G4doubl 304 G4double G4ForwardXrayTR::AngleDensity(G4double energy, G4double varAngle) const
305 { 305 {
306 G4double x, x2, c, d, f, a2, b2, a4, b4; 306 G4double x, x2, c, d, f, a2, b2, a4, b4;
307 G4double cof1, cof2, cof3; 307 G4double cof1, cof2, cof3;
308 x = 1.0 / energy; 308 x = 1.0 / energy;
309 x2 = x * x; 309 x2 = x * x;
310 c = 1.0 / fSigma1; 310 c = 1.0 / fSigma1;
311 d = 1.0 / fSigma2; 311 d = 1.0 / fSigma2;
312 f = (varAngle + 1.0 / (fGamma * fGamma)); 312 f = (varAngle + 1.0 / (fGamma * fGamma));
313 a2 = c * f; 313 a2 = c * f;
314 b2 = d * f; 314 b2 = d * f;
315 a4 = a2 * a2; 315 a4 = a2 * a2;
316 b4 = b2 * b2; 316 b4 = b2 * b2;
317 cof1 = c * c * (0.5 / (a2 * (x2 + a2)) + 0.5 317 cof1 = c * c * (0.5 / (a2 * (x2 + a2)) + 0.5 * std::log(x2 / (x2 + a2)) / a4);
318 cof3 = d * d * (0.5 / (b2 * (x2 + b2)) + 0.5 318 cof3 = d * d * (0.5 / (b2 * (x2 + b2)) + 0.5 * std::log(x2 / (x2 + b2)) / b4);
319 cof2 = -c * d * 319 cof2 = -c * d *
320 (std::log(x2 / (x2 + b2)) / b2 - std: 320 (std::log(x2 / (x2 + b2)) / b2 - std::log(x2 / (x2 + a2)) / a2) /
321 (a2 - b2); 321 (a2 - b2);
322 return -varAngle * (cof1 + cof2 + cof3); 322 return -varAngle * (cof1 + cof2 + cof3);
323 } 323 }
324 324
325 ////////////////////////////////////////////// 325 /////////////////////////////////////////////////////////////////////
326 // Definite integral of X-ray TR spectral-angl 326 // Definite integral of X-ray TR spectral-angle density from energy1
327 // to energy2 327 // to energy2
328 G4double G4ForwardXrayTR::EnergyInterval(G4dou 328 G4double G4ForwardXrayTR::EnergyInterval(G4double energy1, G4double energy2,
329 G4dou 329 G4double varAngle) const
330 { 330 {
331 return AngleDensity(energy2, varAngle) - Ang 331 return AngleDensity(energy2, varAngle) - AngleDensity(energy1, varAngle);
332 } 332 }
333 333
334 ////////////////////////////////////////////// 334 //////////////////////////////////////////////////////////////////////
335 // Integral angle distribution of X-ray TR pho 335 // Integral angle distribution of X-ray TR photons based on analytical
336 // formula for angle density 336 // formula for angle density
337 G4double G4ForwardXrayTR::AngleSum(G4double va 337 G4double G4ForwardXrayTR::AngleSum(G4double varAngle1, G4double varAngle2) const
338 { 338 {
339 G4int i; 339 G4int i;
340 G4double h, sumEven = 0.0, sumOdd = 0.0; 340 G4double h, sumEven = 0.0, sumOdd = 0.0;
341 h = 0.5 * (varAngle2 - varAngle1) / fSympson 341 h = 0.5 * (varAngle2 - varAngle1) / fSympsonNumber;
342 for(i = 1; i < fSympsonNumber; ++i) 342 for(i = 1; i < fSympsonNumber; ++i)
343 { 343 {
344 sumEven += 344 sumEven +=
345 EnergyInterval(fMinEnergyTR, fMaxEnergyT 345 EnergyInterval(fMinEnergyTR, fMaxEnergyTR, varAngle1 + 2 * i * h);
346 sumOdd += 346 sumOdd +=
347 EnergyInterval(fMinEnergyTR, fMaxEnergyT 347 EnergyInterval(fMinEnergyTR, fMaxEnergyTR, varAngle1 + (2 * i - 1) * h);
348 } 348 }
349 sumOdd += EnergyInterval(fMinEnergyTR, fMaxE 349 sumOdd += EnergyInterval(fMinEnergyTR, fMaxEnergyTR,
350 varAngle1 + (2 * fS 350 varAngle1 + (2 * fSympsonNumber - 1) * h);
351 351
352 return h * 352 return h *
353 (EnergyInterval(fMinEnergyTR, fMaxEne 353 (EnergyInterval(fMinEnergyTR, fMaxEnergyTR, varAngle1) +
354 EnergyInterval(fMinEnergyTR, fMaxEne 354 EnergyInterval(fMinEnergyTR, fMaxEnergyTR, varAngle2) + 4.0 * sumOdd +
355 2.0 * sumEven) / 355 2.0 * sumEven) /
356 3.0; 356 3.0;
357 } 357 }
358 358
359 ////////////////////////////////////////////// 359 /////////////////////////////////////////////////////////////////////
360 // Analytical Expression for spectral densit 360 // Analytical Expression for spectral density of Xray TR photons
361 // x = 2*(1 - std::cos(Theta)) ~ Theta^2 361 // x = 2*(1 - std::cos(Theta)) ~ Theta^2
362 G4double G4ForwardXrayTR::SpectralDensity(G4do 362 G4double G4ForwardXrayTR::SpectralDensity(G4double energy, G4double x) const
363 { 363 {
364 G4double a, b; 364 G4double a, b;
365 a = 1.0 / (fGamma * fGamma) + fSigma1 / (ene 365 a = 1.0 / (fGamma * fGamma) + fSigma1 / (energy * energy);
366 b = 1.0 / (fGamma * fGamma) + fSigma2 / (ene 366 b = 1.0 / (fGamma * fGamma) + fSigma2 / (energy * energy);
367 return ((a + b) * std::log((x + b) / (x + a) 367 return ((a + b) * std::log((x + b) / (x + a)) / (a - b) + a / (x + a) +
368 b / (x + b)) / 368 b / (x + b)) /
369 energy; 369 energy;
370 } 370 }
371 371
372 ////////////////////////////////////////////// 372 ////////////////////////////////////////////////////////////////////
373 // The spectral density in some angle interva 373 // The spectral density in some angle interval from varAngle1 to
374 // varAngle2 374 // varAngle2
375 G4double G4ForwardXrayTR::AngleInterval(G4doub 375 G4double G4ForwardXrayTR::AngleInterval(G4double energy, G4double varAngle1,
376 G4doub 376 G4double varAngle2) const
377 { 377 {
378 return SpectralDensity(energy, varAngle2) - 378 return SpectralDensity(energy, varAngle2) -
379 SpectralDensity(energy, varAngle1); 379 SpectralDensity(energy, varAngle1);
380 } 380 }
381 381
382 ////////////////////////////////////////////// 382 ////////////////////////////////////////////////////////////////////
383 // Integral spectral distribution of X-ray TR 383 // Integral spectral distribution of X-ray TR photons based on
384 // analytical formula for spectral density 384 // analytical formula for spectral density
385 G4double G4ForwardXrayTR::EnergySum(G4double e 385 G4double G4ForwardXrayTR::EnergySum(G4double energy1, G4double energy2) const
386 { 386 {
387 G4int i; 387 G4int i;
388 G4double h, sumEven = 0.0, sumOdd = 0.0; 388 G4double h, sumEven = 0.0, sumOdd = 0.0;
389 h = 0.5 * (energy2 - energy1) / fSympsonNumb 389 h = 0.5 * (energy2 - energy1) / fSympsonNumber;
390 for(i = 1; i < fSympsonNumber; ++i) 390 for(i = 1; i < fSympsonNumber; ++i)
391 { 391 {
392 sumEven += AngleInterval(energy1 + 2 * i * 392 sumEven += AngleInterval(energy1 + 2 * i * h, 0.0, fMaxThetaTR);
393 sumOdd += AngleInterval(energy1 + (2 * i - 393 sumOdd += AngleInterval(energy1 + (2 * i - 1) * h, 0.0, fMaxThetaTR);
394 } 394 }
395 sumOdd += 395 sumOdd +=
396 AngleInterval(energy1 + (2 * fSympsonNumbe 396 AngleInterval(energy1 + (2 * fSympsonNumber - 1) * h, 0.0, fMaxThetaTR);
397 397
398 return h * 398 return h *
399 (AngleInterval(energy1, 0.0, fMaxThet 399 (AngleInterval(energy1, 0.0, fMaxThetaTR) +
400 AngleInterval(energy2, 0.0, fMaxThet 400 AngleInterval(energy2, 0.0, fMaxThetaTR) + 4.0 * sumOdd +
401 2.0 * sumEven) / 401 2.0 * sumEven) /
402 3.0; 402 3.0;
403 } 403 }
404 404
405 ////////////////////////////////////////////// 405 /////////////////////////////////////////////////////////////////////////
406 // PostStepDoIt function for creation of forwa 406 // PostStepDoIt function for creation of forward X-ray photons in TR process
407 // on boundary between two materials with real 407 // on boundary between two materials with really different plasma energies
408 G4VParticleChange* G4ForwardXrayTR::PostStepDo 408 G4VParticleChange* G4ForwardXrayTR::PostStepDoIt(const G4Track& aTrack,
409 409 const G4Step& aStep)
410 { 410 {
411 aParticleChange.Initialize(aTrack); 411 aParticleChange.Initialize(aTrack);
412 G4int iMat, jMat, iTkin, iPlace, numOfTR, iT 412 G4int iMat, jMat, iTkin, iPlace, numOfTR, iTR, iTransfer;
413 413
414 G4double energyPos, anglePos, energyTR, thet 414 G4double energyPos, anglePos, energyTR, theta, phi, dirX, dirY, dirZ;
415 G4double W, W1, W2, E1, E2; 415 G4double W, W1, W2, E1, E2;
416 416
417 G4StepPoint* pPreStepPoint = aStep.GetPreSt 417 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
418 G4StepPoint* pPostStepPoint = aStep.GetPostS 418 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
419 G4double tol = 419 G4double tol =
420 0.5 * G4GeometryTolerance::GetInstance()-> 420 0.5 * G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
421 421
422 if(pPostStepPoint->GetStepStatus() != fGeomB 422 if(pPostStepPoint->GetStepStatus() != fGeomBoundary)
423 { 423 {
424 return G4VDiscreteProcess::PostStepDoIt(aT 424 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
425 } 425 }
426 if(aTrack.GetStepLength() <= tol) 426 if(aTrack.GetStepLength() <= tol)
427 { 427 {
428 return G4VDiscreteProcess::PostStepDoIt(aT 428 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
429 } 429 }
430 // Arrived at boundary, so begin to try TR 430 // Arrived at boundary, so begin to try TR
431 431
432 const G4MaterialCutsCouple* iCouple = pPreSt 432 const G4MaterialCutsCouple* iCouple = pPreStepPoint->GetPhysicalVolume()
433 ->Ge 433 ->GetLogicalVolume()
434 ->Ge 434 ->GetMaterialCutsCouple();
435 const G4MaterialCutsCouple* jCouple = pPostS 435 const G4MaterialCutsCouple* jCouple = pPostStepPoint->GetPhysicalVolume()
436 ->Ge 436 ->GetLogicalVolume()
437 ->Ge 437 ->GetMaterialCutsCouple();
438 const G4Material* iMaterial = iCouple->GetMa 438 const G4Material* iMaterial = iCouple->GetMaterial();
439 const G4Material* jMaterial = jCouple->GetMa 439 const G4Material* jMaterial = jCouple->GetMaterial();
440 iMat = iCouple->GetIn 440 iMat = iCouple->GetIndex();
441 jMat = jCouple->GetIn 441 jMat = jCouple->GetIndex();
442 442
443 // The case of equal or approximate (in term 443 // The case of equal or approximate (in terms of plasma energy) materials
444 // No TR photons ?! 444 // No TR photons ?!
445 445
446 if(iMat == jMat || 446 if(iMat == jMat ||
447 ((fMatIndex1 >= 0 && fMatIndex2 >= 0) && 447 ((fMatIndex1 >= 0 && fMatIndex2 >= 0) &&
448 (iMat != fMatIndex1 && iMat != fMatIndex 448 (iMat != fMatIndex1 && iMat != fMatIndex2) &&
449 (jMat != fMatIndex1 && jMat != fMatIndex 449 (jMat != fMatIndex1 && jMat != fMatIndex2))
450 450
451 || iMaterial->GetState() == jMaterial->Ge 451 || iMaterial->GetState() == jMaterial->GetState()
452 452
453 || (iMaterial->GetState() == kStateSolid 453 || (iMaterial->GetState() == kStateSolid &&
454 jMaterial->GetState() == kStateLiquid 454 jMaterial->GetState() == kStateLiquid)
455 455
456 || (iMaterial->GetState() == kStateLiquid 456 || (iMaterial->GetState() == kStateLiquid &&
457 jMaterial->GetState() == kStateSolid) 457 jMaterial->GetState() == kStateSolid))
458 { 458 {
459 return G4VDiscreteProcess::PostStepDoIt(aT 459 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
460 } 460 }
461 461
462 const G4DynamicParticle* aParticle = aTrack. 462 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
463 G4double charge = aParticle->GetDefinition() 463 G4double charge = aParticle->GetDefinition()->GetPDGCharge();
464 464
465 if(charge == 0.0) // Uncharged particle doe 465 if(charge == 0.0) // Uncharged particle doesn't Generate TR photons
466 { 466 {
467 return G4VDiscreteProcess::PostStepDoIt(aT 467 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
468 } 468 }
469 // Now we are ready to Generate TR photons 469 // Now we are ready to Generate TR photons
470 470
471 G4double chargeSq = charge * charge; 471 G4double chargeSq = charge * charge;
472 G4double kinEnergy = aParticle->GetKineticEn 472 G4double kinEnergy = aParticle->GetKineticEnergy();
473 G4double massRatio = 473 G4double massRatio =
474 proton_mass_c2 / aParticle->GetDefinition( 474 proton_mass_c2 / aParticle->GetDefinition()->GetPDGMass();
475 G4double TkinScaled = kinEnergy * massRatio; 475 G4double TkinScaled = kinEnergy * massRatio;
476 for(iTkin = 0; iTkin < fTotBin; ++iTkin) 476 for(iTkin = 0; iTkin < fTotBin; ++iTkin)
477 { 477 {
478 if(TkinScaled < fProtonEnergyVector->GetLo 478 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
479 { 479 {
480 break; 480 break;
481 } 481 }
482 } 482 }
483 if(jMat < iMat) 483 if(jMat < iMat)
484 { 484 {
485 iPlace = fTotBin + iTkin - 1; 485 iPlace = fTotBin + iTkin - 1;
486 } 486 }
487 else 487 else
488 { 488 {
489 iPlace = iTkin - 1; 489 iPlace = iTkin - 1;
490 } 490 }
491 491
492 G4ParticleMomentum particleDir = aParticle-> 492 G4ParticleMomentum particleDir = aParticle->GetMomentumDirection();
493 493
494 if(iTkin == fTotBin) // TR plato, try from 494 if(iTkin == fTotBin) // TR plato, try from left
495 { 495 {
496 numOfTR = (G4int)G4Poisson( 496 numOfTR = (G4int)G4Poisson(
497 ((*(*fEnergyDistrTable)(iPlace))(0) + (* 497 ((*(*fEnergyDistrTable)(iPlace))(0) + (*(*fAngleDistrTable)(iPlace))(0)) *
498 chargeSq * 0.5); 498 chargeSq * 0.5);
499 if(numOfTR == 0) 499 if(numOfTR == 0)
500 { 500 {
501 return G4VDiscreteProcess::PostStepDoIt( 501 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
502 } 502 }
503 else 503 else
504 { 504 {
505 aParticleChange.SetNumberOfSecondaries(n 505 aParticleChange.SetNumberOfSecondaries(numOfTR);
506 506
507 for(iTR = 0; iTR < numOfTR; ++iTR) 507 for(iTR = 0; iTR < numOfTR; ++iTR)
508 { 508 {
509 energyPos = (*(*fEnergyDistrTable)(iPl 509 energyPos = (*(*fEnergyDistrTable)(iPlace))(0) * G4UniformRand();
510 for(iTransfer = 0; iTransfer < fBinTR 510 for(iTransfer = 0; iTransfer < fBinTR - 1; ++iTransfer)
511 { 511 {
512 if(energyPos >= (*(*fEnergyDistrTabl 512 if(energyPos >= (*(*fEnergyDistrTable)(iPlace))(iTransfer))
513 break; 513 break;
514 } 514 }
515 energyTR = (*fEnergyDistrTable)(iPlace 515 energyTR = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
516 516
517 kinEnergy -= energyTR; 517 kinEnergy -= energyTR;
518 aParticleChange.ProposeEnergy(kinEnerg 518 aParticleChange.ProposeEnergy(kinEnergy);
519 519
520 anglePos = (*(*fAngleDistrTable)(iPlac 520 anglePos = (*(*fAngleDistrTable)(iPlace))(0) * G4UniformRand();
521 for(iTransfer = 0; iTransfer < fBinTR 521 for(iTransfer = 0; iTransfer < fBinTR - 1; ++iTransfer)
522 { 522 {
523 if(anglePos > (*(*fAngleDistrTable)( 523 if(anglePos > (*(*fAngleDistrTable)(iPlace))(iTransfer))
524 break; 524 break;
525 } 525 }
526 theta = std::sqrt( 526 theta = std::sqrt(
527 (*fAngleDistrTable)(iPlace)->GetLowE 527 (*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1));
528 528
529 phi = twopi * G4UniformRand(); 529 phi = twopi * G4UniformRand();
530 dirX = std::sin(theta) * std::cos(phi) 530 dirX = std::sin(theta) * std::cos(phi);
531 dirY = std::sin(theta) * std::sin(phi) 531 dirY = std::sin(theta) * std::sin(phi);
532 dirZ = std::cos(theta); 532 dirZ = std::cos(theta);
533 G4ThreeVector directionTR(dirX, dirY, 533 G4ThreeVector directionTR(dirX, dirY, dirZ);
534 directionTR.rotateUz(particleDir); 534 directionTR.rotateUz(particleDir);
535 auto aPhotonTR = new G4DynamicParticle 535 auto aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(), directionTR, energyTR);
536 536
537 // Create the G4Track 537 // Create the G4Track
538 auto aSecondaryTrack = new G4Track(aPhotonTR 538 auto aSecondaryTrack = new G4Track(aPhotonTR, aTrack.GetGlobalTime(), aTrack.GetPosition());
539 aSecondaryTrack->SetTouchableHandle(aStep.Ge 539 aSecondaryTrack->SetTouchableHandle(aStep.GetPostStepPoint()->GetTouchableHandle());
540 aSecondaryTrack->SetParentID(aTrack.GetTrack 540 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
541 aSecondaryTrack->SetCreatorModelID(secID); 541 aSecondaryTrack->SetCreatorModelID(secID);
542 aParticleChange.AddSecondary(aSecondaryTrack 542 aParticleChange.AddSecondary(aSecondaryTrack);
543 } 543 }
544 } 544 }
545 } 545 }
546 else 546 else
547 { 547 {
548 if(iTkin == 0) // Tkin is too small, negl 548 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
549 { 549 {
550 return G4VDiscreteProcess::PostStepDoIt( 550 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
551 } 551 }
552 else // general case: Tkin between two ve 552 else // general case: Tkin between two vectors of the material
553 { 553 {
554 E1 = fProtonEnergyVector->GetLowEdgeEner 554 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
555 E2 = fProtonEnergyVector->GetLowEdgeEner 555 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin);
556 W = 1.0 / (E2 - E1); 556 W = 1.0 / (E2 - E1);
557 W1 = (E2 - TkinScaled) * W; 557 W1 = (E2 - TkinScaled) * W;
558 W2 = (TkinScaled - E1) * W; 558 W2 = (TkinScaled - E1) * W;
559 559
560 numOfTR = (G4int)G4Poisson((((*(*fEnergy 560 numOfTR = (G4int)G4Poisson((((*(*fEnergyDistrTable)(iPlace))(0) +
561 (*(*fAngleD 561 (*(*fAngleDistrTable)(iPlace))(0)) *
562 W1 + 562 W1 +
563 ((*(*fEnergy 563 ((*(*fEnergyDistrTable)(iPlace + 1))(0) +
564 (*(*fAngleD 564 (*(*fAngleDistrTable)(iPlace + 1))(0)) *
565 W2) * 565 W2) *
566 chargeSq * 0. 566 chargeSq * 0.5);
567 if(numOfTR == 0) 567 if(numOfTR == 0)
568 { 568 {
569 return G4VDiscreteProcess::PostStepDoI 569 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
570 } 570 }
571 else 571 else
572 { 572 {
573 aParticleChange.SetNumberOfSecondaries 573 aParticleChange.SetNumberOfSecondaries(numOfTR);
574 for(iTR = 0; iTR < numOfTR; ++iTR) 574 for(iTR = 0; iTR < numOfTR; ++iTR)
575 { 575 {
576 energyPos = ((*(*fEnergyDistrTable)( 576 energyPos = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 +
577 (*(*fEnergyDistrTable)( 577 (*(*fEnergyDistrTable)(iPlace + 1))(0) * W2) *
578 G4UniformRand(); 578 G4UniformRand();
579 for(iTransfer = 0; iTransfer < fBinT 579 for(iTransfer = 0; iTransfer < fBinTR - 1; ++iTransfer)
580 { 580 {
581 if(energyPos >= 581 if(energyPos >=
582 ((*(*fEnergyDistrTable)(iPlace) 582 ((*(*fEnergyDistrTable)(iPlace))(iTransfer) *W1 +
583 (*(*fEnergyDistrTable)(iPlace 583 (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer) *W2))
584 break; 584 break;
585 } 585 }
586 energyTR = 586 energyTR =
587 ((*fEnergyDistrTable)(iPlace)->Get 587 ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) * W1 +
588 ((*fEnergyDistrTable)(iPlace + 1)- 588 ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer)) *
589 W2; 589 W2;
590 590
591 kinEnergy -= energyTR; 591 kinEnergy -= energyTR;
592 aParticleChange.ProposeEnergy(kinEne 592 aParticleChange.ProposeEnergy(kinEnergy);
593 593
594 anglePos = ((*(*fAngleDistrTable)(iP 594 anglePos = ((*(*fAngleDistrTable)(iPlace))(0) * W1 +
595 (*(*fAngleDistrTable)(iP 595 (*(*fAngleDistrTable)(iPlace + 1))(0) * W2) *
596 G4UniformRand(); 596 G4UniformRand();
597 for(iTransfer = 0; iTransfer < fBinT 597 for(iTransfer = 0; iTransfer < fBinTR - 1; ++iTransfer)
598 { 598 {
599 if(anglePos > ((*(*fAngleDistrTabl 599 if(anglePos > ((*(*fAngleDistrTable)(iPlace))(iTransfer) *W1 +
600 (*(*fAngleDistrTabl 600 (*(*fAngleDistrTable)(iPlace + 1))(iTransfer) *W2))
601 break; 601 break;
602 } 602 }
603 theta = std::sqrt( 603 theta = std::sqrt(
604 ((*fAngleDistrTable)(iPlace)->GetL 604 ((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1)) *
605 W1 + 605 W1 +
606 ((*fAngleDistrTable)(iPlace + 1)-> 606 ((*fAngleDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer - 1)) *
607 W2); 607 W2);
608 608
609 phi = twopi * G4UniformRand(); 609 phi = twopi * G4UniformRand();
610 dirX = std::sin(theta) * std::cos(ph 610 dirX = std::sin(theta) * std::cos(phi);
611 dirY = std::sin(theta) * std::sin(ph 611 dirY = std::sin(theta) * std::sin(phi);
612 dirZ = std::cos(theta); 612 dirZ = std::cos(theta);
613 G4ThreeVector directionTR(dirX, dirY 613 G4ThreeVector directionTR(dirX, dirY, dirZ);
614 directionTR.rotateUz(particleDir); 614 directionTR.rotateUz(particleDir);
615 auto aPhotonTR = 615 auto aPhotonTR =
616 new G4DynamicParticle(G4Gamma::Gam 616 new G4DynamicParticle(G4Gamma::Gamma(), directionTR, energyTR);
617 617
618 // Create the G4Track 618 // Create the G4Track
619 G4Track* aSecondaryTrack = new G4Track(aPh 619 G4Track* aSecondaryTrack = new G4Track(aPhotonTR, aTrack.GetGlobalTime(), aTrack.GetPosition());
620 aSecondaryTrack->SetTouchableHandle(aStep. 620 aSecondaryTrack->SetTouchableHandle(aStep.GetPostStepPoint()->GetTouchableHandle());
621 aSecondaryTrack->SetParentID(aTrack.GetTra 621 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
622 aSecondaryTrack->SetCreatorModelID(secID); 622 aSecondaryTrack->SetCreatorModelID(secID);
623 aParticleChange.AddSecondary(aSecondaryTra 623 aParticleChange.AddSecondary(aSecondaryTrack);
624 } 624 }
625 } 625 }
626 } 626 }
627 } 627 }
628 return &aParticleChange; 628 return &aParticleChange;
629 } 629 }
630 630
631 ////////////////////////////////////////////// 631 ////////////////////////////////////////////////////////////////////////////
632 // Test function for checking of PostStepDoIt 632 // Test function for checking of PostStepDoIt random preparation of TR photon
633 // energy 633 // energy
634 G4double G4ForwardXrayTR::GetEnergyTR(G4int iM 634 G4double G4ForwardXrayTR::GetEnergyTR(G4int iMat, G4int jMat, G4int iTkin) const
635 { 635 {
636 G4int iPlace, numOfTR, iTR, iTransfer; 636 G4int iPlace, numOfTR, iTR, iTransfer;
637 G4double energyTR = 0.0; // return this val 637 G4double energyTR = 0.0; // return this value for no TR photons
638 G4double energyPos; 638 G4double energyPos;
639 G4double W1, W2; 639 G4double W1, W2;
640 640
641 const G4ProductionCutsTable* theCoupleTable 641 const G4ProductionCutsTable* theCoupleTable =
642 G4ProductionCutsTable::GetProductionCutsTa 642 G4ProductionCutsTable::GetProductionCutsTable();
643 G4int numOfCouples = (G4int)theCoupleTable-> 643 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
644 644
645 // The case of equal or approximate (in term 645 // The case of equal or approximate (in terms of plasma energy) materials
646 // No TR photons ?! 646 // No TR photons ?!
647 647
648 const G4MaterialCutsCouple* iCouple = 648 const G4MaterialCutsCouple* iCouple =
649 theCoupleTable->GetMaterialCutsCouple(iMat 649 theCoupleTable->GetMaterialCutsCouple(iMat);
650 const G4MaterialCutsCouple* jCouple = 650 const G4MaterialCutsCouple* jCouple =
651 theCoupleTable->GetMaterialCutsCouple(jMat 651 theCoupleTable->GetMaterialCutsCouple(jMat);
652 const G4Material* iMaterial = iCouple->GetMa 652 const G4Material* iMaterial = iCouple->GetMaterial();
653 const G4Material* jMaterial = jCouple->GetMa 653 const G4Material* jMaterial = jCouple->GetMaterial();
654 654
655 if(iMat == jMat 655 if(iMat == jMat
656 656
657 || iMaterial->GetState() == jMaterial->Ge 657 || iMaterial->GetState() == jMaterial->GetState()
658 658
659 || (iMaterial->GetState() == kStateSolid 659 || (iMaterial->GetState() == kStateSolid &&
660 jMaterial->GetState() == kStateLiquid 660 jMaterial->GetState() == kStateLiquid)
661 661
662 || (iMaterial->GetState() == kStateLiquid 662 || (iMaterial->GetState() == kStateLiquid &&
663 jMaterial->GetState() == kStateSolid) 663 jMaterial->GetState() == kStateSolid))
664 664
665 { 665 {
666 return energyTR; 666 return energyTR;
667 } 667 }
668 668
669 if(jMat < iMat) 669 if(jMat < iMat)
670 { 670 {
671 iPlace = (iMat * (numOfCouples - 1) + jMat 671 iPlace = (iMat * (numOfCouples - 1) + jMat) * fTotBin + iTkin - 1;
672 } 672 }
673 else 673 else
674 { 674 {
675 iPlace = (iMat * (numOfCouples - 1) + jMat 675 iPlace = (iMat * (numOfCouples - 1) + jMat - 1) * fTotBin + iTkin - 1;
676 } 676 }
677 G4PhysicsVector* energyVector1 = (*fEnergyDi 677 G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace);
678 G4PhysicsVector* energyVector2 = (*fEnergyDi 678 G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1);
679 679
680 if(iTkin == fTotBin) // TR plato, try from 680 if(iTkin == fTotBin) // TR plato, try from left
681 { 681 {
682 numOfTR = (G4int)G4Poisson((*energyVector1 682 numOfTR = (G4int)G4Poisson((*energyVector1)(0));
683 if(numOfTR == 0) 683 if(numOfTR == 0)
684 { 684 {
685 return energyTR; 685 return energyTR;
686 } 686 }
687 else 687 else
688 { 688 {
689 for(iTR = 0; iTR < numOfTR; ++iTR) 689 for(iTR = 0; iTR < numOfTR; ++iTR)
690 { 690 {
691 energyPos = (*energyVector1)(0) * G4Un 691 energyPos = (*energyVector1)(0) * G4UniformRand();
692 for(iTransfer = 0; iTransfer < fBinTR 692 for(iTransfer = 0; iTransfer < fBinTR - 1; ++iTransfer)
693 { 693 {
694 if(energyPos >= (*energyVector1)(iTr 694 if(energyPos >= (*energyVector1)(iTransfer))
695 break; 695 break;
696 } 696 }
697 energyTR += energyVector1->GetLowEdgeE 697 energyTR += energyVector1->GetLowEdgeEnergy(iTransfer);
698 } 698 }
699 } 699 }
700 } 700 }
701 else 701 else
702 { 702 {
703 if(iTkin == 0) // Tkin is too small, negl 703 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
704 { 704 {
705 return energyTR; 705 return energyTR;
706 } 706 }
707 else // general case: Tkin between two ve 707 else // general case: Tkin between two vectors of the material
708 { // use trivial mean half/half 708 { // use trivial mean half/half
709 W1 = 0.5; 709 W1 = 0.5;
710 W2 = 0.5; 710 W2 = 0.5;
711 numOfTR = (G4int)G4Poisson((*energyVecto 711 numOfTR = (G4int)G4Poisson((*energyVector1)(0) * W1 + (*energyVector2)(0) * W2);
712 if(numOfTR == 0) 712 if(numOfTR == 0)
713 { 713 {
714 return energyTR; 714 return energyTR;
715 } 715 }
716 else 716 else
717 { 717 {
718 G4cout << "It is still OK in GetEnergy 718 G4cout << "It is still OK in GetEnergyTR(int,int,int)" << G4endl;
719 for(iTR = 0; iTR < numOfTR; ++iTR) 719 for(iTR = 0; iTR < numOfTR; ++iTR)
720 { 720 {
721 energyPos = ((*energyVector1)(0) * W 721 energyPos = ((*energyVector1)(0) * W1 + (*energyVector2)(0) * W2) *
722 G4UniformRand(); 722 G4UniformRand();
723 for(iTransfer = 0; iTransfer < fBinT 723 for(iTransfer = 0; iTransfer < fBinTR - 1; ++iTransfer)
724 { 724 {
725 if(energyPos >= ((*energyVector1)( 725 if(energyPos >= ((*energyVector1)(iTransfer) *W1 +
726 (*energyVector2)( 726 (*energyVector2)(iTransfer) *W2))
727 break; 727 break;
728 } 728 }
729 energyTR += (energyVector1->GetLowEd 729 energyTR += (energyVector1->GetLowEdgeEnergy(iTransfer)) * W1 +
730 (energyVector2->GetLowEd 730 (energyVector2->GetLowEdgeEnergy(iTransfer)) * W2;
731 } 731 }
732 } 732 }
733 } 733 }
734 } 734 }
735 735
736 return energyTR; 736 return energyTR;
737 } 737 }
738 738
739 ////////////////////////////////////////////// 739 ////////////////////////////////////////////////////////////////////////////
740 // Test function for checking of PostStepDoIt 740 // Test function for checking of PostStepDoIt random preparation of TR photon
741 // theta angle relative to particle direction 741 // theta angle relative to particle direction
742 G4double G4ForwardXrayTR::GetThetaTR(G4int, G4 742 G4double G4ForwardXrayTR::GetThetaTR(G4int, G4int, G4int) const { return 0.0; }
743 743
744 G4int G4ForwardXrayTR::GetSympsonNumber() { re 744 G4int G4ForwardXrayTR::GetSympsonNumber() { return fSympsonNumber; }
745 745
746 G4int G4ForwardXrayTR::GetBinTR() { return fBi 746 G4int G4ForwardXrayTR::GetBinTR() { return fBinTR; }
747 747
748 G4double G4ForwardXrayTR::GetMinProtonTkin() { 748 G4double G4ForwardXrayTR::GetMinProtonTkin() { return fMinProtonTkin; }
749 749
750 G4double G4ForwardXrayTR::GetMaxProtonTkin() { 750 G4double G4ForwardXrayTR::GetMaxProtonTkin() { return fMaxProtonTkin; }
751 751
752 G4int G4ForwardXrayTR::GetTotBin() { return fT 752 G4int G4ForwardXrayTR::GetTotBin() { return fTotBin; }
753 753