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