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 ]

Diff markup

Differences between /processes/electromagnetic/xrays/src/G4VXTRenergyLoss.cc (Version 11.3.0) and /processes/electromagnetic/xrays/src/G4VXTRenergyLoss.cc (Version 6.2.p2)


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