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