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