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