| 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: G4VXTRenergyLoss.cc,v 1.21 2005/10/11 14:24:34 grichine Exp $
>> 25 // GEANT4 tag $Name: geant4-08-00 $
>> 26 //
26 // History: 27 // History:
27 // 2001-2002 R&D by V.Grichine 28 // 2001-2002 R&D by V.Grichine
28 // 19.06.03 V. Grichine, modifications in Buil << 29 // 19.06.03 V. Grichine, modifications in BuildTable for the integration
29 // in respect of angle: 30 // in respect of angle: range is increased, accuracy is
30 // improved 31 // improved
31 // 28.07.05, P.Gumplinger add G4ProcessType to 32 // 28.07.05, P.Gumplinger add G4ProcessType to constructor
32 // 28.09.07, V.Ivanchenko general cleanup with <<
33 // 33 //
34 34
35 #include "G4VXTRenergyLoss.hh" << 35 #include "G4Timer.hh"
36 36
37 #include "G4AffineTransform.hh" << 37 #include "G4VXTRenergyLoss.hh"
38 #include "G4DynamicParticle.hh" << 38 #include "G4Poisson.hh"
39 #include "G4EmProcessSubType.hh" <<
40 #include "G4Integrator.hh" <<
41 #include "G4MaterialTable.hh" 39 #include "G4MaterialTable.hh"
42 #include "G4ParticleMomentum.hh" <<
43 #include "G4PhysicalConstants.hh" <<
44 #include "G4PhysicsFreeVector.hh" <<
45 #include "G4PhysicsLinearVector.hh" <<
46 #include "G4PhysicsLogVector.hh" <<
47 #include "G4RotationMatrix.hh" <<
48 #include "G4SandiaTable.hh" <<
49 #include "G4SystemOfUnits.hh" <<
50 #include "G4ThreeVector.hh" <<
51 #include "G4Timer.hh" <<
52 #include "G4VDiscreteProcess.hh" 40 #include "G4VDiscreteProcess.hh"
53 #include "G4VParticleChange.hh" 41 #include "G4VParticleChange.hh"
54 #include "G4VSolid.hh" 42 #include "G4VSolid.hh"
55 #include "G4PhysicsModelCatalog.hh" << 43
>> 44 #include "G4RotationMatrix.hh"
>> 45 #include "G4ThreeVector.hh"
>> 46 #include "G4AffineTransform.hh"
>> 47
>> 48 #include "G4PhysicsVector.hh"
>> 49 #include "G4PhysicsLinearVector.hh"
>> 50
>> 51 using namespace std;
>> 52
>> 53 // Initialization of local constants
>> 54
>> 55 G4double G4XTRenergyLoss::fTheMinEnergyTR = 1.0*keV ;
>> 56 G4double G4XTRenergyLoss::fTheMaxEnergyTR = 100.0*keV ;
>> 57 G4double G4XTRenergyLoss::fTheMaxAngle = 1.0e-3 ;
>> 58 G4double G4XTRenergyLoss::fTheMinAngle = 5.0e-6 ;
>> 59 G4int G4XTRenergyLoss::fBinTR = 50 ;
>> 60
>> 61 G4double G4XTRenergyLoss::fMinProtonTkin = 100.0*GeV ;
>> 62 G4double G4XTRenergyLoss::fMaxProtonTkin = 100.0*TeV ;
>> 63 G4int G4XTRenergyLoss::fTotBin = 50 ;
>> 64 // Proton energy vector initialization
>> 65
>> 66 G4PhysicsLogVector* G4XTRenergyLoss::
>> 67 fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin,
>> 68 fMaxProtonTkin,
>> 69 fTotBin ) ;
>> 70
>> 71 G4double G4XTRenergyLoss::fPlasmaCof = 4.0*pi*fine_structure_const*
>> 72 hbarc*hbarc*hbarc/electron_mass_c2 ;
>> 73
>> 74 G4double G4XTRenergyLoss::fCofTR = fine_structure_const/pi ;
>> 75
>> 76
>> 77
>> 78
56 79
57 ////////////////////////////////////////////// 80 ////////////////////////////////////////////////////////////////////////////
>> 81 //
58 // Constructor, destructor 82 // Constructor, destructor
59 G4VXTRenergyLoss::G4VXTRenergyLoss(G4LogicalVo <<
60 G4Material* <<
61 G4double a, <<
62 const G4Str <<
63 G4ProcessTy <<
64 : G4VDiscreteProcess(processName, type) <<
65 , fGammaCutInKineticEnergy(nullptr) <<
66 , fAngleDistrTable(nullptr) <<
67 , fEnergyDistrTable(nullptr) <<
68 , fAngleForEnergyTable(nullptr) <<
69 , fPlatePhotoAbsCof(nullptr) <<
70 , fGasPhotoAbsCof(nullptr) <<
71 , fGammaTkinCut(0.0) <<
72 { <<
73 verboseLevel = 1; <<
74 secID = G4PhysicsModelCatalog::GetModelID("m <<
75 SetProcessSubType(fTransitionRadiation); <<
76 <<
77 fPtrGamma = nullptr; <<
78 fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = <<
79 fVarAngle = fLambda = fTotalDist = fPlateThi <<
80 fAlphaPlate = 100.; <<
81 fAlphaGas = 40.; <<
82 <<
83 fTheMinEnergyTR = CLHEP::keV * 1.; // 1.; / <<
84 fTheMaxEnergyTR = CLHEP::keV * 100.; // 40.; <<
85 <<
86 fTheMinAngle = 1.e-8; // <<
87 fTheMaxAngle = 4.e-4; <<
88 <<
89 fTotBin = 50; // number of bins in log sca <<
90 fBinTR = 100; // number of bins in TR vec <<
91 fKrange = 229; <<
92 // min/max angle2 in log-vectors <<
93 <<
94 fMinThetaTR = 3.0e-9; <<
95 fMaxThetaTR = 1.0e-4; <<
96 83
97 << 84 G4XTRenergyLoss::G4XTRenergyLoss(G4LogicalVolume *anEnvelope,
98 // Proton energy vector initialization << 85 G4Material* foilMat,G4Material* gasMat,
99 fProtonEnergyVector = << 86 G4double a, G4double b,
100 new G4PhysicsLogVector(fMinProtonTkin, fMa << 87 G4int n,const G4String& processName,
101 << 88 G4ProcessType type) :
102 fXTREnergyVector = << 89 G4VDiscreteProcess(processName, type)
103 new G4PhysicsLogVector(fTheMinEnergyTR, fT << 90 // G4VContinuousProcess(processName, type)
104 << 91 {
105 fEnvelope = anEnvelope; << 92 fEnvelope = anEnvelope ;
106 << 93 // fPlateNumber = fEnvelope->GetNoDaughters() ;
107 fPlateNumber = n; << 94 fPlateNumber = n ;
108 if(verboseLevel > 0) << 95 G4cout<<"the number of TR radiator plates = "<<fPlateNumber<<G4endl ;
109 G4cout << "### G4VXTRenergyLoss: the numbe <<
110 << fPlateNumber << G4endl; <<
111 if(fPlateNumber == 0) 96 if(fPlateNumber == 0)
112 { 97 {
113 G4Exception("G4VXTRenergyLoss::G4VXTRenerg << 98 G4Exception("No plates in X-ray TR radiator") ;
114 FatalException, "No plates in <<
115 } 99 }
116 // default is XTR dEdx, not flux after radia 100 // default is XTR dEdx, not flux after radiator
117 fExitFlux = false; << 101 fExitFlux = false;
118 // default angle distribution according nume <<
119 fFastAngle = false; // no angle accordin <<
120 fAngleRadDistr = true; <<
121 fCompton = false; <<
122 <<
123 fLambda = DBL_MAX; 102 fLambda = DBL_MAX;
124 <<
125 // Mean thicknesses of plates and gas gaps 103 // Mean thicknesses of plates and gas gaps
126 fPlateThick = a; << 104
127 fGasThick = b; << 105 fPlateThick = a ;
128 fTotalDist = fPlateNumber * (fPlateThick + << 106 fGasThick = b ;
129 if(verboseLevel > 0) << 107 fTotalDist = fPlateNumber*(fPlateThick+fGasThick) ;
130 G4cout << "total radiator thickness = " << << 108 G4cout<<"total radiator thickness = "<<fTotalDist/cm<<" cm"<<G4endl ;
131 << G4endl; <<
132 109
133 // index of plate material 110 // index of plate material
134 fMatIndex1 = (G4int)foilMat->GetIndex(); << 111 fMatIndex1 = foilMat->GetIndex() ;
135 if(verboseLevel > 0) << 112 G4cout<<"plate material = "<<foilMat->GetName()<<G4endl ;
136 G4cout << "plate material = " << foilMat-> <<
137 113
138 // index of gas material 114 // index of gas material
139 fMatIndex2 = (G4int)gasMat->GetIndex(); << 115 fMatIndex2 = gasMat->GetIndex() ;
140 if(verboseLevel > 0) << 116 G4cout<<"gas material = "<<gasMat->GetName()<<G4endl ;
141 G4cout << "gas material = " << gasMat->Get <<
142 117
143 // plasma energy squared for plate material 118 // plasma energy squared for plate material
144 fSigma1 = fPlasmaCof * foilMat->GetElectronD << 119
145 if(verboseLevel > 0) << 120 fSigma1 = fPlasmaCof*foilMat->GetElectronDensity() ;
146 G4cout << "plate plasma energy = " << std: << 121 // fSigma1 = (20.9*eV)*(20.9*eV) ;
147 << G4endl; << 122 G4cout<<"plate plasma energy = "<<sqrt(fSigma1)/eV<<" eV"<<G4endl ;
148 123
149 // plasma energy squared for gas material 124 // plasma energy squared for gas material
150 fSigma2 = fPlasmaCof * gasMat->GetElectronDe << 125
151 if(verboseLevel > 0) << 126 fSigma2 = fPlasmaCof*gasMat->GetElectronDensity() ;
152 G4cout << "gas plasma energy = " << std::s << 127 G4cout<<"gas plasma energy = "<<sqrt(fSigma2)/eV<<" eV"<<G4endl ;
153 << G4endl; <<
154 128
155 // Compute cofs for preparation of linear ph 129 // Compute cofs for preparation of linear photo absorption
156 ComputePlatePhotoAbsCof(); << 130
157 ComputeGasPhotoAbsCof(); << 131 ComputePlatePhotoAbsCof() ;
>> 132 ComputeGasPhotoAbsCof() ;
158 133
159 pParticleChange = &fParticleChange; 134 pParticleChange = &fParticleChange;
>> 135
160 } 136 }
161 137
162 ////////////////////////////////////////////// 138 ///////////////////////////////////////////////////////////////////////////
163 G4VXTRenergyLoss::~G4VXTRenergyLoss() <<
164 { <<
165 delete fProtonEnergyVector; <<
166 delete fXTREnergyVector; <<
167 if(fEnergyDistrTable) <<
168 { <<
169 fEnergyDistrTable->clearAndDestroy(); <<
170 delete fEnergyDistrTable; <<
171 } <<
172 if(fAngleRadDistr) <<
173 { <<
174 fAngleDistrTable->clearAndDestroy(); <<
175 delete fAngleDistrTable; <<
176 } <<
177 if(fAngleForEnergyTable) <<
178 { <<
179 fAngleForEnergyTable->clearAndDestroy(); <<
180 delete fAngleForEnergyTable; <<
181 } <<
182 } <<
183 139
184 void G4VXTRenergyLoss::ProcessDescription(std: << 140 G4XTRenergyLoss::~G4XTRenergyLoss()
185 { 141 {
186 out << "Base class for 'fast' parameterisati << 142 G4int i ;
187 "transition\n" << 143
188 "radiation. Angular distribution is v << 144 if(fEnvelope) delete fEnvelope;
>> 145
>> 146 for(i=0;i<fGasIntervalNumber;i++)
>> 147 {
>> 148 delete[] fGasPhotoAbsCof[i] ;
>> 149 }
>> 150 delete[] fGasPhotoAbsCof ;
>> 151
>> 152 for(i=0;i<fPlateIntervalNumber;i++)
>> 153 {
>> 154 delete[] fPlatePhotoAbsCof[i] ;
>> 155 }
>> 156 delete[] fPlatePhotoAbsCof ;
189 } 157 }
190 158
191 ////////////////////////////////////////////// 159 ///////////////////////////////////////////////////////////////////////////////
>> 160 //
192 // Returns condition for application of the mo 161 // Returns condition for application of the model depending on particle type
193 G4bool G4VXTRenergyLoss::IsApplicable(const G4 << 162
>> 163
>> 164 G4bool G4XTRenergyLoss::IsApplicable(const G4ParticleDefinition& particle)
>> 165 {
>> 166 return ( particle.GetPDGCharge() != 0.0 ) ;
>> 167 }
>> 168
>> 169 //////////////////////////////////////////////////////////////////////////////////
>> 170 //
>> 171 // GetContinuousStepLimit
>> 172 //
>> 173
>> 174 G4double
>> 175 G4XTRenergyLoss::GetContinuousStepLimit(const G4Track& ,
>> 176 G4double ,
>> 177 G4double ,
>> 178 G4double& )
194 { 179 {
195 return (particle.GetPDGCharge() != 0.0); << 180 G4double StepLimit = DBL_MAX;
>> 181
>> 182 return StepLimit;
196 } 183 }
197 184
198 ////////////////////////////////////////////// 185 /////////////////////////////////////////////////////////////////////////////////
>> 186 //
199 // Calculate step size for XTR process inside 187 // Calculate step size for XTR process inside raaditor
200 G4double G4VXTRenergyLoss::GetMeanFreePath(con << 188
201 G4F << 189 G4double G4XTRenergyLoss::GetMeanFreePath(const G4Track& aTrack,
>> 190 G4double, // previousStepSize,
>> 191 G4ForceCondition* condition)
202 { 192 {
203 G4int iTkin, iPlace; 193 G4int iTkin, iPlace;
204 G4double lambda, sigma, kinEnergy, mass, gam 194 G4double lambda, sigma, kinEnergy, mass, gamma;
205 G4double charge, chargeSq, massRatio, TkinSc 195 G4double charge, chargeSq, massRatio, TkinScaled;
206 G4double E1, E2, W, W1, W2; << 196 G4double E1,E2,W,W1,W2;
207 <<
208 *condition = NotForced; <<
209 197
210 if(aTrack.GetVolume()->GetLogicalVolume() != << 198 *condition = NotForced;
211 lambda = DBL_MAX; << 199
>> 200 if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope ) lambda = DBL_MAX;
212 else 201 else
213 { 202 {
214 const G4DynamicParticle* aParticle = aTrac 203 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
215 kinEnergy = aPart << 204 kinEnergy = aParticle->GetKineticEnergy();
216 mass = aParticle->GetDefinition()->GetPDG << 205 mass = aParticle->GetDefinition()->GetPDGMass();
217 gamma = 1.0 + kinEnergy / mass; << 206 gamma = 1.0 + kinEnergy/mass;
218 if(verboseLevel > 1) << 207 if(verboseLevel)
219 { 208 {
220 G4cout << " gamma = " << gamma << "; f << 209 G4cout<<" gamma = "<<gamma<<"; fGamma = "<<fGamma<<G4endl;
221 } 210 }
222 211
223 if(std::fabs(gamma - fGamma) < 0.05 * gamm << 212 if ( fabs( gamma - fGamma ) < 0.05*gamma ) lambda = fLambda;
224 lambda = fLambda; <<
225 else 213 else
226 { 214 {
227 charge = aParticle->GetDefinition()- << 215 charge = aParticle->GetDefinition()->GetPDGCharge();
228 chargeSq = charge * charge; << 216 chargeSq = charge*charge;
229 massRatio = proton_mass_c2 / mass; << 217 massRatio = proton_mass_c2/mass;
230 TkinScaled = kinEnergy * massRatio; << 218 TkinScaled = kinEnergy*massRatio;
231 219
232 for(iTkin = 0; iTkin < fTotBin; ++iTkin) << 220 for(iTkin = 0; iTkin < fTotBin; iTkin++)
233 { 221 {
234 if(TkinScaled < fProtonEnergyVector->G << 222 if( TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;
235 break; <<
236 } 223 }
237 iPlace = iTkin - 1; << 224 iPlace = iTkin - 1 ;
238 225
239 if(iTkin == 0) << 226 if(iTkin == 0) lambda = DBL_MAX; // Tkin is too small, neglect of TR photon generation
240 lambda = DBL_MAX; // Tkin is too smal << 227 else // general case: Tkin between two vectors of the material
241 else // general case: Tkin between two <<
242 { 228 {
243 if(iTkin == fTotBin) << 229 if(iTkin == fTotBin)
244 { 230 {
245 sigma = (*(*fEnergyDistrTable)(iPlac << 231 sigma = (*(*fEnergyDistrTable)(iPlace))(0)*chargeSq;
246 } 232 }
247 else 233 else
248 { 234 {
249 E1 = fProtonEnergyVector->GetLowE << 235 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
250 E2 = fProtonEnergyVector->GetLowE << 236 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
251 W = 1.0 / (E2 - E1); << 237 W = 1.0/(E2 - E1) ;
252 W1 = (E2 - TkinScaled) * W; << 238 W1 = (E2 - TkinScaled)*W ;
253 W2 = (TkinScaled - E1) * W; << 239 W2 = (TkinScaled - E1)*W ;
254 sigma = ((*(*fEnergyDistrTable)(iPla << 240 sigma = ( (*(*fEnergyDistrTable)(iPlace ))(0)*W1 +
255 (*(*fEnergyDistrTable)(iPla << 241 (*(*fEnergyDistrTable)(iPlace+1))(0)*W2 )*chargeSq;
256 chargeSq; << 242
257 } 243 }
258 if(sigma < DBL_MIN) << 244 if (sigma < DBL_MIN) lambda = DBL_MAX;
259 lambda = DBL_MAX; << 245 else lambda = 1./sigma;
260 else <<
261 lambda = 1. / sigma; <<
262 fLambda = lambda; 246 fLambda = lambda;
263 fGamma = gamma; << 247 fGamma = gamma;
264 if(verboseLevel > 1) << 248 if(verboseLevel)
265 { 249 {
266 G4cout << " lambda = " << lambda / m << 250 G4cout<<" lambda = "<<lambda/mm<<" mm"<<G4endl;
267 } 251 }
268 } 252 }
269 } 253 }
270 } << 254 }
271 return lambda; 255 return lambda;
272 } 256 }
273 257
274 ////////////////////////////////////////////// <<
275 // Interface for build table from physics list <<
276 void G4VXTRenergyLoss::BuildPhysicsTable(const <<
277 { <<
278 if(pd.GetPDGCharge() == 0.) <<
279 { <<
280 G4Exception("G4VXTRenergyLoss::BuildPhysic <<
281 JustWarning, "XTR initialisati <<
282 } <<
283 BuildEnergyTable(); <<
284 <<
285 if(fAngleRadDistr) <<
286 { <<
287 if(verboseLevel > 0) <<
288 { <<
289 G4cout <<
290 << "Build angle for energy distributio <<
291 << G4endl; <<
292 } <<
293 BuildAngleForEnergyBank(); <<
294 } <<
295 } <<
296 258
297 ////////////////////////////////////////////// 259 //////////////////////////////////////////////////////////////////////////
>> 260 //
298 // Build integral energy distribution of XTR p 261 // Build integral energy distribution of XTR photons
299 void G4VXTRenergyLoss::BuildEnergyTable() << 262
>> 263 void G4XTRenergyLoss::BuildTable()
300 { 264 {
301 G4int iTkin, iTR, iPlace; 265 G4int iTkin, iTR, iPlace;
302 G4double radiatorCof = 1.0; // for tuning o << 266 G4double radiatorCof = 1.0; // for tuning of XTR yield
303 G4double energySum = 0.0; <<
304 267
305 fEnergyDistrTable = new G4PhysicsTable(fTotB 268 fEnergyDistrTable = new G4PhysicsTable(fTotBin);
306 if(fAngleRadDistr) <<
307 fAngleDistrTable = new G4PhysicsTable(fTot <<
308 269
309 fGammaTkinCut = 0.0; 270 fGammaTkinCut = 0.0;
>> 271
>> 272 // setting of min/max TR energies
>> 273
>> 274 if(fGammaTkinCut > fTheMinEnergyTR) fMinEnergyTR = fGammaTkinCut ;
>> 275 else fMinEnergyTR = fTheMinEnergyTR ;
>> 276
>> 277 if(fGammaTkinCut > fTheMaxEnergyTR) fMaxEnergyTR = 2.0*fGammaTkinCut ;
>> 278 else fMaxEnergyTR = fTheMaxEnergyTR ;
310 279
311 // setting of min/max TR energies << 280 G4cout.precision(4) ;
312 if(fGammaTkinCut > fTheMinEnergyTR) << 281 G4Timer timer ;
313 fMinEnergyTR = fGammaTkinCut; << 282 timer.Start() ;
314 else << 283 G4cout<<G4endl;
315 fMinEnergyTR = fTheMinEnergyTR; << 284 G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
316 << 285 G4cout<<G4endl;
317 if(fGammaTkinCut > fTheMaxEnergyTR) << 286
318 fMaxEnergyTR = 2.0 * fGammaTkinCut; << 287 for( iTkin = 0 ; iTkin < fTotBin ; iTkin++ ) // Lorentz factor loop
319 else <<
320 fMaxEnergyTR = fTheMaxEnergyTR; <<
321 <<
322 G4Integrator<G4VXTRenergyLoss, G4double (G4V <<
323 integral; <<
324 <<
325 G4cout.precision(4); <<
326 G4Timer timer; <<
327 timer.Start(); <<
328 <<
329 if(verboseLevel > 0) <<
330 { <<
331 G4cout << G4endl; <<
332 G4cout << "Lorentz Factor" <<
333 << "\t" <<
334 << "XTR photon number" << G4endl; <<
335 G4cout << G4endl; <<
336 } <<
337 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // <<
338 { 288 {
339 auto energyVector = << 289 G4PhysicsLogVector* energyVector = new G4PhysicsLogVector( fMinEnergyTR,
340 new G4PhysicsLogVector(fMinEnergyTR, fMa << 290 fMaxEnergyTR,
341 << 291 fBinTR ) ;
342 fGamma = <<
343 1.0 + (fProtonEnergyVector->GetLowEdgeEn <<
344 <<
345 // if(fMaxThetaTR > fTheMaxAngle) fMax <<
346 // else if(fMaxThetaTR < fTheMinAngle) <<
347 292
348 energySum = 0.0; << 293 fGamma = 1.0 + (fProtonEnergyVector->
>> 294 GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
349 295
350 energyVector->PutValue(fBinTR - 1, energyS << 296 fMaxThetaTR = 25.0/(fGamma*fGamma) ; // theta^2
351 297
352 for(iTR = fBinTR - 2; iTR >= 0; --iTR) << 298 fTheMinAngle = 1.0e-3 ; // was 5.e-6, e-6 !!!, e-5, e-4
353 { << 299
354 // Legendre96 or Legendre10 << 300 if( fMaxThetaTR > fTheMaxAngle ) fMaxThetaTR = fTheMaxAngle;
355 << 301 else
356 energySum += radiatorCof * fCofTR * << 302 {
357 << 303 if( fMaxThetaTR < fTheMinAngle ) fMaxThetaTR = fTheMinAngle;
358 // integral.Legendre10(this, &G4VXTRenergyLo << 304 }
359 << 305 G4PhysicsLinearVector* angleVector = new G4PhysicsLinearVector(0.0,
360 integral.Legendre96(this, & << 306 fMaxThetaTR,
361 << 307 fBinTR );
362 energyV << 308
363 energyV << 309 G4double energySum = 0.0;
364 << 310 G4double angleSum = 0.0;
365 energyVector->PutValue(iTR, energySum / << 311
366 } << 312 G4Integrator<G4XTRenergyLoss,G4double(G4XTRenergyLoss::*)(G4double)> integral;
367 iPlace = iTkin; << 313
368 fEnergyDistrTable->insertAt(iPlace, energy << 314 energyVector->PutValue(fBinTR-1,energySum);
369 << 315 angleVector->PutValue(fBinTR-1,angleSum);
370 if(verboseLevel > 0) << 316
371 { << 317 for( iTR = fBinTR - 2 ; iTR >= 0 ; iTR-- )
372 G4cout << fGamma << "\t" << energySum << << 318 {
373 } << 319 energySum += radiatorCof*fCofTR*integral.Legendre10(
374 } << 320 this,&G4XTRenergyLoss::SpectralXTRdEdx,
>> 321 energyVector->GetLowEdgeEnergy(iTR),
>> 322 energyVector->GetLowEdgeEnergy(iTR+1) );
>> 323
>> 324 // angleSum += fCofTR*integral.Legendre96(
>> 325 // this,&G4XTRenergyLoss::AngleXTRdEdx,
>> 326 // angleVector->GetLowEdgeEnergy(iTR),
>> 327 // angleVector->GetLowEdgeEnergy(iTR+1) );
>> 328
>> 329 energyVector->PutValue(iTR,energySum/fTotalDist);
>> 330 // angleVector ->PutValue(iTR,angleSum);
>> 331 }
>> 332 G4cout
>> 333 // <<iTkin<<"\t"
>> 334 // <<"fGamma = "
>> 335 <<fGamma<<"\t" // <<" fMaxThetaTR = "<<fMaxThetaTR
>> 336 // <<"sumN = "
>> 337 <<energySum // <<" ; sumA = "<<angleSum
>> 338 <<G4endl;
>> 339 iPlace = iTkin;
>> 340 fEnergyDistrTable->insertAt(iPlace,energyVector);
>> 341 // fAngleDistrTable->insertAt(iPlace,angleVector);
>> 342 }
375 timer.Stop(); 343 timer.Stop();
376 G4cout.precision(6); 344 G4cout.precision(6);
377 if(verboseLevel > 0) << 345 G4cout<<G4endl;
378 { << 346 G4cout<<"total time for build X-ray TR energy loss tables = "
379 G4cout << G4endl; << 347 <<timer.GetUserElapsed()<<" s"<<G4endl;
380 G4cout << "total time for build X-ray TR e <<
381 << timer.GetUserElapsed() << " s" < <<
382 } <<
383 fGamma = 0.; 348 fGamma = 0.;
384 return; << 349 return ;
385 } 350 }
386 351
387 ////////////////////////////////////////////// 352 //////////////////////////////////////////////////////////////////////////
388 // Bank of angle distributions for given energ << 353 //
>> 354 //
389 355
390 void G4VXTRenergyLoss::BuildAngleForEnergyBank << 356 void G4XTRenergyLoss::BuildEnergyTable()
391 { 357 {
392 << 358 return ;
393 if( ( this->GetProcessName() == "TranspRegXT << 359 }
394 this->GetProcessName() == "TranspRegXT <<
395 this->GetProcessName() == "RegularXTRa <<
396 this->GetProcessName() == "RegularXTRmodel" <<
397 { <<
398 BuildAngleTable(); // by sum of delta-func <<
399 return; <<
400 } <<
401 G4int i, iTkin, iTR; <<
402 G4double angleSum = 0.0; <<
403 <<
404 fGammaTkinCut = 0.0; <<
405 <<
406 // setting of min/max TR energies <<
407 if(fGammaTkinCut > fTheMinEnergyTR) <<
408 fMinEnergyTR = fGammaTkinCut; <<
409 else <<
410 fMinEnergyTR = fTheMinEnergyTR; <<
411 <<
412 if(fGammaTkinCut > fTheMaxEnergyTR) <<
413 fMaxEnergyTR = 2.0 * fGammaTkinCut; <<
414 else <<
415 fMaxEnergyTR = fTheMaxEnergyTR; <<
416 <<
417 auto energyVector = <<
418 new G4PhysicsLogVector(fMinEnergyTR, fMaxE <<
419 <<
420 G4Integrator<G4VXTRenergyLoss, G4double (G4V <<
421 integral; <<
422 360
423 G4cout.precision(4); << 361 ////////////////////////////////////////////////////////////////////////
424 G4Timer timer; << 362 //
425 timer.Start(); << 363 //
426 364
427 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // << 365 void G4XTRenergyLoss::BuildAngleTable()
428 { << 366 {
429 fGamma = << 367 return ;
430 1.0 + (fProtonEnergyVector->GetLowEdgeEn << 368 }
431 369
432 if(fMaxThetaTR > fTheMaxAngle) <<
433 fMaxThetaTR = fTheMaxAngle; <<
434 else if(fMaxThetaTR < fTheMinAngle) <<
435 fMaxThetaTR = fTheMinAngle; <<
436 370
437 fAngleForEnergyTable = new G4PhysicsTable( << 371 //////////////////////////////////////////////////////////////////////////////
>> 372 //
>> 373 // The main function which is responsible for the treatment of a particle passage
>> 374 // trough G4Envelope with discrete generation of G4Gamma
438 375
439 for(iTR = 0; iTR < fBinTR; ++iTR) << 376 G4VParticleChange* G4XTRenergyLoss::PostStepDoIt( const G4Track& aTrack,
440 { << 377 const G4Step& aStep )
441 angleSum = 0.0; << 378 {
442 fEnergy = energyVector->GetLowEdgeEnerg << 379 G4int iTkin, iPlace;
443 << 380 G4double energyTR, theta, phi, dirX, dirY, dirZ;
444 // log-vector to increase number of thin <<
445 auto angleVector = new G4PhysicsLogVecto <<
446 381
447 <<
448 <<
449 angleVector->PutValue(fBinTR - 1, angleS <<
450 382
451 for(i = fBinTR - 2; i >= 0; --i) << 383 fParticleChange.Initialize(aTrack);
452 { <<
453 // Legendre96 or Legendre10 <<
454 <<
455 angleSum += <<
456 integral.Legendre10(this, &G4VXTRene <<
457 angleVector->Get <<
458 angleVector->Get <<
459 384
460 angleVector->PutValue(i, angleSum); << 385 if(verboseLevel)
461 } <<
462 fAngleForEnergyTable->insertAt(iTR, angl <<
463 } <<
464 fAngleBank.push_back(fAngleForEnergyTable) <<
465 } <<
466 timer.Stop(); <<
467 G4cout.precision(6); <<
468 if(verboseLevel > 0) <<
469 { 386 {
470 G4cout << G4endl; << 387 G4cout<<"Start of G4XTRenergyLoss::PostStepDoIt "<<G4endl ;
471 G4cout << "total time for build X-ray TR a << 388 G4cout<<"name of current material = "
472 << timer.GetUserElapsed() << " s" < << 389 <<aTrack.GetVolume()->GetLogicalVolume()->GetMaterial()->GetName()<<G4endl ;
473 } 390 }
474 fGamma = 0.; << 391 if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope )
475 delete energyVector; <<
476 } <<
477 <<
478 ////////////////////////////////////////////// <<
479 // Build XTR angular distribution at given ene <<
480 // of transparent regular radiator <<
481 void G4VXTRenergyLoss::BuildAngleTable() <<
482 { <<
483 G4int iTkin, iTR; <<
484 G4double energy; <<
485 <<
486 fGammaTkinCut = 0.0; <<
487 <<
488 // setting of min/max TR energies <<
489 if(fGammaTkinCut > fTheMinEnergyTR) <<
490 fMinEnergyTR = fGammaTkinCut; <<
491 else <<
492 fMinEnergyTR = fTheMinEnergyTR; <<
493 <<
494 if(fGammaTkinCut > fTheMaxEnergyTR) <<
495 fMaxEnergyTR = 2.0 * fGammaTkinCut; <<
496 else <<
497 fMaxEnergyTR = fTheMaxEnergyTR; <<
498 <<
499 G4cout.precision(4); <<
500 G4Timer timer; <<
501 timer.Start(); <<
502 if(verboseLevel > 0) <<
503 { 392 {
504 G4cout << G4endl << "Lorentz Factor" << "\ << 393 if(verboseLevel)
505 << "XTR photon number" << G4endl << << 394 {
>> 395 G4cout<<"Go out from G4XTRenergyLoss::PostStepDoIt: wrong volume "<<G4endl;
>> 396 }
>> 397 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
506 } 398 }
507 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // << 399 else
508 { 400 {
509 fGamma = << 401 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
510 1.0 + (fProtonEnergyVector->GetLowEdgeEn << 402 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
>> 403
>> 404 // Now we are ready to Generate one TR photon
511 405
512 // fMaxThetaTR = 25. * 2500.0 / (fGamma * << 406 G4double kinEnergy = aParticle->GetKineticEnergy() ;
>> 407 G4double mass = aParticle->GetDefinition()->GetPDGMass() ;
>> 408 G4double gamma = 1.0 + kinEnergy/mass ;
513 409
514 if(fMaxThetaTR > fTheMaxAngle) << 410 if(verboseLevel > 0 )
515 fMaxThetaTR = fTheMaxAngle; <<
516 else <<
517 { 411 {
518 if(fMaxThetaTR < fTheMinAngle) << 412 G4cout<<"gamma = "<<gamma<<G4endl ;
519 fMaxThetaTR = fTheMinAngle; <<
520 } 413 }
>> 414 G4double massRatio = proton_mass_c2/mass ;
>> 415 G4double TkinScaled = kinEnergy*massRatio ;
>> 416 G4ThreeVector position = pPostStepPoint->GetPosition();
>> 417 G4ParticleMomentum direction = aParticle->GetMomentumDirection();
>> 418 G4double startTime = pPostStepPoint->GetGlobalTime();
521 419
522 fAngleForEnergyTable = new G4PhysicsTable( << 420 for( iTkin = 0; iTkin < fTotBin; iTkin++ )
>> 421 {
>> 422 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
>> 423 }
>> 424 iPlace = iTkin - 1;
523 425
524 for(iTR = 0; iTR < fBinTR; ++iTR) << 426 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
525 { 427 {
526 energy = fXTREnergyVector->GetLowEdgeEne << 428 if( verboseLevel )
>> 429 {
>> 430 G4cout<<"Go out from G4XTRenergyLoss::PostStepDoIt:iTkin = "<<iTkin<<G4endl;
>> 431 }
>> 432 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
>> 433 }
>> 434 else // general case: Tkin between two vectors of the material
>> 435 {
>> 436 fParticleChange.SetNumberOfSecondaries(1);
527 437
528 G4PhysicsFreeVector* angleVector = GetAn << 438 energyTR = GetXTRrandomEnergy(TkinScaled,iTkin);
529 439
530 fAngleForEnergyTable->insertAt(iTR, angl << 440 if( verboseLevel )
531 } << 441 {
532 fAngleBank.push_back(fAngleForEnergyTable) << 442 G4cout<<"energyTR = "<<energyTR/keV<<" keV"<<G4endl;
533 } << 443 }
534 timer.Stop(); << 444 theta = fabs(G4RandGauss::shoot(0.0,pi/gamma));
535 G4cout.precision(6); <<
536 if(verboseLevel > 0) <<
537 { <<
538 G4cout << G4endl; <<
539 G4cout << "total time for build XTR angle <<
540 << timer.GetUserElapsed() << " s" < <<
541 } <<
542 fGamma = 0.; <<
543 445
544 return; << 446 if( theta >= 0.1 ) theta = 0.1;
545 } <<
546 447
547 ////////////////////////////////////////////// << 448 // G4cout<<" : theta = "<<theta<<endl ;
548 // Vector of angles and angle integral distrib <<
549 G4PhysicsFreeVector* G4VXTRenergyLoss::GetAngl <<
550 { <<
551 G4double theta = 0., result, tmp = 0., cof1, <<
552 angleSum = 0.; <<
553 G4int iTheta, k, kMin; <<
554 449
555 auto angleVector = new G4PhysicsFreeVector(n << 450 phi = twopi*G4UniformRand();
556 451
557 cofPHC = 4. * pi * hbarc; << 452 dirX = sin(theta)*cos(phi);
558 tmp = (fSigma1 - fSigma2) / cofPHC / ener << 453 dirY = sin(theta)*sin(phi);
559 cof1 = fPlateThick * tmp; << 454 dirZ = cos(theta);
560 cof2 = fGasThick * tmp; <<
561 455
562 cofMin = energy * (fPlateThick + fGasThick) << 456 G4ThreeVector directionTR(dirX,dirY,dirZ);
563 cofMin += (fPlateThick * fSigma1 + fGasThick << 457 directionTR.rotateUz(direction);
564 cofMin /= cofPHC; << 458 directionTR.unit();
565 459
566 kMin = G4int(cofMin); << 460 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
567 if(cofMin > kMin) << 461 directionTR, energyTR);
568 kMin++; <<
569 462
570 if(verboseLevel > 2) << 463 // A XTR photon is set on the particle track inside the radiator
571 { << 464 // and is moved to the G4Envelope surface for standard X-ray TR models
572 G4cout << "n-1 = " << n - 1 << 465 // only. The case of fExitFlux=true
573 << "; theta = " << std::sqrt(fMaxTh <<
574 << "; tmp = " << 0. << "; angleS <<
575 } <<
576 466
577 for(iTheta = n - 1; iTheta >= 1; --iTheta) << 467 if( fExitFlux )
578 { << 468 {
579 k = iTheta - 1 + kMin; << 469 const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
580 tmp = pi * fPlateThick * (k + cof2) / ( << 470 G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
581 result = (k - cof1) * (k - cof1) * (k + co << 471 G4AffineTransform transform = G4AffineTransform(rotM,transl);
582 tmp = std::sin(tmp) * std::sin(tmp) * s << 472 transform.Invert();
>> 473 G4ThreeVector localP = transform.TransformPoint(position);
>> 474 G4ThreeVector localV = transform.TransformAxis(directionTR);
583 475
584 if(k == kMin && kMin == G4int(cofMin)) << 476 G4double distance = fEnvelope->GetSolid()->DistanceToOut(localP, localV);
585 { << 477 if(verboseLevel)
586 // angleSum += 0.5 * tmp; << 478 {
587 angleSum += tmp; // ATLAS TB << 479 G4cout<<"distance to exit = "<<distance/mm<<" mm"<<G4endl;
588 } << 480 }
589 else if(iTheta == n - 1) << 481 position += distance*directionTR;
590 ; << 482 startTime += distance/c_light;
591 else << 483 }
592 { << 484 G4Track* aSecondaryTrack = new G4Track( aPhotonTR,
593 angleSum += tmp; << 485 startTime, position );
594 } << 486 aSecondaryTrack->SetTouchableHandle(
595 theta = std::abs(k - cofMin) * cofPHC / en << 487 aStep.GetPostStepPoint()->GetTouchableHandle());
>> 488 aSecondaryTrack->SetParentID( aTrack.GetTrackID() );
596 489
597 if(verboseLevel > 2) << 490 fParticleChange.AddSecondary(aSecondaryTrack);
598 { << 491 fParticleChange.ProposeEnergy(kinEnergy);
599 G4cout << "iTheta = " << iTheta << "; k <<
600 << "; theta = " << std::sqrt(thet <<
601 << "; angleSum = " << angleSum <<
602 } 492 }
603 angleVector->PutValue(iTheta, theta, angle <<
604 } <<
605 if(theta > 0.) <<
606 { <<
607 // angleSum += 0.5 * tmp; <<
608 angleSum += 0.; // ATLAS TB <<
609 theta = 0.; <<
610 } 493 }
611 if(verboseLevel > 2) << 494 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
612 { <<
613 G4cout << "iTheta = " << iTheta << "; thet <<
614 << "; tmp = " << tmp << "; angle <<
615 } <<
616 angleVector->PutValue(iTheta, theta, angleSu <<
617 <<
618 return angleVector; <<
619 } 495 }
620 496
621 ////////////////////////////////////////////// <<
622 // Build XTR angular distribution based on the <<
623 // radiator <<
624 void G4VXTRenergyLoss::BuildGlobalAngleTable() <<
625 { <<
626 G4int iTkin, iTR, iPlace; <<
627 G4double radiatorCof = 1.0; // for tuning o <<
628 G4double angleSum; <<
629 fAngleDistrTable = new G4PhysicsTable(fTotBi <<
630 497
631 fGammaTkinCut = 0.0; <<
632 498
633 // setting of min/max TR energies << 499 //////////////////////////////////////////////////////////////////////////////
634 if(fGammaTkinCut > fTheMinEnergyTR) << 500 //
635 fMinEnergyTR = fGammaTkinCut; << 501 // The main function which is responsible for the treatment of a particle passage
636 else << 502 // trough G4Envelope
637 fMinEnergyTR = fTheMinEnergyTR; <<
638 503
639 if(fGammaTkinCut > fTheMaxEnergyTR) << 504 G4VParticleChange* G4XTRenergyLoss::AlongStepDoIt( const G4Track& aTrack,
640 fMaxEnergyTR = 2.0 * fGammaTkinCut; << 505 const G4Step& aStep )
641 else << 506 {
642 fMaxEnergyTR = fTheMaxEnergyTR; << 507 G4int iTkin, iPlace, numOfTR, iTR ;
>> 508 G4double energyTR, meanNumOfTR, theta, phi, dirX, dirY, dirZ, rand ;
>> 509 G4double W, W1, W2, E1, E2 ;
643 510
644 G4cout.precision(4); << 511 aParticleChange.Initialize(aTrack);
645 G4Timer timer; <<
646 timer.Start(); <<
647 if(verboseLevel > 0) <<
648 { <<
649 G4cout << G4endl; <<
650 G4cout << "Lorentz Factor" <<
651 << "\t" <<
652 << "XTR photon number" << G4endl; <<
653 G4cout << G4endl; <<
654 } <<
655 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // <<
656 { <<
657 fGamma = <<
658 1.0 + (fProtonEnergyVector->GetLowEdgeEn <<
659 512
660 // fMaxThetaTR = 25.0 / (fGamma * fGamma); << 513 if(verboseLevel)
661 // fMaxThetaTR = 1.e-4; // theta^2 << 514 {
>> 515 G4cout<<"Start of G4XTRenergyLoss::AlongStepDoIt "<<G4endl ;
>> 516 G4cout<<"name of current material = "
>> 517 <<aTrack.GetVolume()->GetLogicalVolume()->GetMaterial()->GetName()<<G4endl ;
>> 518 }
>> 519 // if(aStep.GetPreStepPoint()->GetPhysicalVolume()->GetLogicalVolume() != fEnvelope)
662 520
663 if(fMaxThetaTR > fTheMaxAngle) << 521 if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope )
664 fMaxThetaTR = fTheMaxAngle; << 522 {
665 else << 523 if(verboseLevel)
666 { 524 {
667 if(fMaxThetaTR < fTheMinAngle) << 525 G4cout<<"Go out from G4XTRenergyLoss::AlongStepDoIt: wrong volume "<<G4endl;
668 fMaxThetaTR = fTheMinAngle; <<
669 } 526 }
670 auto angleVector = << 527 // return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);
671 // G4PhysicsLogVector* angleVector = << 528 }
672 new G4PhysicsLinearVector(0.0, fMaxTheta << 529 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
673 // new G4PhysicsLogVector(1.e-8, fMaxThet << 530 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
674 << 531
675 angleSum = 0.0; << 532 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
676 << 533 G4double charge = aParticle->GetDefinition()->GetPDGCharge();
677 G4Integrator<G4VXTRenergyLoss, G4double (G << 534
678 integral; << 535
679 << 536 // Now we are ready to Generate TR photons
680 angleVector->PutValue(fBinTR - 1, angleSum <<
681 537
682 for(iTR = fBinTR - 2; iTR >= 0; --iTR) << 538 G4double chargeSq = charge*charge ;
683 { << 539 G4double kinEnergy = aParticle->GetKineticEnergy() ;
684 angleSum += radiatorCof * fCofTR * << 540 G4double mass = aParticle->GetDefinition()->GetPDGMass() ;
685 integral.Legendre96(this, &G << 541 G4double gamma = 1.0 + kinEnergy/mass ;
686 angleVec <<
687 angleVec <<
688 542
689 angleVector->PutValue(iTR, angleSum); << 543 if(verboseLevel > 0 )
690 } <<
691 if(verboseLevel > 1) <<
692 { <<
693 G4cout << fGamma << "\t" << angleSum << <<
694 } <<
695 iPlace = iTkin; <<
696 fAngleDistrTable->insertAt(iPlace, angleVe <<
697 } <<
698 timer.Stop(); <<
699 G4cout.precision(6); <<
700 if(verboseLevel > 0) <<
701 { 544 {
702 G4cout << G4endl; << 545 G4cout<<"gamma = "<<gamma<<G4endl ;
703 G4cout << "total time for build X-ray TR a <<
704 << timer.GetUserElapsed() << " s" < <<
705 } 546 }
706 fGamma = 0.; << 547 G4double massRatio = proton_mass_c2/mass ;
>> 548 G4double TkinScaled = kinEnergy*massRatio ;
707 549
708 return; << 550 G4ThreeVector startPos = pPreStepPoint->GetPosition();
709 } << 551 G4double startTime = pPreStepPoint->GetGlobalTime();
710 552
711 ////////////////////////////////////////////// << 553 G4ParticleMomentum direction = aParticle->GetMomentumDirection();
712 // The main function which is responsible for << 554
713 // passage through G4Envelope with discrete ge << 555 G4double distance = aStep.GetStepLength() ;
714 G4VParticleChange* G4VXTRenergyLoss::PostStepD <<
715 <<
716 { <<
717 G4int iTkin; <<
718 G4double energyTR, theta, theta2, phi, dirX, <<
719 556
720 fParticleChange.Initialize(aTrack); <<
721 557
722 if(verboseLevel > 1) << 558 for(iTkin=0;iTkin<fTotBin;iTkin++)
723 { 559 {
724 G4cout << "Start of G4VXTRenergyLoss::Post << 560 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break;
725 G4cout << "name of current material = " <<
726 << aTrack.GetVolume()->GetLogicalVo <<
727 << G4endl; <<
728 } 561 }
729 if(aTrack.GetVolume()->GetLogicalVolume() != << 562 iPlace = iTkin - 1 ;
>> 563
>> 564 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
730 { 565 {
731 if(verboseLevel > 0) << 566 if(verboseLevel)
732 { 567 {
733 G4cout << "Go out from G4VXTRenergyLoss: << 568 G4cout<<"Go out from G4XTRenergyLoss::AlongStepDoIt:iTkin = "<<iTkin<<G4endl;
734 << G4endl; <<
735 } 569 }
736 return G4VDiscreteProcess::PostStepDoIt(aT << 570 // return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);
737 } << 571 }
738 else << 572 else // general case: Tkin between two vectors of the material
739 { 573 {
740 G4StepPoint* pPostStepPoint = aStep << 574 if(iTkin == fTotBin)
741 const G4DynamicParticle* aParticle = aTrac <<
742 <<
743 // Now we are ready to Generate one TR pho <<
744 G4double kinEnergy = aParticle->GetKinetic <<
745 G4double mass = aParticle->GetDefinit <<
746 G4double gamma = 1.0 + kinEnergy / mas <<
747 <<
748 if(verboseLevel > 1) <<
749 { 575 {
750 G4cout << "gamma = " << gamma << G4endl; << 576 meanNumOfTR = (*(*fEnergyDistrTable)(iPlace))(0)*chargeSq*distance ;
>> 577 numOfTR = G4Poisson(meanNumOfTR) ;
751 } 578 }
752 G4double massRatio = proton_mass << 579 else
753 G4double TkinScaled = kinEnergy * <<
754 G4ThreeVector position = pPostStepPo <<
755 G4ParticleMomentum direction = aParticle-> <<
756 G4double startTime = pPostStepPo <<
757 <<
758 for(iTkin = 0; iTkin < fTotBin; ++iTkin) <<
759 { 580 {
760 if(TkinScaled < fProtonEnergyVector->Get << 581 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
761 break; << 582 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
>> 583 W = 1.0/(E2 - E1) ;
>> 584 W1 = (E2 - TkinScaled)*W ;
>> 585 W2 = (TkinScaled - E1)*W ;
>> 586 meanNumOfTR = ( (*(*fEnergyDistrTable)(iPlace ))(0)*W1+
>> 587 (*(*fEnergyDistrTable)(iPlace+1))(0)*W2 )*chargeSq*distance ;
>> 588
>> 589 if(verboseLevel > 0 )
>> 590 {
>> 591 G4cout<<iTkin<<" mean TR number = "<<meanNumOfTR
>> 592 <<" or mean over energy-angle tables "
>> 593 <<(((*(*fEnergyDistrTable)(iPlace))(0)+
>> 594 (*(*fAngleDistrTable)(iPlace))(0))*W1 +
>> 595 ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
>> 596 (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)*chargeSq*0.5
>> 597 <<G4endl ;
>> 598 }
>> 599 numOfTR = G4Poisson( meanNumOfTR ) ;
762 } 600 }
763 << 601 if( numOfTR == 0 ) // no change, return
764 if(iTkin == 0) // Tkin is too small, negl <<
765 { 602 {
766 if(verboseLevel > 0) << 603 aParticleChange.SetNumberOfSecondaries(0);
>> 604 if(verboseLevel)
767 { 605 {
768 G4cout << "Go out from G4VXTRenergyLos << 606 G4cout<<"Go out from G4XTRenergyLoss::AlongStepDoIt: numOfTR = "
769 << G4endl; << 607 <<numOfTR<<G4endl ;
770 } 608 }
771 return G4VDiscreteProcess::PostStepDoIt( << 609 // return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);
772 } 610 }
773 else // general case: Tkin between two ve << 611 else
774 { 612 {
775 fParticleChange.SetNumberOfSecondaries(1 << 613 if(verboseLevel)
>> 614 {
>> 615 G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl ;
>> 616 }
>> 617 aParticleChange.SetNumberOfSecondaries(numOfTR);
776 618
777 energyTR = GetXTRrandomEnergy(TkinScaled << 619 G4double sumEnergyTR = 0.0 ;
778 620
779 if(verboseLevel > 1) << 621 for(iTR=0;iTR<numOfTR;iTR++)
780 { 622 {
781 G4cout << "energyTR = " << energyTR / << 623
782 } << 624 // energyPos = ((*(*fEnergyDistrTable)(iPlace))(0)*W1+
783 if(fAngleRadDistr) << 625 // (*(*fEnergyDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand() ;
>> 626 // for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
>> 627 // {
>> 628 // if(energyPos >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1+
>> 629 // (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer)*W2)) break ;
>> 630 // }
>> 631 // energyTR = ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer))*W1+
>> 632 // ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer))*W2 ;
>> 633
>> 634 energyTR = GetXTRrandomEnergy(TkinScaled,iTkin) ;
>> 635
>> 636 if(verboseLevel)
784 { 637 {
785 theta2 = GetRandomAngle(energyTR, iTki << 638 G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl ;
786 if(theta2 > 0.) <<
787 theta = std::sqrt(theta2); <<
788 else <<
789 theta = 0.; <<
790 } 639 }
791 else << 640 sumEnergyTR += energyTR ;
792 theta = std::fabs(G4RandGauss::shoot(0 <<
793 641
794 if(theta >= 0.1) << 642 theta = fabs(G4RandGauss::shoot(0.0,pi/gamma)) ;
795 theta = 0.1; <<
796 643
797 phi = twopi * G4UniformRand(); << 644 if( theta >= 0.1 ) theta = 0.1 ;
798 645
799 dirX = std::sin(theta) * std::cos(phi); << 646 // G4cout<<" : theta = "<<theta<<endl ;
800 dirY = std::sin(theta) * std::sin(phi); <<
801 dirZ = std::cos(theta); <<
802 647
803 G4ThreeVector directionTR(dirX, dirY, di << 648 phi = twopi*G4UniformRand() ;
804 directionTR.rotateUz(direction); <<
805 directionTR.unit(); <<
806 649
807 auto aPhotonTR = << 650 dirX = sin(theta)*cos(phi) ;
808 new G4DynamicParticle(G4Gamma::Gamma() << 651 dirY = sin(theta)*sin(phi) ;
>> 652 dirZ = cos(theta) ;
809 653
810 // A XTR photon is set on the particle t << 654 G4ThreeVector directionTR(dirX,dirY,dirZ) ;
811 // and is moved to the G4Envelope surfac << 655 directionTR.rotateUz(direction) ;
812 // only. The case of fExitFlux=true << 656 directionTR.unit() ;
813 657
814 if(fExitFlux) << 658 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
815 { << 659 directionTR,energyTR) ;
816 const G4RotationMatrix* rotM = <<
817 pPostStepPoint->GetTouchable()->GetR <<
818 G4ThreeVector transl = pPostStepPoint- <<
819 G4AffineTransform transform = G4Affine <<
820 transform.Invert(); <<
821 G4ThreeVector localP = transform.Trans <<
822 G4ThreeVector localV = transform.Trans <<
823 660
824 G4double distance = << 661 // A XTR photon is set along the particle track and is not moved to
825 fEnvelope->GetSolid()->DistanceToOut << 662 // the G4Envelope surface as in standard X-ray TR models
826 if(verboseLevel > 1) << 663
>> 664 rand = G4UniformRand();
>> 665 G4double delta = rand*distance ;
>> 666 G4double deltaTime = delta /
>> 667 ((pPreStepPoint->GetVelocity()+
>> 668 pPostStepPoint->GetVelocity())/2.);
>> 669
>> 670 G4double aSecondaryTime = startTime + deltaTime;
>> 671
>> 672 G4ThreeVector positionTR = startPos + delta*direction ;
>> 673
>> 674 if( fExitFlux )
827 { 675 {
828 G4cout << "distance to exit = " << d << 676 const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
>> 677 G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
>> 678 G4AffineTransform transform = G4AffineTransform(rotM,transl);
>> 679 transform.Invert();
>> 680 G4ThreeVector localP = transform.TransformPoint(positionTR);
>> 681 G4ThreeVector localV = transform.TransformAxis(directionTR);
>> 682
>> 683 G4double distance = fEnvelope->GetSolid()->DistanceToOut(localP, localV);
>> 684 if(verboseLevel)
>> 685 {
>> 686 G4cout<<"distance to exit = "<<distance/mm<<" mm"<<G4endl;
>> 687 }
>> 688 positionTR += distance*directionTR;
>> 689 aSecondaryTime += distance/c_light;
829 } 690 }
830 position += distance * directionTR; << 691
831 startTime += distance / c_light; << 692 G4Track* aSecondaryTrack = new G4Track( aPhotonTR,
832 } << 693 aSecondaryTime,positionTR ) ;
833 G4Track* aSecondaryTrack = new G4Track(a << 694 aSecondaryTrack->SetTouchableHandle(aStep.GetPostStepPoint()
834 aSecondaryTrack->SetTouchableHandle( << 695 ->GetTouchableHandle());
835 aStep.GetPostStepPoint()->GetTouchable << 696 aSecondaryTrack->SetParentID(aTrack.GetTrackID());
836 aSecondaryTrack->SetParentID(aTrack.GetT <<
837 697
838 fParticleChange.AddSecondary(aSecondaryT << 698 aParticleChange.AddSecondary(aSecondaryTrack);
839 fParticleChange.ProposeEnergy(kinEnergy) << 699 }
>> 700 kinEnergy -= sumEnergyTR ;
>> 701 aParticleChange.ProposeEnergy(kinEnergy) ;
840 } 702 }
841 } 703 }
842 return G4VDiscreteProcess::PostStepDoIt(aTra << 704 // return G4VContinuousProcess::AlongStepDoIt(aTrack, aStep);
>> 705 return &aParticleChange;
843 } 706 }
844 707
>> 708
845 ////////////////////////////////////////////// 709 ///////////////////////////////////////////////////////////////////////
>> 710 //
846 // This function returns the spectral and angl 711 // This function returns the spectral and angle density of TR quanta
847 // in X-ray energy region generated forward wh 712 // in X-ray energy region generated forward when a relativistic
848 // charged particle crosses interface between 713 // charged particle crosses interface between two materials.
849 // The high energy small theta approximation i 714 // The high energy small theta approximation is applied.
850 // (matter1 -> matter2, or 2->1) 715 // (matter1 -> matter2, or 2->1)
851 // varAngle =2* (1 - std::cos(theta)) or appro << 716 // varAngle =2* (1 - cos(theta)) or approximately = theta*theta
852 G4complex G4VXTRenergyLoss::OneInterfaceXTRdEd << 717 //
853 << 718
854 { << 719 G4complex G4XTRenergyLoss::OneInterfaceXTRdEdx( G4double energy,
855 G4complex Z1 = GetPlateComplexFZ(energy, gam << 720 G4double gamma,
856 G4complex Z2 = GetGasComplexFZ(energy, gamma << 721 G4double varAngle )
>> 722 {
>> 723 G4complex Z1 = GetPlateComplexFZ(energy,gamma,varAngle) ;
>> 724 G4complex Z2 = GetGasComplexFZ(energy,gamma,varAngle) ;
>> 725
>> 726 G4complex zOut = (Z1 - Z2)*(Z1 - Z2)
>> 727 * (varAngle*energy/hbarc/hbarc) ;
>> 728 return zOut ;
857 729
858 G4complex zOut = (Z1 - Z2) * (Z1 - Z2) * (va <<
859 return zOut; <<
860 } 730 }
861 731
>> 732
862 ////////////////////////////////////////////// 733 //////////////////////////////////////////////////////////////////////////////
>> 734 //
863 // For photon energy distribution tables. Inte 735 // For photon energy distribution tables. Integrate first over angle
864 G4double G4VXTRenergyLoss::SpectralAngleXTRdEd << 736 //
>> 737
>> 738 G4double G4XTRenergyLoss::SpectralAngleXTRdEdx(G4double varAngle)
865 { 739 {
866 G4double result = GetStackFactor(fEnergy, fG << 740 G4double result = GetStackFactor(fEnergy,fGamma,varAngle) ;
867 if(result < 0.0) << 741 if(result < 0.0) result = 0.0 ;
868 result = 0.0; << 742 return result ;
869 return result; <<
870 } 743 }
871 744
872 ////////////////////////////////////////////// 745 /////////////////////////////////////////////////////////////////////////
>> 746 //
873 // For second integration over energy 747 // For second integration over energy
874 G4double G4VXTRenergyLoss::SpectralXTRdEdx(G4d << 748
>> 749 G4double G4XTRenergyLoss::SpectralXTRdEdx(G4double energy)
875 { 750 {
876 G4int i; << 751 fEnergy = energy ;
877 static constexpr G4int iMax = 8; << 752 G4Integrator<G4XTRenergyLoss,G4double(G4XTRenergyLoss::*)(G4double)> integral ;
878 G4double angleSum = 0.0; << 753 return integral.Legendre96(this,&G4XTRenergyLoss::SpectralAngleXTRdEdx,
879 << 754 0.0,0.1*fMaxThetaTR) +
880 G4double lim[iMax] = { 0.0, 0.01, 0.02, 0.05 << 755 integral.Legendre96(this,&G4XTRenergyLoss::SpectralAngleXTRdEdx,
881 << 756 0.1*fMaxThetaTR,0.2*fMaxThetaTR) +
882 for(i = 0; i < iMax; ++i) << 757 integral.Legendre96(this,&G4XTRenergyLoss::SpectralAngleXTRdEdx,
883 lim[i] *= fMaxThetaTR; << 758 0.2*fMaxThetaTR,0.4*fMaxThetaTR) +
884 << 759 integral.Legendre96(this,&G4XTRenergyLoss::SpectralAngleXTRdEdx,
885 G4Integrator<G4VXTRenergyLoss, G4double (G4V << 760 0.4*fMaxThetaTR,0.7*fMaxThetaTR) +
886 integral; << 761 integral.Legendre96(this,&G4XTRenergyLoss::SpectralAngleXTRdEdx,
887 << 762 0.7*fMaxThetaTR,fMaxThetaTR) ;
888 fEnergy = energy; << 763 }
889 { << 764
890 for(i = 0; i < iMax - 1; ++i) <<
891 { <<
892 angleSum += integral.Legendre96( <<
893 this, &G4VXTRenergyLoss::SpectralAngle <<
894 } <<
895 } <<
896 return angleSum; <<
897 } <<
898 <<
899 ////////////////////////////////////////////// 765 //////////////////////////////////////////////////////////////////////////
>> 766 //
900 // for photon angle distribution tables 767 // for photon angle distribution tables
901 G4double G4VXTRenergyLoss::AngleSpectralXTRdEd << 768 //
902 { <<
903 G4double result = GetStackFactor(energy, fGa <<
904 if(result < 0) <<
905 result = 0.0; <<
906 return result; <<
907 } <<
908 769
909 ////////////////////////////////////////////// << 770 G4double G4XTRenergyLoss::AngleSpectralXTRdEdx(G4double energy)
910 // The XTR angular distribution based on trans <<
911 G4double G4VXTRenergyLoss::AngleXTRdEdx(G4doub <<
912 { 771 {
913 G4double result; << 772 G4double result = GetStackFactor(energy,fGamma,fVarAngle) ;
914 G4double sum = 0., tmp1, tmp2, tmp = 0., cof << 773 if(result < 0) result = 0.0 ;
915 energy2; << 774 return result ;
916 G4int k, kMax, kMin, i; << 775 }
917 776
918 cofPHC = twopi * hbarc; << 777 ///////////////////////////////////////////////////////////////////////////
919 << 778 //
920 cof1 = (fPlateThick + fGasThick) * (1. / fGa << 779 //
921 cof2 = fPlateThick * fSigma1 + fGasThick * f <<
922 <<
923 cofMin = std::sqrt(cof1 * cof2); <<
924 cofMin /= cofPHC; <<
925 <<
926 kMin = G4int(cofMin); <<
927 if(cofMin > kMin) <<
928 kMin++; <<
929 <<
930 kMax = kMin + 9; <<
931 <<
932 for(k = kMin; k <= kMax; ++k) <<
933 { <<
934 tmp1 = cofPHC * k; <<
935 tmp2 = std::sqrt(tmp1 * tmp1 - cof1 * c <<
936 energy1 = (tmp1 + tmp2) / cof1; <<
937 energy2 = (tmp1 - tmp2) / cof1; <<
938 <<
939 for(i = 0; i < 2; ++i) <<
940 { <<
941 if(i == 0) <<
942 { <<
943 if(energy1 > fTheMaxEnergyTR || energy <<
944 continue; <<
945 <<
946 tmp1 = <<
947 (energy1 * energy1 * (1. / fGamma / <<
948 fPlateThick / (4 * hbarc * energy1); <<
949 tmp2 = std::sin(tmp1); <<
950 tmp = energy1 * tmp2 * tmp2; <<
951 tmp2 = fPlateThick / (4. * tmp1); <<
952 tmp1 = <<
953 hbarc * energy1 / <<
954 (energy1 * energy1 * (1. / fGamma / <<
955 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2); <<
956 tmp1 = cof1 / (4. * hbarc) - cof2 / (4 <<
957 tmp2 = std::abs(tmp1); <<
958 <<
959 if(tmp2 > 0.) <<
960 tmp /= tmp2; <<
961 else <<
962 continue; <<
963 } <<
964 else <<
965 { <<
966 if(energy2 > fTheMaxEnergyTR || energy <<
967 continue; <<
968 <<
969 tmp1 = <<
970 (energy2 * energy2 * (1. / fGamma / <<
971 fPlateThick / (4. * hbarc * energy2) <<
972 tmp2 = std::sin(tmp1); <<
973 tmp = energy2 * tmp2 * tmp2; <<
974 tmp2 = fPlateThick / (4. * tmp1); <<
975 tmp1 = <<
976 hbarc * energy2 / <<
977 (energy2 * energy2 * (1. / fGamma / <<
978 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2); <<
979 tmp1 = cof1 / (4. * hbarc) - cof2 / (4 <<
980 tmp2 = std::abs(tmp1); <<
981 <<
982 if(tmp2 > 0.) <<
983 tmp /= tmp2; <<
984 else <<
985 continue; <<
986 } <<
987 sum += tmp; <<
988 } <<
989 } <<
990 result = 4. * pi * fPlateNumber * sum * varA <<
991 result /= hbarc * hbarc; <<
992 780
993 return result; << 781 G4double G4XTRenergyLoss::AngleXTRdEdx(G4double varAngle)
>> 782 {
>> 783 fVarAngle = varAngle ;
>> 784 G4Integrator<G4XTRenergyLoss,G4double(G4XTRenergyLoss::*)(G4double)> integral ;
>> 785 return integral.Legendre10(this,&G4XTRenergyLoss::AngleSpectralXTRdEdx,
>> 786 fMinEnergyTR,fMaxEnergyTR) ;
994 } 787 }
995 788
>> 789
>> 790 //////////////////////////////////////////////////////////////////////
996 ////////////////////////////////////////////// 791 //////////////////////////////////////////////////////////////////////
>> 792 //////////////////////////////////////////////////////////////////////
>> 793 //
997 // Calculates formation zone for plates. Omega 794 // Calculates formation zone for plates. Omega is energy !!!
998 G4double G4VXTRenergyLoss::GetPlateFormationZo << 795
999 << 796 G4double G4XTRenergyLoss::GetPlateFormationZone( G4double omega ,
1000 { << 797 G4double gamma ,
1001 G4double cof, lambda; << 798 G4double varAngle )
1002 lambda = 1.0 / gamma / gamma + varAngle + f << 799 {
1003 cof = 2.0 * hbarc / omega / lambda; << 800 G4double cof, lambda ;
1004 return cof; << 801 lambda = 1.0/gamma/gamma + varAngle + fSigma1/omega/omega ;
>> 802 cof = 2.0*hbarc/omega/lambda ;
>> 803 return cof ;
1005 } 804 }
1006 805
1007 ///////////////////////////////////////////// 806 //////////////////////////////////////////////////////////////////////
>> 807 //
1008 // Calculates complex formation zone for plat 808 // Calculates complex formation zone for plates. Omega is energy !!!
1009 G4complex G4VXTRenergyLoss::GetPlateComplexFZ << 809
1010 << 810 G4complex G4XTRenergyLoss::GetPlateComplexFZ( G4double omega ,
>> 811 G4double gamma ,
>> 812 G4double varAngle )
1011 { 813 {
1012 G4double cof, length, delta, real_v, image_ << 814 G4double cof, length,delta, real, image ;
1013 815
1014 length = 0.5 * GetPlateFormationZone(omega, << 816 length = 0.5*GetPlateFormationZone(omega,gamma,varAngle) ;
1015 delta = length * GetPlateLinearPhotoAbs(om << 817 delta = length*GetPlateLinearPhotoAbs(omega) ;
1016 cof = 1.0 / (1.0 + delta * delta); << 818 cof = 1.0/(1.0 + delta*delta) ;
1017 819
1018 real_v = length * cof; << 820 real = length*cof ;
1019 image_v = real_v * delta; << 821 image = real*delta ;
1020 822
1021 G4complex zone(real_v, image_v); << 823 G4complex zone(real,image);
1022 return zone; << 824 return zone ;
1023 } 825 }
1024 826
1025 ///////////////////////////////////////////// 827 ////////////////////////////////////////////////////////////////////////
>> 828 //
1026 // Computes matrix of Sandia photo absorption 829 // Computes matrix of Sandia photo absorption cross section coefficients for
1027 // plate material 830 // plate material
1028 void G4VXTRenergyLoss::ComputePlatePhotoAbsCo <<
1029 { <<
1030 const G4MaterialTable* theMaterialTable = G <<
1031 const G4Material* mat = ( <<
1032 fPlatePhotoAbsCof = m <<
1033 831
1034 return; << 832 void G4XTRenergyLoss::ComputePlatePhotoAbsCof()
>> 833 {
>> 834 G4int i, j, numberOfElements ;
>> 835 static const G4MaterialTable*
>> 836 theMaterialTable = G4Material::GetMaterialTable();
>> 837
>> 838 G4SandiaTable thisMaterialSandiaTable(fMatIndex1) ;
>> 839 numberOfElements = (*theMaterialTable)[fMatIndex1]->GetNumberOfElements() ;
>> 840 G4int* thisMaterialZ = new G4int[numberOfElements] ;
>> 841
>> 842 for(i=0;i<numberOfElements;i++)
>> 843 {
>> 844 thisMaterialZ[i] = (G4int)(*theMaterialTable)[fMatIndex1]->
>> 845 GetElement(i)->GetZ() ;
>> 846 }
>> 847 fPlateIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
>> 848 (thisMaterialZ,numberOfElements) ;
>> 849
>> 850 fPlateIntervalNumber = thisMaterialSandiaTable.SandiaMixing
>> 851 ( thisMaterialZ ,
>> 852 (*theMaterialTable)[fMatIndex1]->GetFractionVector() ,
>> 853 numberOfElements,fPlateIntervalNumber) ;
>> 854
>> 855 fPlatePhotoAbsCof = new G4double*[fPlateIntervalNumber] ;
>> 856
>> 857 for(i=0;i<fPlateIntervalNumber;i++)
>> 858 {
>> 859 fPlatePhotoAbsCof[i] = new G4double[5] ;
>> 860 }
>> 861 for(i=0;i<fPlateIntervalNumber;i++)
>> 862 {
>> 863 fPlatePhotoAbsCof[i][0] = thisMaterialSandiaTable.
>> 864 GetPhotoAbsorpCof(i+1,0) ;
>> 865
>> 866 for(j=1;j<5;j++)
>> 867 {
>> 868 fPlatePhotoAbsCof[i][j] = thisMaterialSandiaTable.
>> 869 GetPhotoAbsorpCof(i+1,j)*
>> 870 (*theMaterialTable)[fMatIndex1]->GetDensity() ;
>> 871 }
>> 872 }
>> 873 delete[] thisMaterialZ ;
>> 874 return ;
1035 } 875 }
1036 876
1037 ///////////////////////////////////////////// 877 //////////////////////////////////////////////////////////////////////
1038 // Returns the value of linear photo absorpti << 878 //
>> 879 // Returns the value of linear photo absorption coefficient (in reciprocal
1039 // length) for plate for given energy of X-ra 880 // length) for plate for given energy of X-ray photon omega
1040 G4double G4VXTRenergyLoss::GetPlateLinearPhot << 881
>> 882 G4double G4XTRenergyLoss::GetPlateLinearPhotoAbs(G4double omega)
1041 { 883 {
1042 G4double omega2, omega3, omega4; << 884 G4int i ;
>> 885 G4double omega2, omega3, omega4 ;
>> 886
>> 887 omega2 = omega*omega ;
>> 888 omega3 = omega2*omega ;
>> 889 omega4 = omega2*omega2 ;
1043 890
1044 omega2 = omega * omega; << 891 for(i=0;i<fPlateIntervalNumber;i++)
1045 omega3 = omega2 * omega; << 892 {
1046 omega4 = omega2 * omega2; << 893 if( omega < fPlatePhotoAbsCof[i][0] ) break ;
1047 << 894 }
1048 const G4double* SandiaCof = fPlatePhotoAbsC << 895 if( i == 0 )
1049 G4double cross = SandiaCof[0] / << 896 {
1050 SandiaCof[2] / omega3 + Sa << 897 G4Exception("Invalid (<I1) energy in G4XTRenergyLoss::GetPlateLinearPhotoAbs");
1051 return cross; << 898 }
>> 899 else i-- ;
>> 900
>> 901 return fPlatePhotoAbsCof[i][1]/omega + fPlatePhotoAbsCof[i][2]/omega2 +
>> 902 fPlatePhotoAbsCof[i][3]/omega3 + fPlatePhotoAbsCof[i][4]/omega4 ;
1052 } 903 }
1053 904
1054 ///////////////////////////////////////////// 905 //////////////////////////////////////////////////////////////////////
>> 906 //
1055 // Calculates formation zone for gas. Omega i 907 // Calculates formation zone for gas. Omega is energy !!!
1056 G4double G4VXTRenergyLoss::GetGasFormationZon << 908
1057 << 909 G4double G4XTRenergyLoss::GetGasFormationZone( G4double omega ,
1058 { << 910 G4double gamma ,
1059 G4double cof, lambda; << 911 G4double varAngle )
1060 lambda = 1.0 / gamma / gamma + varAngle + f << 912 {
1061 cof = 2.0 * hbarc / omega / lambda; << 913 G4double cof, lambda ;
1062 return cof; << 914 lambda = 1.0/gamma/gamma + varAngle + fSigma2/omega/omega ;
>> 915 cof = 2.0*hbarc/omega/lambda ;
>> 916 return cof ;
>> 917
1063 } 918 }
1064 919
>> 920
1065 ///////////////////////////////////////////// 921 //////////////////////////////////////////////////////////////////////
>> 922 //
1066 // Calculates complex formation zone for gas 923 // Calculates complex formation zone for gas gaps. Omega is energy !!!
1067 G4complex G4VXTRenergyLoss::GetGasComplexFZ(G << 924
1068 G << 925 G4complex G4XTRenergyLoss::GetGasComplexFZ( G4double omega ,
>> 926 G4double gamma ,
>> 927 G4double varAngle )
1069 { 928 {
1070 G4double cof, length, delta, real_v, image_ << 929 G4double cof, length,delta, real, image ;
1071 930
1072 length = 0.5 * GetGasFormationZone(omega, g << 931 length = 0.5*GetGasFormationZone(omega,gamma,varAngle) ;
1073 delta = length * GetGasLinearPhotoAbs(omeg << 932 delta = length*GetGasLinearPhotoAbs(omega) ;
1074 cof = 1.0 / (1.0 + delta * delta); << 933 cof = 1.0/(1.0 + delta*delta) ;
1075 934
1076 real_v = length * cof; << 935 real = length*cof ;
1077 image_v = real_v * delta; << 936 image = real*delta ;
1078 937
1079 G4complex zone(real_v, image_v); << 938 G4complex zone(real,image);
1080 return zone; << 939 return zone ;
1081 } 940 }
1082 941
>> 942
>> 943
1083 ///////////////////////////////////////////// 944 ////////////////////////////////////////////////////////////////////////
>> 945 //
1084 // Computes matrix of Sandia photo absorption 946 // Computes matrix of Sandia photo absorption cross section coefficients for
1085 // gas material 947 // gas material
1086 void G4VXTRenergyLoss::ComputeGasPhotoAbsCof( << 948
>> 949 void G4XTRenergyLoss::ComputeGasPhotoAbsCof()
1087 { 950 {
1088 const G4MaterialTable* theMaterialTable = G << 951 G4int i, j, numberOfElements ;
1089 const G4Material* mat = ( << 952 static const G4MaterialTable*
1090 fGasPhotoAbsCof = m << 953 theMaterialTable = G4Material::GetMaterialTable();
1091 return; << 954
>> 955 G4SandiaTable thisMaterialSandiaTable(fMatIndex2) ;
>> 956 numberOfElements = (*theMaterialTable)[fMatIndex2]->GetNumberOfElements() ;
>> 957 G4int* thisMaterialZ = new G4int[numberOfElements] ;
>> 958
>> 959 for(i=0;i<numberOfElements;i++)
>> 960 {
>> 961 thisMaterialZ[i] = (G4int)(*theMaterialTable)[fMatIndex2]->
>> 962 GetElement(i)->GetZ() ;
>> 963 }
>> 964 fGasIntervalNumber = thisMaterialSandiaTable.SandiaIntervals
>> 965 (thisMaterialZ,numberOfElements) ;
>> 966
>> 967 fGasIntervalNumber = thisMaterialSandiaTable.SandiaMixing
>> 968 ( thisMaterialZ ,
>> 969 (*theMaterialTable)[fMatIndex2]->GetFractionVector() ,
>> 970 numberOfElements,fGasIntervalNumber) ;
>> 971
>> 972 fGasPhotoAbsCof = new G4double*[fGasIntervalNumber] ;
>> 973
>> 974 for(i=0;i<fGasIntervalNumber;i++)
>> 975 {
>> 976 fGasPhotoAbsCof[i] = new G4double[5] ;
>> 977 }
>> 978 for(i=0;i<fGasIntervalNumber;i++)
>> 979 {
>> 980 fGasPhotoAbsCof[i][0] = thisMaterialSandiaTable.
>> 981 GetPhotoAbsorpCof(i+1,0) ;
>> 982
>> 983 for(j=1;j<5;j++)
>> 984 {
>> 985 fGasPhotoAbsCof[i][j] = thisMaterialSandiaTable.
>> 986 GetPhotoAbsorpCof(i+1,j)*
>> 987 (*theMaterialTable)[fMatIndex2]->GetDensity() ;
>> 988 }
>> 989 }
>> 990 delete[] thisMaterialZ ;
>> 991 return ;
1092 } 992 }
1093 993
1094 ///////////////////////////////////////////// 994 //////////////////////////////////////////////////////////////////////
1095 // Returns the value of linear photo absorpti << 995 //
>> 996 // Returns the value of linear photo absorption coefficient (in reciprocal
1096 // length) for gas 997 // length) for gas
1097 G4double G4VXTRenergyLoss::GetGasLinearPhotoA << 998
>> 999 G4double G4XTRenergyLoss::GetGasLinearPhotoAbs(G4double omega)
1098 { 1000 {
1099 G4double omega2, omega3, omega4; << 1001 G4int i ;
>> 1002 G4double omega2, omega3, omega4 ;
>> 1003
>> 1004 omega2 = omega*omega ;
>> 1005 omega3 = omega2*omega ;
>> 1006 omega4 = omega2*omega2 ;
>> 1007
>> 1008 for(i=0;i<fGasIntervalNumber;i++)
>> 1009 {
>> 1010 if( omega < fGasPhotoAbsCof[i][0] ) break ;
>> 1011 }
>> 1012 if( i == 0 )
>> 1013 {
>> 1014 G4Exception("Invalid (<I1) energy in G4XTRenergyLoss::GetGasLinearPhotoAbs");
>> 1015 }
>> 1016 else i-- ;
>> 1017
>> 1018 return fGasPhotoAbsCof[i][1]/omega + fGasPhotoAbsCof[i][2]/omega2 +
>> 1019 fGasPhotoAbsCof[i][3]/omega3 + fGasPhotoAbsCof[i][4]/omega4 ;
1100 1020
1101 omega2 = omega * omega; <<
1102 omega3 = omega2 * omega; <<
1103 omega4 = omega2 * omega2; <<
1104 <<
1105 const G4double* SandiaCof = fGasPhotoAbsCof <<
1106 G4double cross = SandiaCof[0] / <<
1107 SandiaCof[2] / omega3 + Sa <<
1108 return cross; <<
1109 } 1021 }
1110 1022
1111 ///////////////////////////////////////////// 1023 //////////////////////////////////////////////////////////////////////
1112 // Calculates the product of linear cof by fo << 1024 //
>> 1025 // Calculates the product of linear cof by formation zone for plate.
1113 // Omega is energy !!! 1026 // Omega is energy !!!
1114 G4double G4VXTRenergyLoss::GetPlateZmuProduct << 1027
1115 << 1028 G4double G4XTRenergyLoss::GetPlateZmuProduct( G4double omega ,
>> 1029 G4double gamma ,
>> 1030 G4double varAngle )
1116 { 1031 {
1117 return GetPlateFormationZone(omega, gamma, << 1032 return GetPlateFormationZone(omega,gamma,varAngle)*GetPlateLinearPhotoAbs(omega) ;
1118 GetPlateLinearPhotoAbs(omega); <<
1119 } 1033 }
1120 ///////////////////////////////////////////// 1034 //////////////////////////////////////////////////////////////////////
1121 // Calculates the product of linear cof by fo << 1035 //
>> 1036 // Calculates the product of linear cof by formation zone for plate.
1122 // G4cout and output in file in some energy r 1037 // G4cout and output in file in some energy range.
1123 void G4VXTRenergyLoss::GetPlateZmuProduct() << 1038
>> 1039 void G4XTRenergyLoss::GetPlateZmuProduct()
1124 { 1040 {
1125 std::ofstream outPlate("plateZmu.dat", std: << 1041 ofstream outPlate("plateZmu.dat", ios::out ) ;
1126 outPlate.setf(std::ios::scientific, std::io << 1042 outPlate.setf( ios::scientific, ios::floatfield );
1127 1043
1128 G4int i; << 1044 G4int i ;
1129 G4double omega, varAngle, gamma; << 1045 G4double omega, varAngle, gamma ;
1130 gamma = 10000.; << 1046 gamma = 10000. ;
1131 varAngle = 1 / gamma / gamma; << 1047 varAngle = 1/gamma/gamma ;
1132 if(verboseLevel > 0) << 1048 G4cout<<"energy, keV"<<"\t"<<"Zmu for plate"<<G4endl ;
1133 G4cout << "energy, keV" << "\t" << "Zmu f << 1049 for(i=0;i<100;i++)
1134 for(i = 0; i < 100; ++i) << 1050 {
1135 { << 1051 omega = (1.0 + i)*keV ;
1136 omega = (1.0 + i) * keV; << 1052 G4cout<<omega/keV<<"\t"<<GetPlateZmuProduct(omega,gamma,varAngle)<<"\t" ;
1137 if(verboseLevel > 1) << 1053 outPlate<<omega/keV<<"\t\t"<<GetPlateZmuProduct(omega,gamma,varAngle)<<G4endl ;
1138 G4cout << omega / keV << "\t" <<
1139 << GetPlateZmuProduct(omega, gam <<
1140 if(verboseLevel > 0) <<
1141 outPlate << omega / keV << "\t\t" <<
1142 << GetPlateZmuProduct(omega, g <<
1143 } 1054 }
1144 return; << 1055 return ;
1145 } 1056 }
1146 1057
1147 ///////////////////////////////////////////// 1058 //////////////////////////////////////////////////////////////////////
1148 // Calculates the product of linear cof by fo << 1059 //
>> 1060 // Calculates the product of linear cof by formation zone for gas.
1149 // Omega is energy !!! 1061 // Omega is energy !!!
1150 G4double G4VXTRenergyLoss::GetGasZmuProduct(G << 1062
1151 G << 1063 G4double G4XTRenergyLoss::GetGasZmuProduct( G4double omega ,
>> 1064 G4double gamma ,
>> 1065 G4double varAngle )
1152 { 1066 {
1153 return GetGasFormationZone(omega, gamma, va << 1067 return GetGasFormationZone(omega,gamma,varAngle)*GetGasLinearPhotoAbs(omega) ;
1154 GetGasLinearPhotoAbs(omega); <<
1155 } 1068 }
1156 <<
1157 ///////////////////////////////////////////// 1069 //////////////////////////////////////////////////////////////////////
1158 // Calculates the product of linear cof by fo << 1070 //
>> 1071 // Calculates the product of linear cof byformation zone for gas.
1159 // G4cout and output in file in some energy r 1072 // G4cout and output in file in some energy range.
1160 void G4VXTRenergyLoss::GetGasZmuProduct() <<
1161 { <<
1162 std::ofstream outGas("gasZmu.dat", std::ios <<
1163 outGas.setf(std::ios::scientific, std::ios: <<
1164 G4int i; <<
1165 G4double omega, varAngle, gamma; <<
1166 gamma = 10000.; <<
1167 varAngle = 1 / gamma / gamma; <<
1168 if(verboseLevel > 0) <<
1169 G4cout << "energy, keV" << "\t" << "Zmu f <<
1170 for(i = 0; i < 100; ++i) <<
1171 { <<
1172 omega = (1.0 + i) * keV; <<
1173 if(verboseLevel > 1) <<
1174 G4cout << omega / keV << "\t" << GetGas <<
1175 << "\t"; <<
1176 if(verboseLevel > 0) <<
1177 outGas << omega / keV << "\t\t" <<
1178 << GetGasZmuProduct(omega, gamma <<
1179 } <<
1180 return; <<
1181 } <<
1182 1073
1183 ///////////////////////////////////////////// << 1074 void G4XTRenergyLoss::GetGasZmuProduct()
1184 // Computes Compton cross section for plate m <<
1185 G4double G4VXTRenergyLoss::GetPlateCompton(G4 <<
1186 { 1075 {
1187 G4int i, numberOfElements; << 1076 ofstream outGas("gasZmu.dat", ios::out ) ;
1188 G4double xSection = 0., nowZ, sumZ = 0.; << 1077 outGas.setf( ios::scientific, ios::floatfield );
1189 << 1078 G4int i ;
1190 const G4MaterialTable* theMaterialTable = G << 1079 G4double omega, varAngle, gamma ;
1191 numberOfElements = (G4int)(*theMaterialTabl << 1080 gamma = 10000. ;
1192 << 1081 varAngle = 1/gamma/gamma ;
1193 for(i = 0; i < numberOfElements; ++i) << 1082 G4cout<<"energy, keV"<<"\t"<<"Zmu for gas"<<G4endl ;
1194 { << 1083 for(i=0;i<100;i++)
1195 nowZ = (*theMaterialTable)[fMatIndex1]->G << 1084 {
1196 sumZ += nowZ; << 1085 omega = (1.0 + i)*keV ;
1197 xSection += GetComptonPerAtom(omega, nowZ << 1086 G4cout<<omega/keV<<"\t"<<GetGasZmuProduct(omega,gamma,varAngle)<<"\t" ;
1198 } << 1087 outGas<<omega/keV<<"\t\t"<<GetGasZmuProduct(omega,gamma,varAngle)<<G4endl ;
1199 xSection /= sumZ; <<
1200 xSection *= (*theMaterialTable)[fMatIndex1] <<
1201 return xSection; <<
1202 } <<
1203 <<
1204 ///////////////////////////////////////////// <<
1205 // Computes Compton cross section for gas mat <<
1206 G4double G4VXTRenergyLoss::GetGasCompton(G4do <<
1207 { <<
1208 G4double xSection = 0., sumZ = 0.; <<
1209 <<
1210 const G4MaterialTable* theMaterialTable = G <<
1211 G4int numberOfElements = (G4int)(*theMateri <<
1212 <<
1213 for (G4int i = 0; i < numberOfElements; ++i <<
1214 { <<
1215 G4double nowZ = (*theMaterialTable)[fMatI <<
1216 sumZ += nowZ; <<
1217 xSection += GetComptonPerAtom(omega, nowZ <<
1218 } 1088 }
1219 if (sumZ > 0.0) { xSection /= sumZ; } << 1089 return ;
1220 xSection *= (*theMaterialTable)[fMatIndex2] <<
1221 return xSection; <<
1222 } <<
1223 <<
1224 ///////////////////////////////////////////// <<
1225 // Computes Compton cross section per atom wi <<
1226 // the energy GammaEnergy <<
1227 G4double G4VXTRenergyLoss::GetComptonPerAtom( <<
1228 { <<
1229 G4double CrossSection = 0.0; <<
1230 if(Z < 0.9999) <<
1231 return CrossSection; <<
1232 if(GammaEnergy < 0.1 * keV) <<
1233 return CrossSection; <<
1234 if(GammaEnergy > (100. * GeV / Z)) <<
1235 return CrossSection; <<
1236 <<
1237 static constexpr G4double a = 20.0; <<
1238 static constexpr G4double b = 230.0; <<
1239 static constexpr G4double c = 440.0; <<
1240 <<
1241 static constexpr G4double d1 = 2.7965e-1 * <<
1242 d3 = 6.7527 * bar <<
1243 e1 = 1.9756e-5 * <<
1244 e3 = -7.3913e-2 * <<
1245 f1 = -3.9178e-7 * <<
1246 f3 = 6.0480e-5 * <<
1247 <<
1248 G4double p1Z = Z * (d1 + e1 * Z + f1 * Z * <<
1249 G4double p2Z = Z * (d2 + e2 * Z + f2 * Z * <<
1250 G4double p3Z = Z * (d3 + e3 * Z + f3 * Z * <<
1251 G4double p4Z = Z * (d4 + e4 * Z + f4 * Z * <<
1252 <<
1253 G4double T0 = 15.0 * keV; <<
1254 if(Z < 1.5) <<
1255 T0 = 40.0 * keV; <<
1256 <<
1257 G4double X = std::max(GammaEnergy, T0) / el <<
1258 CrossSection = <<
1259 p1Z * std::log(1. + 2. * X) / X + <<
1260 (p2Z + p3Z * X + p4Z * X * X) / (1. + a * <<
1261 <<
1262 // modification for low energy. (special c <<
1263 if(GammaEnergy < T0) <<
1264 { <<
1265 G4double dT0 = 1. * keV; <<
1266 X = (T0 + dT0) / electron_mass <<
1267 G4double sigma = <<
1268 p1Z * std::log(1. + 2. * X) / X + <<
1269 (p2Z + p3Z * X + p4Z * X * X) / (1. + a <<
1270 G4double c1 = -T0 * (sigma - CrossSection <<
1271 G4double c2 = 0.150; <<
1272 if(Z > 1.5) <<
1273 c2 = 0.375 - 0.0556 * std::log(Z); <<
1274 G4double y = std::log(GammaEnergy / T0); <<
1275 CrossSection *= std::exp(-y * (c1 + c2 * <<
1276 } <<
1277 return CrossSection; <<
1278 } 1090 }
1279 1091
1280 ///////////////////////////////////////////// 1092 ///////////////////////////////////////////////////////////////////////
>> 1093 //
1281 // This function returns the spectral and ang 1094 // This function returns the spectral and angle density of TR quanta
1282 // in X-ray energy region generated forward w 1095 // in X-ray energy region generated forward when a relativistic
1283 // charged particle crosses interface between 1096 // charged particle crosses interface between two materials.
1284 // The high energy small theta approximation 1097 // The high energy small theta approximation is applied.
1285 // (matter1 -> matter2, or 2->1) 1098 // (matter1 -> matter2, or 2->1)
1286 // varAngle =2* (1 - std::cos(theta)) or appr << 1099 // varAngle =2* (1 - cos(theta)) or approximately = theta*theta
1287 G4double G4VXTRenergyLoss::OneBoundaryXTRNden << 1100 //
1288 << 1101
1289 << 1102 G4double
1290 { << 1103 G4XTRenergyLoss::OneBoundaryXTRNdensity( G4double energy,G4double gamma,
1291 G4double formationLength1, formationLength2 << 1104 G4double varAngle ) const
1292 formationLength1 = << 1105 {
1293 1.0 / (1.0 / (gamma * gamma) + fSigma1 / << 1106 G4double formationLength1, formationLength2 ;
1294 formationLength2 = << 1107 formationLength1 = 1.0/
1295 1.0 / (1.0 / (gamma * gamma) + fSigma2 / << 1108 (1.0/(gamma*gamma)
1296 return (varAngle / energy) * (formationLeng << 1109 + fSigma1/(energy*energy)
1297 (formationLength1 - formationLength2 << 1110 + varAngle) ;
>> 1111 formationLength2 = 1.0/
>> 1112 (1.0/(gamma*gamma)
>> 1113 + fSigma2/(energy*energy)
>> 1114 + varAngle) ;
>> 1115 return (varAngle/energy)*(formationLength1 - formationLength2)
>> 1116 *(formationLength1 - formationLength2) ;
>> 1117
1298 } 1118 }
1299 1119
1300 G4double G4VXTRenergyLoss::GetStackFactor(G4d << 1120 G4double G4XTRenergyLoss::GetStackFactor( G4double energy, G4double gamma,
1301 G4d << 1121 G4double varAngle )
1302 { 1122 {
1303 // return stack factor corresponding to one 1123 // return stack factor corresponding to one interface
1304 return std::real(OneInterfaceXTRdEdx(energy << 1124
>> 1125 return std::real( OneInterfaceXTRdEdx(energy,gamma,varAngle) );
1305 } 1126 }
1306 1127
1307 ///////////////////////////////////////////// 1128 //////////////////////////////////////////////////////////////////////////////
>> 1129 //
1308 // For photon energy distribution tables. Int 1130 // For photon energy distribution tables. Integrate first over angle
1309 G4double G4VXTRenergyLoss::XTRNSpectralAngleD << 1131 //
>> 1132
>> 1133 G4double G4XTRenergyLoss::XTRNSpectralAngleDensity(G4double varAngle)
1310 { 1134 {
1311 return OneBoundaryXTRNdensity(fEnergy, fGam << 1135 return OneBoundaryXTRNdensity(fEnergy,fGamma,varAngle)*
1312 GetStackFactor(fEnergy, fGamma, varA << 1136 GetStackFactor(fEnergy,fGamma,varAngle) ;
1313 } 1137 }
1314 1138
1315 ///////////////////////////////////////////// 1139 /////////////////////////////////////////////////////////////////////////
>> 1140 //
1316 // For second integration over energy 1141 // For second integration over energy
1317 G4double G4VXTRenergyLoss::XTRNSpectralDensit << 1142
>> 1143 G4double G4XTRenergyLoss::XTRNSpectralDensity(G4double energy)
1318 { 1144 {
1319 fEnergy = energy; << 1145 fEnergy = energy ;
1320 G4Integrator<G4VXTRenergyLoss, G4double (G4 << 1146 G4Integrator<G4XTRenergyLoss,G4double(G4XTRenergyLoss::*)(G4double)> integral ;
1321 integral; << 1147 return integral.Legendre96(this,&G4XTRenergyLoss::XTRNSpectralAngleDensity,
1322 return integral.Legendre96(this, &G4VXTRene << 1148 0.0,0.2*fMaxThetaTR) +
1323 0.0, 0.2 * fMaxT << 1149 integral.Legendre10(this,&G4XTRenergyLoss::XTRNSpectralAngleDensity,
1324 integral.Legendre10(this, &G4VXTRene << 1150 0.2*fMaxThetaTR,fMaxThetaTR) ;
1325 0.2 * fMaxThetaT << 1151 }
1326 } << 1152
1327 <<
1328 ///////////////////////////////////////////// 1153 //////////////////////////////////////////////////////////////////////////
>> 1154 //
1329 // for photon angle distribution tables 1155 // for photon angle distribution tables
1330 G4double G4VXTRenergyLoss::XTRNAngleSpectralD << 1156 //
>> 1157
>> 1158 G4double G4XTRenergyLoss::XTRNAngleSpectralDensity(G4double energy)
1331 { 1159 {
1332 return OneBoundaryXTRNdensity(energy, fGamm << 1160 return OneBoundaryXTRNdensity(energy,fGamma,fVarAngle)*
1333 GetStackFactor(energy, fGamma, fVarA << 1161 GetStackFactor(energy,fGamma,fVarAngle) ;
1334 } << 1162 }
1335 1163
1336 ///////////////////////////////////////////// 1164 ///////////////////////////////////////////////////////////////////////////
1337 G4double G4VXTRenergyLoss::XTRNAngleDensity(G << 1165 //
>> 1166 //
>> 1167
>> 1168 G4double G4XTRenergyLoss::XTRNAngleDensity(G4double varAngle)
1338 { 1169 {
1339 fVarAngle = varAngle; << 1170 fVarAngle = varAngle ;
1340 G4Integrator<G4VXTRenergyLoss, G4double (G4 << 1171 G4Integrator<G4XTRenergyLoss,G4double(G4XTRenergyLoss::*)(G4double)> integral ;
1341 integral; << 1172 return integral.Legendre96(this,&G4XTRenergyLoss::XTRNAngleSpectralDensity,
1342 return integral.Legendre96(this, &G4VXTRene << 1173 fMinEnergyTR,fMaxEnergyTR) ;
1343 fMinEnergyTR, fM <<
1344 } 1174 }
1345 1175
1346 ///////////////////////////////////////////// 1176 //////////////////////////////////////////////////////////////////////////////
1347 // Check number of photons for a range of Lor << 1177 //
>> 1178 // Check number of photons for a range of Lorentz factors from both energy
1348 // and angular tables 1179 // and angular tables
1349 void G4VXTRenergyLoss::GetNumberOfPhotons() << 1180
>> 1181 void G4XTRenergyLoss::GetNumberOfPhotons()
1350 { 1182 {
1351 G4int iTkin; << 1183 G4int iTkin ;
1352 G4double gamma, numberE; << 1184 G4double gamma, numberE ;
1353 1185
1354 std::ofstream outEn("numberE.dat", std::ios << 1186 ofstream outEn("numberE.dat", ios::out ) ;
1355 outEn.setf(std::ios::scientific, std::ios:: << 1187 outEn.setf( ios::scientific, ios::floatfield );
1356 1188
1357 std::ofstream outAng("numberAng.dat", std:: << 1189 ofstream outAng("numberAng.dat", ios::out ) ;
1358 outAng.setf(std::ios::scientific, std::ios: << 1190 outAng.setf( ios::scientific, ios::floatfield );
1359 << 1191
1360 for(iTkin = 0; iTkin < fTotBin; ++iTkin) / << 1192 for(iTkin=0;iTkin<fTotBin;iTkin++) // Lorentz factor loop
1361 { << 1193 {
1362 gamma = << 1194 gamma = 1.0 + (fProtonEnergyVector->
1363 1.0 + (fProtonEnergyVector->GetLowEdgeE << 1195 GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
1364 numberE = (*(*fEnergyDistrTable)(iTkin))( << 1196 numberE = (*(*fEnergyDistrTable)(iTkin))(0) ;
1365 if(verboseLevel > 1) << 1197 // numberA = (*(*fAngleDistrTable)(iTkin))(0) ;
1366 G4cout << gamma << "\t\t" << numberE << << 1198 G4cout<<gamma<<"\t\t"<<numberE<<"\t" // <<numberA
1367 if(verboseLevel > 0) << 1199 <<G4endl ;
1368 outEn << gamma << "\t\t" << numberE << << 1200 outEn<<gamma<<"\t\t"<<numberE<<G4endl ;
>> 1201 // outAng<<gamma<<"\t\t"<<numberA<<G4endl ;
1369 } 1202 }
1370 return; << 1203 return ;
1371 } << 1204 }
1372 1205
1373 ///////////////////////////////////////////// 1206 /////////////////////////////////////////////////////////////////////////
1374 // Returns random energy of a X-ray TR photon << 1207 //
>> 1208 // Returns randon energy of a X-ray TR photon for given scaled kinetic energy
1375 // of a charged particle 1209 // of a charged particle
1376 G4double G4VXTRenergyLoss::GetXTRrandomEnergy << 1210
>> 1211 G4double G4XTRenergyLoss::GetXTRrandomEnergy( G4double scaledTkin, G4int iTkin )
1377 { 1212 {
1378 G4int iTransfer, iPlace; << 1213 G4int iTransfer, iPlace ;
1379 G4double transfer = 0.0, position, E1, E2, << 1214 G4double transfer = 0.0, position, E1, E2, W1, W2, W ;
>> 1215
>> 1216 iPlace = iTkin - 1 ;
1380 1217
1381 iPlace = iTkin - 1; << 1218 // G4cout<<"iPlace = "<<iPlace<<endl ;
1382 1219
1383 if(iTkin == fTotBin) // relativistic plato << 1220 if(iTkin == fTotBin) // relativistic plato, try from left
1384 { 1221 {
1385 position = (*(*fEnergyDistrTable)(iPlace) << 1222 position = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand() ;
1386 1223
1387 for(iTransfer = 0;; ++iTransfer) << 1224 for(iTransfer=0;;iTransfer++)
1388 { << 1225 {
1389 if(position >= (*(*fEnergyDistrTable)(i << 1226 if(position >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break ;
1390 break; << 1227 }
1391 } << 1228 transfer = GetXTRenergy(iPlace,position,iTransfer);
1392 transfer = GetXTRenergy(iPlace, position, <<
1393 } 1229 }
1394 else 1230 else
1395 { 1231 {
1396 E1 = fProtonEnergyVector->GetLowEdgeEnerg << 1232 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
1397 E2 = fProtonEnergyVector->GetLowEdgeEnerg << 1233 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
1398 W = 1.0 / (E2 - E1); << 1234 W = 1.0/(E2 - E1) ;
1399 W1 = (E2 - scaledTkin) * W; << 1235 W1 = (E2 - scaledTkin)*W ;
1400 W2 = (scaledTkin - E1) * W; << 1236 W2 = (scaledTkin - E1)*W ;
1401 << 1237
1402 position = ((*(*fEnergyDistrTable)(iPlace << 1238 position =( (*(*fEnergyDistrTable)(iPlace))(0)*W1 +
1403 (*(*fEnergyDistrTable)(iPlace << 1239 (*(*fEnergyDistrTable)(iPlace+1))(0)*W2 )*G4UniformRand() ;
1404 G4UniformRand(); << 1240
1405 << 1241 // G4cout<<position<<"\t" ;
1406 for(iTransfer = 0;; ++iTransfer) << 1242
1407 { << 1243 for(iTransfer=0;;iTransfer++)
1408 if(position >= ((*(*fEnergyDistrTable)( << 1244 {
1409 (*(*fEnergyDistrTable)( << 1245 if( position >=
1410 break; << 1246 ( (*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1 +
1411 } << 1247 (*(*fEnergyDistrTable)(iPlace+1))(iTransfer)*W2) ) break ;
1412 transfer = GetXTRenergy(iPlace, position, << 1248 }
1413 } << 1249 transfer = GetXTRenergy(iPlace,position,iTransfer);
1414 if(transfer < 0.0) << 1250
1415 transfer = 0.0; << 1251 }
1416 return transfer; << 1252 // G4cout<<"XTR transfer = "<<transfer/keV<<" keV"<<endl ;
>> 1253 if(transfer < 0.0 ) transfer = 0.0 ;
>> 1254 return transfer ;
1417 } 1255 }
1418 1256
1419 ///////////////////////////////////////////// 1257 ////////////////////////////////////////////////////////////////////////
>> 1258 //
1420 // Returns approximate position of X-ray phot 1259 // Returns approximate position of X-ray photon energy during random sampling
1421 // over integral energy distribution 1260 // over integral energy distribution
1422 G4double G4VXTRenergyLoss::GetXTRenergy(G4int << 1261
>> 1262 G4double G4XTRenergyLoss::GetXTRenergy( G4int iPlace,
>> 1263 G4double position,
>> 1264 G4int iTransfer )
1423 { 1265 {
1424 G4double x1, x2, y1, y2, result; << 1266 G4double x1, x2, y1, y2, result ;
1425 1267
1426 if(iTransfer == 0) 1268 if(iTransfer == 0)
1427 { 1269 {
1428 result = (*fEnergyDistrTable)(iPlace)->Ge << 1270 result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1429 } << 1271 }
1430 else 1272 else
1431 { 1273 {
1432 y1 = (*(*fEnergyDistrTable)(iPlace))(iTra << 1274 y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer-1) ;
1433 y2 = (*(*fEnergyDistrTable)(iPlace))(iTra << 1275 y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer) ;
1434 1276
1435 x1 = (*fEnergyDistrTable)(iPlace)->GetLow << 1277 x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
1436 x2 = (*fEnergyDistrTable)(iPlace)->GetLow << 1278 x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1437 1279
1438 if(x1 == x2) << 1280 if ( x1 == x2 ) result = x2 ;
1439 result = x2; <<
1440 else 1281 else
1441 { 1282 {
1442 if(y1 == y2) << 1283 if ( y1 == y2 ) result = x1 + (x2 - x1)*G4UniformRand() ;
1443 result = x1 + (x2 - x1) * G4UniformRa <<
1444 else 1284 else
1445 { 1285 {
1446 result = x1 + (x2 - x1) * G4UniformRa << 1286 result = x1 + (position - y1)*(x2 - x1)/(y2 - y1) ;
1447 } 1287 }
1448 } 1288 }
1449 } 1289 }
1450 return result; << 1290 return result ;
1451 } <<
1452 <<
1453 ///////////////////////////////////////////// <<
1454 // Get XTR photon angle at given energy and <<
1455 <<
1456 G4double G4VXTRenergyLoss::GetRandomAngle(G4d <<
1457 { <<
1458 G4int iTR, iAngle; <<
1459 G4double position, angle; <<
1460 <<
1461 if(iTkin == fTotBin) <<
1462 --iTkin; <<
1463 <<
1464 fAngleForEnergyTable = fAngleBank[iTkin]; <<
1465 <<
1466 for(iTR = 0; iTR < fBinTR; ++iTR) <<
1467 { <<
1468 if(energyXTR < fXTREnergyVector->GetLowEd <<
1469 break; <<
1470 } <<
1471 if(iTR == fBinTR) <<
1472 --iTR; <<
1473 <<
1474 position = (*(*fAngleForEnergyTable)(iTR))( <<
1475 // position = (*(*fAngleForEnergyTable)(iTR <<
1476 <<
1477 for(iAngle = 0;; ++iAngle) <<
1478 // for(iAngle = 1;; ++iAngle) // ATLAS TB <<
1479 { <<
1480 if(position >= (*(*fAngleForEnergyTable)( <<
1481 break; <<
1482 } <<
1483 angle = GetAngleXTR(iTR, position, iAngle); <<
1484 return angle; <<
1485 } 1291 }
1486 1292
1487 ///////////////////////////////////////////// <<
1488 // Returns approximate position of X-ray phot <<
1489 // random sampling over integral energy distr <<
1490 1293
1491 G4double G4VXTRenergyLoss::GetAngleXTR(G4int <<
1492 G4int <<
1493 { <<
1494 G4double x1, x2, y1, y2, result; <<
1495 1294
1496 if( iTransfer == 0 ) << 1295 //
1497 // if( iTransfer == 1 ) // ATLAS TB << 1296 //
1498 { << 1297 ///////////////////////////////////////////////////////////////////////
1499 result = (*fAngleForEnergyTable)(iPlace)- <<
1500 } <<
1501 else <<
1502 { <<
1503 y1 = (*(*fAngleForEnergyTable)(iPlace))(i <<
1504 y2 = (*(*fAngleForEnergyTable)(iPlace))(i <<
1505 <<
1506 x1 = (*fAngleForEnergyTable)(iPlace)->Get <<
1507 x2 = (*fAngleForEnergyTable)(iPlace)->Get <<
1508 1298
1509 if(x1 == x2) result = x2; <<
1510 else <<
1511 { <<
1512 if( y1 == y2 ) result = x1 + (x2 - x1) <<
1513 else <<
1514 { <<
1515 result = x1 + (position - y1) * (x2 - <<
1516 // result = x1 + 0.1*(position - y1) <<
1517 // result = x1 + 0.05*(position - y1) <<
1518 } <<
1519 } <<
1520 } <<
1521 return result; <<
1522 } <<
1523 1299