| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26
27 #include "G4AdjointBremsstrahlungModel.hh"
28
29 #include "G4AdjointCSManager.hh"
30 #include "G4AdjointElectron.hh"
31 #include "G4AdjointGamma.hh"
32 #include "G4Electron.hh"
33 #include "G4EmModelManager.hh"
34 #include "G4Gamma.hh"
35 #include "G4ParticleChange.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4SeltzerBergerModel.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4TrackStatus.hh"
40
41 //////////////////////////////////////////////
42 G4AdjointBremsstrahlungModel::G4AdjointBremsst
43 : G4VEmAdjointModel("AdjointeBremModel")
44 {
45 fDirectModel = aModel;
46 Initialize();
47 }
48
49 //////////////////////////////////////////////
50 G4AdjointBremsstrahlungModel::G4AdjointBremsst
51 : G4VEmAdjointModel("AdjointeBremModel")
52 {
53 fDirectModel = new G4SeltzerBergerModel();
54 Initialize();
55 }
56
57 //////////////////////////////////////////////
58 void G4AdjointBremsstrahlungModel::Initialize(
59 {
60 SetUseMatrix(false);
61 SetUseMatrixPerElement(false);
62
63 fEmModelManagerForFwdModels = new G4EmModelM
64 fEmModelManagerForFwdModels->AddEmModel(1, f
65 SetApplyCutInRange(true);
66
67 fElectron = G4Electron::Electron();
68 fGamma = G4Gamma::Gamma();
69
70 fAdjEquivDirectPrimPart = G4AdjointElectro
71 fAdjEquivDirectSecondPart = G4AdjointGamma::
72 fDirectPrimaryPart = fElectron;
73 fSecondPartSameType = false;
74
75 fCSManager = G4AdjointCSManager::GetAdjointC
76 }
77
78 //////////////////////////////////////////////
79 G4AdjointBremsstrahlungModel::~G4AdjointBremss
80 {
81 if(fEmModelManagerForFwdModels)
82 delete fEmModelManagerForFwdModels;
83 }
84
85 //////////////////////////////////////////////
86 void G4AdjointBremsstrahlungModel::SampleSecon
87 const G4Track& aTrack, G4bool isScatProjToPr
88 G4ParticleChange* fParticleChange)
89 {
90 if(!fUseMatrix)
91 return RapidSampleSecondaries(aTrack, isSc
92
93 const G4DynamicParticle* theAdjointPrimary =
94 DefineCurrentMaterial(aTrack.GetMaterialCuts
95
96 G4double adjointPrimKinEnergy = theAdjoint
97 G4double adjointPrimTotalEnergy = theAdjoint
98
99 if(adjointPrimKinEnergy > GetHighEnergyLimit
100 {
101 return;
102 }
103
104 G4double projectileKinEnergy =
105 SampleAdjSecEnergyFromCSMatrix(adjointPrim
106
107 // Weight correction
108 CorrectPostStepWeight(fParticleChange, aTrac
109 adjointPrimKinEnergy,
110 isScatProjToProj);
111
112 // Kinematic
113 G4double projectileM0 = fAdjEquivDi
114 G4double projectileTotalEnergy = projectileM
115 G4double projectileP2 =
116 projectileTotalEnergy * projectileTotalEne
117 G4double projectileP = std::sqrt(projectileP
118
119 // Angle of the gamma direction with the pro
120 // G4eBremsstrahlungModel
121 G4double u;
122 if(0.25 > G4UniformRand())
123 u = -std::log(G4UniformRand() * G4UniformR
124 else
125 u = -std::log(G4UniformRand() * G4UniformR
126
127 G4double theta = u * electron_mass_c2 / proj
128 G4double sint = std::sin(theta);
129 G4double cost = std::cos(theta);
130
131 G4double phi = twopi * G4UniformRand();
132
133 G4ThreeVector projectileMomentum =
134 G4ThreeVector(std::cos(phi) * sint, std::s
135 projectileP; // gamma frame
136 if(isScatProjToProj)
137 { // the adjoint primary is the scattered e
138 G4ThreeVector gammaMomentum =
139 (projectileTotalEnergy - adjointPrimTota
140 G4ThreeVector(0., 0., 1.);
141 G4ThreeVector dirProd = projectileMomentum
142 G4double cost1 = std::cos(dirProd.a
143 G4double sint1 = std::sqrt(1. - cos
144 projectileMomentum =
145 G4ThreeVector(std::cos(phi) * sint1, std
146 projectileP;
147 }
148
149 projectileMomentum.rotateUz(theAdjointPrimar
150
151 if(!isScatProjToProj)
152 { // kill the primary and add a secondary
153 fParticleChange->ProposeTrackStatus(fStopA
154 fParticleChange->AddSecondary(
155 new G4DynamicParticle(fAdjEquivDirectPri
156 }
157 else
158 {
159 fParticleChange->ProposeEnergy(projectileK
160 fParticleChange->ProposeMomentumDirection(
161 }
162 }
163
164 //////////////////////////////////////////////
165 void G4AdjointBremsstrahlungModel::RapidSample
166 const G4Track& aTrack, G4bool isScatProjToPr
167 G4ParticleChange* fParticleChange)
168 {
169 const G4DynamicParticle* theAdjointPrimary =
170 DefineCurrentMaterial(aTrack.GetMaterialCuts
171
172 G4double adjointPrimKinEnergy = theAdjoint
173 G4double adjointPrimTotalEnergy = theAdjoint
174
175 if(adjointPrimKinEnergy > GetHighEnergyLimit
176 {
177 return;
178 }
179
180 G4double projectileKinEnergy = 0.;
181 G4double gammaEnergy = 0.;
182 G4double diffCSUsed = 0.;
183 if(!isScatProjToProj)
184 {
185 gammaEnergy = adjointPrimKinEnergy;
186 G4double Emax = GetSecondAdjEnergyMaxForPr
187 G4double Emin = GetSecondAdjEnergyMinForPr
188 if(Emin >= Emax)
189 return;
190 projectileKinEnergy = Emin * std::pow(Emax
191 diffCSUsed = fCsBiasingFactor * f
192 }
193 else
194 {
195 G4double Emax =
196 GetSecondAdjEnergyMaxForScatProjToProj(a
197 G4double Emin =
198 GetSecondAdjEnergyMinForScatProjToProj(a
199 if(Emin >= Emax)
200 return;
201 G4double f1 = (Emin - adjointPrimKinEnergy
202 G4double f2 = (Emax - adjointPrimKinEnergy
203 projectileKinEnergy =
204 adjointPrimKinEnergy / (1. - f1 * std::p
205 gammaEnergy = projectileKinEnergy - adjoin
206 diffCSUsed =
207 fLastCZ * adjointPrimKinEnergy / project
208 }
209
210 // Weight correction:
211 // First w_corr is set to the ratio between
212 // if this has to be done in the model.
213 // For the case of forced interaction this w
214 // the forced interaction. It is important
215 // creation of the secondary
216 G4double w_corr = fOutsideWeightFactor;
217 if(fInModelWeightCorr)
218 {
219 w_corr = fCSManager->GetPostStepWeightCorr
220 }
221
222 // Then another correction is needed due to
223 // differential CS has been used rather than
224 // direct model Here we consider the true di
225 // numerical differentiation over Tcut of th
226 // Migdal term. Basically any other differen
227 // (example Penelope).
228 G4double diffCS = DiffCrossSectionPerVolumeP
229 fCurrentMaterial, projectileKinEnergy, gam
230 w_corr *= diffCS / diffCSUsed;
231
232 G4double new_weight = aTrack.GetWeight() * w
233 fParticleChange->SetParentWeightByProcess(fa
234 fParticleChange->SetSecondaryWeightByProcess
235 fParticleChange->ProposeParentWeight(new_wei
236
237 // Kinematic
238 G4double projectileM0 = fAdjEquivDi
239 G4double projectileTotalEnergy = projectileM
240 G4double projectileP2 =
241 projectileTotalEnergy * projectileTotalEne
242 G4double projectileP = std::sqrt(projectileP
243
244 // Use the angular model of the forward mode
245 // Dummy dynamic particle to use the model
246 G4DynamicParticle* aDynPart =
247 new G4DynamicParticle(fElectron, G4ThreeVe
248
249 // Get the element from the direct model
250 const G4Element* elm = fDirectModel->SelectR
251 fCurrentCouple, fElectron, projectileKinEn
252 G4int Z = elm->GetZasInt();
253 G4double energy = aDynPart->GetTotalEnergy()
254 G4ThreeVector projectileMomentum =
255 fDirectModel->GetAngularDistribution()->Sa
256 fC
257 G4double phi = projectileMomentum.getPhi();
258
259 if(isScatProjToProj)
260 { // the adjoint primary is the scattered e
261 G4ThreeVector gammaMomentum =
262 (projectileTotalEnergy - adjointPrimTota
263 G4ThreeVector(0., 0., 1.);
264 G4ThreeVector dirProd = projectileMomentum
265 G4double cost1 = std::cos(dirProd.a
266 G4double sint1 = std::sqrt(1. - cos
267 projectileMomentum =
268 G4ThreeVector(std::cos(phi) * sint1, std
269 projectileP;
270 }
271
272 projectileMomentum.rotateUz(theAdjointPrimar
273
274 if(!isScatProjToProj)
275 { // kill the primary and add a secondary
276 fParticleChange->ProposeTrackStatus(fStopA
277 fParticleChange->AddSecondary(
278 new G4DynamicParticle(fAdjEquivDirectPri
279 }
280 else
281 {
282 fParticleChange->ProposeEnergy(projectileK
283 fParticleChange->ProposeMomentumDirection(
284 }
285 }
286
287 //////////////////////////////////////////////
288 G4double G4AdjointBremsstrahlungModel::DiffCro
289 const G4Material* aMaterial,
290 G4double kinEnergyProj, // kin energy of pr
291 G4double kinEnergyProd // kinetic energy o
292 )
293 {
294 if(!fIsDirectModelInitialised)
295 {
296 fEmModelManagerForFwdModels->Initialise(fE
297 fIsDirectModelInitialised = true;
298 }
299 return G4VEmAdjointModel::DiffCrossSectionPe
300 aMaterial, kinEnergyProj, kinEnergyProd);
301 }
302
303 //////////////////////////////////////////////
304 G4double G4AdjointBremsstrahlungModel::Adjoint
305 const G4MaterialCutsCouple* aCouple, G4doubl
306 G4bool isScatProjToProj)
307 {
308 static constexpr G4double maxEnergy = 100. *
309 // 2.78.. == std::exp(1.)
310 if(!fIsDirectModelInitialised)
311 {
312 fEmModelManagerForFwdModels->Initialise(fE
313 fIsDirectModelInitialised = true;
314 }
315 if(fUseMatrix)
316 return G4VEmAdjointModel::AdjointCrossSect
317
318 DefineCurrentMaterial(aCouple);
319 G4double Cross = 0.;
320 // this gives the constant above
321 fLastCZ = fDirectModel->CrossSectionPerVolum
322 aCouple->GetMaterial(), fDirectPrimaryPart
323
324 if(!isScatProjToProj)
325 {
326 G4double Emax_proj = GetSecondAdjEnergyMax
327 G4double Emin_proj = GetSecondAdjEnergyMin
328 if(Emax_proj > Emin_proj && primEnergy > f
329 Cross = fCsBiasingFactor * fLastCZ * std
330 }
331 else
332 {
333 G4double Emax_proj = GetSecondAdjEnergyMax
334 G4double Emin_proj =
335 GetSecondAdjEnergyMinForScatProjToProj(p
336 if(Emax_proj > Emin_proj)
337 Cross = fLastCZ * std::log((Emax_proj -
338 Emax_proj / (
339 }
340 return Cross;
341 }
342