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