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