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