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