Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/xrays/src/G4VXTRenergyLoss.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // 2001-2002 R&D by V.Grichine
 28 // 19.06.03 V. Grichine, modifications in BuildTable for the integration
 29 //                       in respect of angle: range is increased, accuracy is
 30 //                       improved
 31 // 28.07.05, P.Gumplinger add G4ProcessType to constructor
 32 // 28.09.07, V.Ivanchenko general cleanup without change of algorithms
 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(G4LogicalVolume* anEnvelope,
 60                                    G4Material* foilMat, G4Material* gasMat,
 61                                    G4double a, G4double b, G4int n,
 62                                    const G4String& processName,
 63                                    G4ProcessType type)
 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("model_XTRenergyLoss");
 75   SetProcessSubType(fTransitionRadiation);
 76 
 77   fPtrGamma    = nullptr;
 78   fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = fGamma = fEnergy = 0.0;
 79   fVarAngle = fLambda = fTotalDist = fPlateThick = fGasThick = 0.0;
 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 scale 
 90   fBinTR  = 100; //   number of bins in TR vectors
 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, fMaxProtonTkin, fTotBin);
101 
102   fXTREnergyVector =
103     new G4PhysicsLogVector(fTheMinEnergyTR, fTheMaxEnergyTR, fBinTR);
104 
105   fEnvelope = anEnvelope;
106 
107   fPlateNumber = n;
108   if(verboseLevel > 0)
109     G4cout << "### G4VXTRenergyLoss: the number of TR radiator plates = "
110            << fPlateNumber << G4endl;
111   if(fPlateNumber == 0)
112   {
113     G4Exception("G4VXTRenergyLoss::G4VXTRenergyLoss()", "VXTRELoss01",
114                 FatalException, "No plates in X-ray TR radiator");
115   }
116   // default is XTR dEdx, not flux after radiator
117   fExitFlux      = false;
118   // default angle distribution according numerical integration
119   fFastAngle     = false; // no angle according sum of delta-functions by default
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 + fGasThick);
129   if(verboseLevel > 0)
130     G4cout << "total radiator thickness = " << fTotalDist / cm << " cm"
131            << G4endl;
132 
133   // index of plate material
134   fMatIndex1 = (G4int)foilMat->GetIndex();
135   if(verboseLevel > 0)
136     G4cout << "plate material = " << foilMat->GetName() << G4endl;
137 
138   // index of gas material
139   fMatIndex2 = (G4int)gasMat->GetIndex();
140   if(verboseLevel > 0)
141     G4cout << "gas material = " << gasMat->GetName() << G4endl;
142 
143   // plasma energy squared for plate material
144   fSigma1 = fPlasmaCof * foilMat->GetElectronDensity();
145   if(verboseLevel > 0)
146     G4cout << "plate plasma energy = " << std::sqrt(fSigma1) / eV << " eV"
147            << G4endl;
148 
149   // plasma energy squared for gas material
150   fSigma2 = fPlasmaCof * gasMat->GetElectronDensity();
151   if(verboseLevel > 0)
152     G4cout << "gas plasma energy = " << std::sqrt(fSigma2) / eV << " eV"
153            << G4endl;
154 
155   // Compute cofs for preparation of linear photo absorption
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::ostream& out) const
185 {
186   out << "Base class for 'fast' parameterisation model describing X-ray "
187          "transition\n"
188          "radiation. Angular distribution is very rough.\n";
189 }
190 
191 ///////////////////////////////////////////////////////////////////////////////
192 // Returns condition for application of the model depending on particle type
193 G4bool G4VXTRenergyLoss::IsApplicable(const G4ParticleDefinition& particle)
194 {
195   return (particle.GetPDGCharge() != 0.0);
196 }
197 
198 /////////////////////////////////////////////////////////////////////////////////
199 // Calculate step size for XTR process inside raaditor
200 G4double G4VXTRenergyLoss::GetMeanFreePath(const G4Track& aTrack, G4double,
201                                            G4ForceCondition* condition)
202 {
203   G4int iTkin, iPlace;
204   G4double lambda, sigma, kinEnergy, mass, gamma;
205   G4double charge, chargeSq, massRatio, TkinScaled;
206   G4double E1, E2, W, W1, W2;
207 
208   *condition = NotForced;
209 
210   if(aTrack.GetVolume()->GetLogicalVolume() != fEnvelope)
211     lambda = DBL_MAX;
212   else
213   {
214     const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
215     kinEnergy                          = aParticle->GetKineticEnergy();
216     mass  = aParticle->GetDefinition()->GetPDGMass();
217     gamma = 1.0 + kinEnergy / mass;
218     if(verboseLevel > 1)
219     {
220       G4cout << " gamma = " << gamma << ";   fGamma = " << fGamma << G4endl;
221     }
222 
223     if(std::fabs(gamma - fGamma) < 0.05 * gamma)
224       lambda = fLambda;
225     else
226     {
227       charge     = aParticle->GetDefinition()->GetPDGCharge();
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->GetLowEdgeEnergy(iTkin))
235           break;
236       }
237       iPlace = iTkin - 1;
238 
239       if(iTkin == 0)
240         lambda = DBL_MAX;  // Tkin is too small, neglect of TR photon generation
241       else  // general case: Tkin between two vectors of the material
242       {
243         if(iTkin == fTotBin)
244         {
245           sigma = (*(*fEnergyDistrTable)(iPlace))(0) * chargeSq;
246         }
247         else
248         {
249           E1    = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
250           E2    = fProtonEnergyVector->GetLowEdgeEnergy(iTkin);
251           W     = 1.0 / (E2 - E1);
252           W1    = (E2 - TkinScaled) * W;
253           W2    = (TkinScaled - E1) * W;
254           sigma = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 +
255                    (*(*fEnergyDistrTable)(iPlace + 1))(0) * W2) *
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 / mm << " mm" << G4endl;
267         }
268       }
269     }
270   }
271   return lambda;
272 }
273 
274 //////////////////////////////////////////////////////////////////////////
275 // Interface for build table from physics list
276 void G4VXTRenergyLoss::BuildPhysicsTable(const G4ParticleDefinition& pd)
277 {
278   if(pd.GetPDGCharge() == 0.)
279   {
280     G4Exception("G4VXTRenergyLoss::BuildPhysicsTable", "Notification",
281                 JustWarning, "XTR initialisation for neutral particle ?!");
282   }
283   BuildEnergyTable();
284 
285   if(fAngleRadDistr)
286   {
287     if(verboseLevel > 0)
288     {
289       G4cout
290         << "Build angle for energy distribution according the current radiator"
291         << G4endl;
292     }
293     BuildAngleForEnergyBank();
294   }
295 }
296 
297 //////////////////////////////////////////////////////////////////////////
298 // Build integral energy distribution of XTR photons
299 void G4VXTRenergyLoss::BuildEnergyTable()
300 {
301   G4int iTkin, iTR, iPlace;
302   G4double radiatorCof = 1.0;  // for tuning of XTR yield
303   G4double energySum   = 0.0;
304 
305   fEnergyDistrTable = new G4PhysicsTable(fTotBin);
306   if(fAngleRadDistr)
307     fAngleDistrTable = new G4PhysicsTable(fTotBin);
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 (G4VXTRenergyLoss::*)(G4double)>
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)  // Lorentz factor loop
338   {
339     auto energyVector =
340       new G4PhysicsLogVector(fMinEnergyTR, fMaxEnergyTR, fBinTR);
341 
342     fGamma =
343       1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2);
344 
345     // if(fMaxThetaTR > fTheMaxAngle)     fMaxThetaTR = fTheMaxAngle;
346     // else if(fMaxThetaTR < fTheMinAngle)     fMaxThetaTR = fTheMinAngle;
347 
348     energySum = 0.0;
349 
350     energyVector->PutValue(fBinTR - 1, energySum);
351 
352     for(iTR = fBinTR - 2; iTR >= 0; --iTR)
353     {
354       // Legendre96 or Legendre10
355 
356       energySum += radiatorCof * fCofTR *
357   
358   // integral.Legendre10(this, &G4VXTRenergyLoss::SpectralXTRdEdx,
359   
360                    integral.Legendre96(this, &G4VXTRenergyLoss::SpectralXTRdEdx,
361                
362                                        energyVector->GetLowEdgeEnergy(iTR),
363                                        energyVector->GetLowEdgeEnergy(iTR + 1));
364 
365       energyVector->PutValue(iTR, energySum / fTotalDist);
366     }
367     iPlace = iTkin;
368     fEnergyDistrTable->insertAt(iPlace, energyVector);
369 
370     if(verboseLevel > 0)
371     {
372       G4cout << fGamma << "\t" << energySum << G4endl;
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 energy loss tables = "
381            << timer.GetUserElapsed() << " s" << G4endl;
382   }
383   fGamma = 0.;
384   return;
385 }
386 
387 //////////////////////////////////////////////////////////////////////////
388 // Bank of angle distributions for given energies (slow!)
389 
390 void G4VXTRenergyLoss::BuildAngleForEnergyBank()
391 {
392   
393   if( ( this->GetProcessName() == "TranspRegXTRadiator" ||
394         this->GetProcessName() == "TranspRegXTRmodel" ||
395         this->GetProcessName() == "RegularXTRadiator" ||
396   this->GetProcessName() == "RegularXTRmodel"  )       && fFastAngle    ) // ffastAngle=true!
397   {
398     BuildAngleTable(); // by sum of delta-functions
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, fMaxEnergyTR, fBinTR);
419 
420   G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)>
421     integral;
422 
423   G4cout.precision(4);
424   G4Timer timer;
425   timer.Start();
426 
427   for(iTkin = 0; iTkin < fTotBin; ++iTkin)  // Lorentz factor loop
428   {
429     fGamma =
430       1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2);
431 
432     if(fMaxThetaTR > fTheMaxAngle)
433       fMaxThetaTR = fTheMaxAngle;
434     else if(fMaxThetaTR < fTheMinAngle)
435       fMaxThetaTR = fTheMinAngle;
436 
437     fAngleForEnergyTable = new G4PhysicsTable(fBinTR);
438 
439     for(iTR = 0; iTR < fBinTR; ++iTR)
440     {
441       angleSum = 0.0;
442       fEnergy  = energyVector->GetLowEdgeEnergy(iTR);
443       
444      // log-vector to increase number of thin bins for small angles
445       auto angleVector = new G4PhysicsLogVector(fMinThetaTR, fMaxThetaTR, fBinTR);
446  
447       
448 
449       angleVector->PutValue(fBinTR - 1, angleSum);
450 
451       for(i = fBinTR - 2; i >= 0; --i)
452       {
453         // Legendre96 or Legendre10
454 
455         angleSum +=
456           integral.Legendre10(this, &G4VXTRenergyLoss::SpectralAngleXTRdEdx,
457                               angleVector->GetLowEdgeEnergy(i),
458                               angleVector->GetLowEdgeEnergy(i + 1));
459 
460         angleVector->PutValue(i, angleSum);
461       }
462       fAngleForEnergyTable->insertAt(iTR, angleVector);
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 angle for energy loss tables = "
472            << timer.GetUserElapsed() << " s" << G4endl;
473   }
474   fGamma = 0.;
475   delete energyVector;
476 }
477 
478 ////////////////////////////////////////////////////////////////////////
479 // Build XTR angular distribution at given energy based on the model
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" << "\t"
505            << "XTR photon number" << G4endl << G4endl;
506   }
507   for(iTkin = 0; iTkin < fTotBin; ++iTkin)  // Lorentz factor loop
508   {
509     fGamma =
510       1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2);
511 
512     // fMaxThetaTR = 25. * 2500.0 / (fGamma * fGamma);  // theta^2
513 
514     if(fMaxThetaTR > fTheMaxAngle)
515       fMaxThetaTR = fTheMaxAngle;
516     else
517     {
518       if(fMaxThetaTR < fTheMinAngle)
519         fMaxThetaTR = fTheMinAngle;
520     }
521 
522     fAngleForEnergyTable = new G4PhysicsTable(fBinTR);
523 
524     for(iTR = 0; iTR < fBinTR; ++iTR)
525     {
526       energy = fXTREnergyVector->GetLowEdgeEnergy(iTR);
527 
528       G4PhysicsFreeVector* angleVector = GetAngleVector(energy, fBinTR);
529 
530       fAngleForEnergyTable->insertAt(iTR, angleVector);
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 for given energy tables = "
540            << timer.GetUserElapsed() << " s" << G4endl;
541   }
542   fGamma = 0.;
543 
544   return;
545 }
546 
547 /////////////////////////////////////////////////////////////////////////
548 // Vector of angles and angle integral distributions
549 G4PhysicsFreeVector* G4VXTRenergyLoss::GetAngleVector(G4double energy, G4int n)
550 {
551   G4double theta = 0., result, tmp = 0., cof1, cof2, cofMin, cofPHC,
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 / energy;
559   cof1   = fPlateThick * tmp;
560   cof2   = fGasThick * tmp;
561 
562   cofMin = energy * (fPlateThick + fGasThick) / fGamma / fGamma;
563   cofMin += (fPlateThick * fSigma1 + fGasThick * fSigma2) / energy;
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(fMaxThetaTR) * fGamma
574            << "; tmp = " << 0. << ";    angleSum = " << angleSum << G4endl;
575   }
576 
577   for(iTheta = n - 1; iTheta >= 1; --iTheta)
578   {
579     k      = iTheta - 1 + kMin;
580     tmp    = pi * fPlateThick * (k + cof2) / (fPlateThick + fGasThick);
581     result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2);
582     tmp    = std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result;
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 / energy / (fPlateThick + fGasThick);
596 
597     if(verboseLevel > 2)
598     {
599       G4cout << "iTheta = " << iTheta << "; k = " << k
600              << "; theta = " << std::sqrt(theta) * fGamma << "; tmp = " << tmp
601              << ";    angleSum = " << angleSum << G4endl;
602     }
603     angleVector->PutValue(iTheta, theta, angleSum);
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 << "; theta = " << std::sqrt(theta) * fGamma
614            << "; tmp = " << tmp << ";    angleSum = " << angleSum << G4endl;
615   }
616   angleVector->PutValue(iTheta, theta, angleSum);
617 
618   return angleVector;
619 }
620 
621 ////////////////////////////////////////////////////////////////////////
622 // Build XTR angular distribution based on the model of transparent regular
623 // radiator
624 void G4VXTRenergyLoss::BuildGlobalAngleTable()
625 {
626   G4int iTkin, iTR, iPlace;
627   G4double radiatorCof = 1.0;  // for tuning of XTR yield
628   G4double angleSum;
629   fAngleDistrTable = new G4PhysicsTable(fTotBin);
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)  // Lorentz factor loop
656   {
657     fGamma =
658       1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2);
659 
660     // fMaxThetaTR = 25.0 / (fGamma * fGamma);  // theta^2
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, fMaxThetaTR, fBinTR);
673     //  new G4PhysicsLogVector(1.e-8, fMaxThetaTR, fBinTR);
674 
675     angleSum = 0.0;
676 
677     G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)>
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, &G4VXTRenergyLoss::AngleXTRdEdx,
686                                       angleVector->GetLowEdgeEnergy(iTR),
687                                       angleVector->GetLowEdgeEnergy(iTR + 1));
688 
689       angleVector->PutValue(iTR, angleSum);
690     }
691     if(verboseLevel > 1)
692     {
693       G4cout << fGamma << "\t" << angleSum << G4endl;
694     }
695     iPlace = iTkin;
696     fAngleDistrTable->insertAt(iPlace, angleVector);
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 angle tables = "
704            << timer.GetUserElapsed() << " s" << G4endl;
705   }
706   fGamma = 0.;
707 
708   return;
709 }
710 
711 //////////////////////////////////////////////////////////////////////////////
712 // The main function which is responsible for the treatment of a particle
713 // passage through G4Envelope with discrete generation of G4Gamma
714 G4VParticleChange* G4VXTRenergyLoss::PostStepDoIt(const G4Track& aTrack,
715                                                   const G4Step& aStep)
716 {
717   G4int iTkin;
718   G4double energyTR, theta, theta2, phi, dirX, dirY, dirZ;
719 
720   fParticleChange.Initialize(aTrack);
721 
722   if(verboseLevel > 1)
723   {
724     G4cout << "Start of G4VXTRenergyLoss::PostStepDoIt " << G4endl;
725     G4cout << "name of current material =  "
726            << aTrack.GetVolume()->GetLogicalVolume()->GetMaterial()->GetName()
727            << G4endl;
728   }
729   if(aTrack.GetVolume()->GetLogicalVolume() != fEnvelope)
730   {
731     if(verboseLevel > 0)
732     {
733       G4cout << "Go out from G4VXTRenergyLoss::PostStepDoIt: wrong volume "
734              << G4endl;
735     }
736     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
737   }
738   else
739   {
740     G4StepPoint* pPostStepPoint        = aStep.GetPostStepPoint();
741     const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
742 
743     // Now we are ready to Generate one TR photon
744     G4double kinEnergy = aParticle->GetKineticEnergy();
745     G4double mass      = aParticle->GetDefinition()->GetPDGMass();
746     G4double gamma     = 1.0 + kinEnergy / mass;
747 
748     if(verboseLevel > 1)
749     {
750       G4cout << "gamma = " << gamma << G4endl;
751     }
752     G4double massRatio           = proton_mass_c2 / mass;
753     G4double TkinScaled          = kinEnergy * massRatio;
754     G4ThreeVector position       = pPostStepPoint->GetPosition();
755     G4ParticleMomentum direction = aParticle->GetMomentumDirection();
756     G4double startTime           = pPostStepPoint->GetGlobalTime();
757 
758     for(iTkin = 0; iTkin < fTotBin; ++iTkin)
759     {
760       if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))
761         break;
762     }
763 
764     if(iTkin == 0)  // Tkin is too small, neglect of TR photon generation
765     {
766       if(verboseLevel > 0)
767       {
768         G4cout << "Go out from G4VXTRenergyLoss::PostStepDoIt:iTkin = " << iTkin
769                << G4endl;
770       }
771       return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
772     }
773     else  // general case: Tkin between two vectors of the material
774     {
775       fParticleChange.SetNumberOfSecondaries(1);
776 
777       energyTR = GetXTRrandomEnergy(TkinScaled, iTkin);
778 
779       if(verboseLevel > 1)
780       {
781         G4cout << "energyTR = " << energyTR / keV << " keV" << G4endl;
782       }
783       if(fAngleRadDistr)
784       {
785         theta2 = GetRandomAngle(energyTR, iTkin);
786         if(theta2 > 0.)
787           theta = std::sqrt(theta2);
788         else
789           theta = 0.;
790       }
791       else
792         theta = std::fabs(G4RandGauss::shoot(0.0, pi / gamma));
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, dirZ);
804       directionTR.rotateUz(direction);
805       directionTR.unit();
806 
807       auto aPhotonTR =
808         new G4DynamicParticle(G4Gamma::Gamma(), directionTR, energyTR);
809 
810       // A XTR photon is set on the particle track inside the radiator
811       // and is moved to the G4Envelope surface for standard X-ray TR models
812       // only. The case of fExitFlux=true
813 
814       if(fExitFlux)
815       {
816         const G4RotationMatrix* rotM =
817           pPostStepPoint->GetTouchable()->GetRotation();
818         G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
819         G4AffineTransform transform = G4AffineTransform(rotM, transl);
820         transform.Invert();
821         G4ThreeVector localP = transform.TransformPoint(position);
822         G4ThreeVector localV = transform.TransformAxis(directionTR);
823 
824         G4double distance =
825           fEnvelope->GetSolid()->DistanceToOut(localP, localV);
826         if(verboseLevel > 1)
827         {
828           G4cout << "distance to exit = " << distance / mm << " mm" << G4endl;
829         }
830         position += distance * directionTR;
831         startTime += distance / c_light;
832       }
833       G4Track* aSecondaryTrack = new G4Track(aPhotonTR, startTime, position);
834       aSecondaryTrack->SetTouchableHandle(
835         aStep.GetPostStepPoint()->GetTouchableHandle());
836       aSecondaryTrack->SetParentID(aTrack.GetTrackID());
837 
838       fParticleChange.AddSecondary(aSecondaryTrack);
839       fParticleChange.ProposeEnergy(kinEnergy);
840     }
841   }
842   return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
843 }
844 
845 ///////////////////////////////////////////////////////////////////////
846 // This function returns the spectral and angle density of TR quanta
847 // in X-ray energy region generated forward when a relativistic
848 // charged particle crosses interface between two materials.
849 // The high energy small theta approximation is applied.
850 // (matter1 -> matter2, or 2->1)
851 // varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta
852 G4complex G4VXTRenergyLoss::OneInterfaceXTRdEdx(G4double energy, G4double gamma,
853                                                 G4double varAngle)
854 {
855   G4complex Z1 = GetPlateComplexFZ(energy, gamma, varAngle);
856   G4complex Z2 = GetGasComplexFZ(energy, gamma, varAngle);
857 
858   G4complex zOut = (Z1 - Z2) * (Z1 - Z2) * (varAngle * energy / hbarc / hbarc);
859   return zOut;
860 }
861 
862 //////////////////////////////////////////////////////////////////////////////
863 // For photon energy distribution tables. Integrate first over angle
864 G4double G4VXTRenergyLoss::SpectralAngleXTRdEdx(G4double varAngle)
865 {
866   G4double result = GetStackFactor(fEnergy, fGamma, varAngle);
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(G4double energy)
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, 0.1, 0.2, 0.5, 1.0 };
881 
882   for(i = 0; i < iMax; ++i)
883     lim[i] *= fMaxThetaTR;
884 
885   G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)>
886     integral;
887 
888   fEnergy = energy;
889   {
890     for(i = 0; i < iMax - 1; ++i)
891     {
892       angleSum += integral.Legendre96(
893         this, &G4VXTRenergyLoss::SpectralAngleXTRdEdx, lim[i], lim[i + 1]);
894     }
895   }
896   return angleSum;
897 }
898 
899 //////////////////////////////////////////////////////////////////////////
900 // for photon angle distribution tables
901 G4double G4VXTRenergyLoss::AngleSpectralXTRdEdx(G4double energy)
902 {
903   G4double result = GetStackFactor(energy, fGamma, fVarAngle);
904   if(result < 0)
905     result = 0.0;
906   return result;
907 }
908 
909 ///////////////////////////////////////////////////////////////////////////
910 // The XTR angular distribution based on transparent regular radiator
911 G4double G4VXTRenergyLoss::AngleXTRdEdx(G4double varAngle)
912 {
913   G4double result;
914   G4double sum = 0., tmp1, tmp2, tmp = 0., cof1, cof2, cofMin, cofPHC, energy1,
915            energy2;
916   G4int k, kMax, kMin, i;
917 
918   cofPHC = twopi * hbarc;
919 
920   cof1 = (fPlateThick + fGasThick) * (1. / fGamma / fGamma + varAngle);
921   cof2 = fPlateThick * fSigma1 + fGasThick * fSigma2;
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 * cof2);
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 || energy1 < fTheMinEnergyTR)
944           continue;
945 
946         tmp1 =
947           (energy1 * energy1 * (1. / fGamma / fGamma + varAngle) + fSigma1) *
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 / fGamma + varAngle) + fSigma2);
955         tmp *= (tmp1 - tmp2) * (tmp1 - tmp2);
956         tmp1 = cof1 / (4. * hbarc) - cof2 / (4. * hbarc * energy1 * energy1);
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 || energy2 < fTheMinEnergyTR)
967           continue;
968 
969         tmp1 =
970           (energy2 * energy2 * (1. / fGamma / fGamma + varAngle) + fSigma1) *
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 / fGamma + varAngle) + fSigma2);
978         tmp *= (tmp1 - tmp2) * (tmp1 - tmp2);
979         tmp1 = cof1 / (4. * hbarc) - cof2 / (4. * hbarc * energy2 * energy2);
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 * varAngle;
991   result /= hbarc * hbarc;
992 
993   return result;
994 }
995 
996 //////////////////////////////////////////////////////////////////////
997 // Calculates formation zone for plates. Omega is energy !!!
998 G4double G4VXTRenergyLoss::GetPlateFormationZone(G4double omega, G4double gamma,
999                                                  G4double varAngle)
1000 {
1001   G4double cof, lambda;
1002   lambda = 1.0 / gamma / gamma + varAngle + fSigma1 / omega / omega;
1003   cof    = 2.0 * hbarc / omega / lambda;
1004   return cof;
1005 }
1006 
1007 //////////////////////////////////////////////////////////////////////
1008 // Calculates complex formation zone for plates. Omega is energy !!!
1009 G4complex G4VXTRenergyLoss::GetPlateComplexFZ(G4double omega, G4double gamma,
1010                                               G4double varAngle)
1011 {
1012   G4double cof, length, delta, real_v, image_v;
1013 
1014   length = 0.5 * GetPlateFormationZone(omega, gamma, varAngle);
1015   delta  = length * GetPlateLinearPhotoAbs(omega);
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 cross section coefficients for
1027 // plate material
1028 void G4VXTRenergyLoss::ComputePlatePhotoAbsCof()
1029 {
1030   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1031   const G4Material* mat                   = (*theMaterialTable)[fMatIndex1];
1032   fPlatePhotoAbsCof                       = mat->GetSandiaTable();
1033 
1034   return;
1035 }
1036 
1037 //////////////////////////////////////////////////////////////////////
1038 // Returns the value of linear photo absorption coefficient (in reciprocal
1039 // length) for plate for given energy of X-ray photon omega
1040 G4double G4VXTRenergyLoss::GetPlateLinearPhotoAbs(G4double omega)
1041 {
1042   G4double omega2, omega3, omega4;
1043 
1044   omega2 = omega * omega;
1045   omega3 = omega2 * omega;
1046   omega4 = omega2 * omega2;
1047 
1048   const G4double* SandiaCof = fPlatePhotoAbsCof->GetSandiaCofForMaterial(omega);
1049   G4double cross            = SandiaCof[0] / omega + SandiaCof[1] / omega2 +
1050                    SandiaCof[2] / omega3 + SandiaCof[3] / omega4;
1051   return cross;
1052 }
1053 
1054 //////////////////////////////////////////////////////////////////////
1055 // Calculates formation zone for gas. Omega is energy !!!
1056 G4double G4VXTRenergyLoss::GetGasFormationZone(G4double omega, G4double gamma,
1057                                                G4double varAngle)
1058 {
1059   G4double cof, lambda;
1060   lambda = 1.0 / gamma / gamma + varAngle + fSigma2 / omega / omega;
1061   cof    = 2.0 * hbarc / omega / lambda;
1062   return cof;
1063 }
1064 
1065 //////////////////////////////////////////////////////////////////////
1066 // Calculates complex formation zone for gas gaps. Omega is energy !!!
1067 G4complex G4VXTRenergyLoss::GetGasComplexFZ(G4double omega, G4double gamma,
1068                                             G4double varAngle)
1069 {
1070   G4double cof, length, delta, real_v, image_v;
1071 
1072   length = 0.5 * GetGasFormationZone(omega, gamma, varAngle);
1073   delta  = length * GetGasLinearPhotoAbs(omega);
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 cross section coefficients for
1085 // gas material
1086 void G4VXTRenergyLoss::ComputeGasPhotoAbsCof()
1087 {
1088   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1089   const G4Material* mat                   = (*theMaterialTable)[fMatIndex2];
1090   fGasPhotoAbsCof                         = mat->GetSandiaTable();
1091   return;
1092 }
1093 
1094 //////////////////////////////////////////////////////////////////////
1095 // Returns the value of linear photo absorption coefficient (in reciprocal
1096 // length) for gas
1097 G4double G4VXTRenergyLoss::GetGasLinearPhotoAbs(G4double omega)
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->GetSandiaCofForMaterial(omega);
1106   G4double cross            = SandiaCof[0] / omega + SandiaCof[1] / omega2 +
1107                    SandiaCof[2] / omega3 + SandiaCof[3] / omega4;
1108   return cross;
1109 }
1110 
1111 //////////////////////////////////////////////////////////////////////
1112 // Calculates the product of linear cof by formation zone for plate.
1113 // Omega is energy !!!
1114 G4double G4VXTRenergyLoss::GetPlateZmuProduct(G4double omega, G4double gamma,
1115                                               G4double varAngle)
1116 {
1117   return GetPlateFormationZone(omega, gamma, varAngle) *
1118          GetPlateLinearPhotoAbs(omega);
1119 }
1120 //////////////////////////////////////////////////////////////////////
1121 // Calculates the product of linear cof by formation zone for plate.
1122 // G4cout and output in file in some energy range.
1123 void G4VXTRenergyLoss::GetPlateZmuProduct()
1124 {
1125   std::ofstream outPlate("plateZmu.dat", std::ios::out);
1126   outPlate.setf(std::ios::scientific, std::ios::floatfield);
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 for plate" << G4endl;
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, gamma, varAngle) << "\t";
1140     if(verboseLevel > 0)
1141       outPlate << omega / keV << "\t\t"
1142                << GetPlateZmuProduct(omega, gamma, varAngle) << G4endl;
1143   }
1144   return;
1145 }
1146 
1147 //////////////////////////////////////////////////////////////////////
1148 // Calculates the product of linear cof by formation zone for gas.
1149 // Omega is energy !!!
1150 G4double G4VXTRenergyLoss::GetGasZmuProduct(G4double omega, G4double gamma,
1151                                             G4double varAngle)
1152 {
1153   return GetGasFormationZone(omega, gamma, varAngle) *
1154          GetGasLinearPhotoAbs(omega);
1155 }
1156 
1157 //////////////////////////////////////////////////////////////////////
1158 // Calculates the product of linear cof by formation zone for gas.
1159 // G4cout and output in file in some energy range.
1160 void G4VXTRenergyLoss::GetGasZmuProduct()
1161 {
1162   std::ofstream outGas("gasZmu.dat", std::ios::out);
1163   outGas.setf(std::ios::scientific, std::ios::floatfield);
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 for gas" << G4endl;
1170   for(i = 0; i < 100; ++i)
1171   {
1172     omega = (1.0 + i) * keV;
1173     if(verboseLevel > 1)
1174       G4cout << omega / keV << "\t" << GetGasZmuProduct(omega, gamma, varAngle)
1175              << "\t";
1176     if(verboseLevel > 0)
1177       outGas << omega / keV << "\t\t"
1178              << GetGasZmuProduct(omega, gamma, varAngle) << G4endl;
1179   }
1180   return;
1181 }
1182 
1183 ////////////////////////////////////////////////////////////////////////
1184 // Computes Compton cross section for plate material in 1/mm
1185 G4double G4VXTRenergyLoss::GetPlateCompton(G4double omega)
1186 {
1187   G4int i, numberOfElements;
1188   G4double xSection = 0., nowZ, sumZ = 0.;
1189 
1190   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1191   numberOfElements = (G4int)(*theMaterialTable)[fMatIndex1]->GetNumberOfElements();
1192 
1193   for(i = 0; i < numberOfElements; ++i)
1194   {
1195     nowZ = (*theMaterialTable)[fMatIndex1]->GetElement(i)->GetZ();
1196     sumZ += nowZ;
1197     xSection += GetComptonPerAtom(omega, nowZ);
1198   }
1199   xSection /= sumZ;
1200   xSection *= (*theMaterialTable)[fMatIndex1]->GetElectronDensity();
1201   return xSection;
1202 }
1203 
1204 ////////////////////////////////////////////////////////////////////////
1205 // Computes Compton cross section for gas material in 1/mm
1206 G4double G4VXTRenergyLoss::GetGasCompton(G4double omega)
1207 {
1208   G4double xSection = 0., sumZ = 0.;
1209 
1210   const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1211   G4int numberOfElements = (G4int)(*theMaterialTable)[fMatIndex2]->GetNumberOfElements();
1212 
1213   for (G4int i = 0; i < numberOfElements; ++i)
1214   {
1215     G4double nowZ = (*theMaterialTable)[fMatIndex2]->GetElement(i)->GetZ();
1216     sumZ += nowZ;
1217     xSection += GetComptonPerAtom(omega, nowZ);
1218   }
1219   if (sumZ > 0.0) { xSection /= sumZ; }
1220   xSection *= (*theMaterialTable)[fMatIndex2]->GetElectronDensity();
1221   return xSection;
1222 }
1223 
1224 ////////////////////////////////////////////////////////////////////////
1225 // Computes Compton cross section per atom with Z electrons for gamma with
1226 // the energy GammaEnergy
1227 G4double G4VXTRenergyLoss::GetComptonPerAtom(G4double GammaEnergy, G4double Z)
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 * barn, d2 = -1.8300e-1 * barn,
1242                             d3 = 6.7527 * barn, d4 = -1.9798e+1 * barn,
1243                             e1 = 1.9756e-5 * barn, e2 = -1.0205e-2 * barn,
1244                             e3 = -7.3913e-2 * barn, e4 = 2.7079e-2 * barn,
1245                             f1 = -3.9178e-7 * barn, f2 = 6.8241e-5 * barn,
1246                             f3 = 6.0480e-5 * barn, f4 = 3.0274e-4 * barn;
1247 
1248   G4double p1Z = Z * (d1 + e1 * Z + f1 * Z * Z);
1249   G4double p2Z = Z * (d2 + e2 * Z + f2 * Z * Z);
1250   G4double p3Z = Z * (d3 + e3 * Z + f3 * Z * Z);
1251   G4double p4Z = Z * (d4 + e4 * Z + f4 * Z * 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) / electron_mass_c2;
1258   CrossSection =
1259     p1Z * std::log(1. + 2. * X) / X +
1260     (p2Z + p3Z * X + p4Z * X * X) / (1. + a * X + b * X * X + c * X * X * X);
1261 
1262   //  modification for low energy. (special case for Hydrogen)
1263   if(GammaEnergy < T0)
1264   {
1265     G4double dT0 = 1. * keV;
1266     X            = (T0 + dT0) / electron_mass_c2;
1267     G4double sigma =
1268       p1Z * std::log(1. + 2. * X) / X +
1269       (p2Z + p3Z * X + p4Z * X * X) / (1. + a * X + b * X * X + c * X * X * X);
1270     G4double c1 = -T0 * (sigma - CrossSection) / (CrossSection * dT0);
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 * y));
1276   }
1277   return CrossSection;
1278 }
1279 
1280 ///////////////////////////////////////////////////////////////////////
1281 // This function returns the spectral and angle density of TR quanta
1282 // in X-ray energy region generated forward when a relativistic
1283 // charged particle crosses interface between two materials.
1284 // The high energy small theta approximation is applied.
1285 // (matter1 -> matter2, or 2->1)
1286 // varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta
1287 G4double G4VXTRenergyLoss::OneBoundaryXTRNdensity(G4double energy,
1288                                                   G4double gamma,
1289                                                   G4double varAngle) const
1290 {
1291   G4double formationLength1, formationLength2;
1292   formationLength1 =
1293     1.0 / (1.0 / (gamma * gamma) + fSigma1 / (energy * energy) + varAngle);
1294   formationLength2 =
1295     1.0 / (1.0 / (gamma * gamma) + fSigma2 / (energy * energy) + varAngle);
1296   return (varAngle / energy) * (formationLength1 - formationLength2) *
1297          (formationLength1 - formationLength2);
1298 }
1299 
1300 G4double G4VXTRenergyLoss::GetStackFactor(G4double energy, G4double gamma,
1301                                           G4double varAngle)
1302 {
1303   // return stack factor corresponding to one interface
1304   return std::real(OneInterfaceXTRdEdx(energy, gamma, varAngle));
1305 }
1306 
1307 //////////////////////////////////////////////////////////////////////////////
1308 // For photon energy distribution tables. Integrate first over angle
1309 G4double G4VXTRenergyLoss::XTRNSpectralAngleDensity(G4double varAngle)
1310 {
1311   return OneBoundaryXTRNdensity(fEnergy, fGamma, varAngle) *
1312          GetStackFactor(fEnergy, fGamma, varAngle);
1313 }
1314 
1315 /////////////////////////////////////////////////////////////////////////
1316 // For second integration over energy
1317 G4double G4VXTRenergyLoss::XTRNSpectralDensity(G4double energy)
1318 {
1319   fEnergy = energy;
1320   G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)>
1321     integral;
1322   return integral.Legendre96(this, &G4VXTRenergyLoss::XTRNSpectralAngleDensity,
1323                              0.0, 0.2 * fMaxThetaTR) +
1324          integral.Legendre10(this, &G4VXTRenergyLoss::XTRNSpectralAngleDensity,
1325                              0.2 * fMaxThetaTR, fMaxThetaTR);
1326 }
1327 
1328 //////////////////////////////////////////////////////////////////////////
1329 // for photon angle distribution tables
1330 G4double G4VXTRenergyLoss::XTRNAngleSpectralDensity(G4double energy)
1331 {
1332   return OneBoundaryXTRNdensity(energy, fGamma, fVarAngle) *
1333          GetStackFactor(energy, fGamma, fVarAngle);
1334 }
1335 
1336 ///////////////////////////////////////////////////////////////////////////
1337 G4double G4VXTRenergyLoss::XTRNAngleDensity(G4double varAngle)
1338 {
1339   fVarAngle = varAngle;
1340   G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)>
1341     integral;
1342   return integral.Legendre96(this, &G4VXTRenergyLoss::XTRNAngleSpectralDensity,
1343                              fMinEnergyTR, fMaxEnergyTR);
1344 }
1345 
1346 //////////////////////////////////////////////////////////////////////////////
1347 // Check number of photons for a range of Lorentz factors from both energy
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::out);
1355   outEn.setf(std::ios::scientific, std::ios::floatfield);
1356 
1357   std::ofstream outAng("numberAng.dat", std::ios::out);
1358   outAng.setf(std::ios::scientific, std::ios::floatfield);
1359 
1360   for(iTkin = 0; iTkin < fTotBin; ++iTkin)  // Lorentz factor loop
1361   {
1362     gamma =
1363       1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2);
1364     numberE = (*(*fEnergyDistrTable)(iTkin))(0);
1365     if(verboseLevel > 1)
1366       G4cout << gamma << "\t\t" << numberE << "\t" << G4endl;
1367     if(verboseLevel > 0)
1368       outEn << gamma << "\t\t" << numberE << G4endl;
1369   }
1370   return;
1371 }
1372 
1373 /////////////////////////////////////////////////////////////////////////
1374 // Returns random energy of a X-ray TR photon for given scaled kinetic energy
1375 // of a charged particle
1376 G4double G4VXTRenergyLoss::GetXTRrandomEnergy(G4double scaledTkin, G4int iTkin)
1377 {
1378   G4int iTransfer, iPlace;
1379   G4double transfer = 0.0, position, E1, E2, W1, W2, W;
1380 
1381   iPlace = iTkin - 1;
1382 
1383   if(iTkin == fTotBin)  // relativistic plato, try from left
1384   {
1385     position = (*(*fEnergyDistrTable)(iPlace))(0) * G4UniformRand();
1386 
1387     for(iTransfer = 0;; ++iTransfer)
1388     {
1389       if(position >= (*(*fEnergyDistrTable)(iPlace))(iTransfer))
1390         break;
1391     }
1392     transfer = GetXTRenergy(iPlace, position, iTransfer);
1393   }
1394   else
1395   {
1396     E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
1397     E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin);
1398     W  = 1.0 / (E2 - E1);
1399     W1 = (E2 - scaledTkin) * W;
1400     W2 = (scaledTkin - E1) * W;
1401 
1402     position = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 +
1403                 (*(*fEnergyDistrTable)(iPlace + 1))(0) * W2) *
1404                G4UniformRand();
1405 
1406     for(iTransfer = 0;; ++iTransfer)
1407     {
1408       if(position >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer) *W1 +
1409                       (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer) *W2))
1410         break;
1411     }
1412     transfer = GetXTRenergy(iPlace, position, iTransfer);
1413   }
1414   if(transfer < 0.0)
1415     transfer = 0.0;
1416   return transfer;
1417 }
1418 
1419 ////////////////////////////////////////////////////////////////////////
1420 // Returns approximate position of X-ray photon energy during random sampling
1421 // over integral energy distribution
1422 G4double G4VXTRenergyLoss::GetXTRenergy(G4int iPlace, G4double, G4int iTransfer)
1423 {
1424   G4double x1, x2, y1, y2, result;
1425 
1426   if(iTransfer == 0)
1427   {
1428     result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1429   }
1430   else
1431   {
1432     y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer - 1);
1433     y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer);
1434 
1435     x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1);
1436     x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1437 
1438     if(x1 == x2)
1439       result = x2;
1440     else
1441     {
1442       if(y1 == y2)
1443         result = x1 + (x2 - x1) * G4UniformRand();
1444       else
1445       {
1446         result = x1 + (x2 - x1) * G4UniformRand();
1447       }
1448     }
1449   }
1450   return result;
1451 }
1452 
1453 /////////////////////////////////////////////////////////////////////////
1454 //  Get XTR photon angle at given energy and Tkin
1455 
1456 G4double G4VXTRenergyLoss::GetRandomAngle(G4double energyXTR, G4int iTkin)
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->GetLowEdgeEnergy(iTR))
1469       break;
1470   }
1471   if(iTR == fBinTR)
1472     --iTR;
1473 
1474   position = (*(*fAngleForEnergyTable)(iTR))(0) * G4UniformRand();
1475   // position = (*(*fAngleForEnergyTable)(iTR))(1) * G4UniformRand(); // ATLAS TB
1476 
1477   for(iAngle = 0;; ++iAngle)
1478   // for(iAngle = 1;; ++iAngle) // ATLAS TB
1479   {
1480     if(position >= (*(*fAngleForEnergyTable)(iTR))(iAngle))
1481       break;
1482   }
1483   angle = GetAngleXTR(iTR, position, iAngle);
1484   return angle;
1485 }
1486 
1487 ////////////////////////////////////////////////////////////////////////
1488 // Returns approximate position of X-ray photon angle at given energy during
1489 // random sampling over integral energy distribution
1490 
1491 G4double G4VXTRenergyLoss::GetAngleXTR(G4int iPlace, G4double position,
1492                                        G4int iTransfer)
1493 {
1494   G4double x1, x2, y1, y2, result;
1495 
1496   if( iTransfer == 0 )
1497   // if( iTransfer == 1 ) // ATLAS TB
1498   {
1499     result = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1500   }
1501   else
1502   {
1503     y1 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer - 1);
1504     y2 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer);
1505 
1506     x1 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1);
1507     x2 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
1508 
1509     if(x1 == x2) result = x2;
1510     else
1511     {
1512       if( y1 == y2 )  result = x1 + (x2 - x1) * G4UniformRand();
1513       else
1514       {
1515         result = x1 + (position - y1) * (x2 - x1) / (y2 - y1);
1516         // result = x1 + 0.1*(position - y1) * (x2 - x1) / (y2 - y1); // ATLAS TB
1517         // result = x1 + 0.05*(position - y1) * (x2 - x1) / (y2 - y1); // ATLAS TB
1518       }
1519     }
1520   }
1521   return result;
1522 }
1523