| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26
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::G4AdjointBremsstrahlungModel(G4VEmModel* aModel)
43 : G4VEmAdjointModel("AdjointeBremModel")
44 {
45 fDirectModel = aModel;
46 Initialize();
47 }
48
49 ////////////////////////////////////////////////////////////////////////////////
50 G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel()
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 G4EmModelManager();
64 fEmModelManagerForFwdModels->AddEmModel(1, fDirectModel, nullptr, nullptr);
65 SetApplyCutInRange(true);
66
67 fElectron = G4Electron::Electron();
68 fGamma = G4Gamma::Gamma();
69
70 fAdjEquivDirectPrimPart = G4AdjointElectron::AdjointElectron();
71 fAdjEquivDirectSecondPart = G4AdjointGamma::AdjointGamma();
72 fDirectPrimaryPart = fElectron;
73 fSecondPartSameType = false;
74
75 fCSManager = G4AdjointCSManager::GetAdjointCSManager();
76 }
77
78 ////////////////////////////////////////////////////////////////////////////////
79 G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel()
80 {
81 if(fEmModelManagerForFwdModels)
82 delete fEmModelManagerForFwdModels;
83 }
84
85 ////////////////////////////////////////////////////////////////////////////////
86 void G4AdjointBremsstrahlungModel::SampleSecondaries(
87 const G4Track& aTrack, G4bool isScatProjToProj,
88 G4ParticleChange* fParticleChange)
89 {
90 if(!fUseMatrix)
91 return RapidSampleSecondaries(aTrack, isScatProjToProj, fParticleChange);
92
93 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
94 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
95
96 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
97 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
98
99 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
100 {
101 return;
102 }
103
104 G4double projectileKinEnergy =
105 SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, isScatProjToProj);
106
107 // Weight correction
108 CorrectPostStepWeight(fParticleChange, aTrack.GetWeight(),
109 adjointPrimKinEnergy, projectileKinEnergy,
110 isScatProjToProj);
111
112 // Kinematic
113 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass();
114 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
115 G4double projectileP2 =
116 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
117 G4double projectileP = std::sqrt(projectileP2);
118
119 // Angle of the gamma direction with the projectile taken from
120 // G4eBremsstrahlungModel
121 G4double u;
122 if(0.25 > G4UniformRand())
123 u = -std::log(G4UniformRand() * G4UniformRand()) / 0.625;
124 else
125 u = -std::log(G4UniformRand() * G4UniformRand()) / 1.875;
126
127 G4double theta = u * electron_mass_c2 / projectileTotalEnergy;
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::sin(phi) * sint, cost) *
135 projectileP; // gamma frame
136 if(isScatProjToProj)
137 { // the adjoint primary is the scattered e-
138 G4ThreeVector gammaMomentum =
139 (projectileTotalEnergy - adjointPrimTotalEnergy) *
140 G4ThreeVector(0., 0., 1.);
141 G4ThreeVector dirProd = projectileMomentum - gammaMomentum;
142 G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
143 G4double sint1 = std::sqrt(1. - cost1 * cost1);
144 projectileMomentum =
145 G4ThreeVector(std::cos(phi) * sint1, std::sin(phi) * sint1, cost1) *
146 projectileP;
147 }
148
149 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
150
151 if(!isScatProjToProj)
152 { // kill the primary and add a secondary
153 fParticleChange->ProposeTrackStatus(fStopAndKill);
154 fParticleChange->AddSecondary(
155 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
156 }
157 else
158 {
159 fParticleChange->ProposeEnergy(projectileKinEnergy);
160 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
161 }
162 }
163
164 ////////////////////////////////////////////////////////////////////////////////
165 void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(
166 const G4Track& aTrack, G4bool isScatProjToProj,
167 G4ParticleChange* fParticleChange)
168 {
169 const G4DynamicParticle* theAdjointPrimary = aTrack.GetDynamicParticle();
170 DefineCurrentMaterial(aTrack.GetMaterialCutsCouple());
171
172 G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy();
173 G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy();
174
175 if(adjointPrimKinEnergy > GetHighEnergyLimit() * 0.999)
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 = GetSecondAdjEnergyMaxForProdToProj(adjointPrimKinEnergy);
187 G4double Emin = GetSecondAdjEnergyMinForProdToProj(adjointPrimKinEnergy);
188 if(Emin >= Emax)
189 return;
190 projectileKinEnergy = Emin * std::pow(Emax / Emin, G4UniformRand());
191 diffCSUsed = fCsBiasingFactor * fLastCZ / projectileKinEnergy;
192 }
193 else
194 {
195 G4double Emax =
196 GetSecondAdjEnergyMaxForScatProjToProj(adjointPrimKinEnergy);
197 G4double Emin =
198 GetSecondAdjEnergyMinForScatProjToProj(adjointPrimKinEnergy, fTcutSecond);
199 if(Emin >= Emax)
200 return;
201 G4double f1 = (Emin - adjointPrimKinEnergy) / Emin;
202 G4double f2 = (Emax - adjointPrimKinEnergy) / Emax / f1;
203 projectileKinEnergy =
204 adjointPrimKinEnergy / (1. - f1 * std::pow(f2, G4UniformRand()));
205 gammaEnergy = projectileKinEnergy - adjointPrimKinEnergy;
206 diffCSUsed =
207 fLastCZ * adjointPrimKinEnergy / projectileKinEnergy / gammaEnergy;
208 }
209
210 // Weight correction:
211 // First w_corr is set to the ratio between adjoint total CS and fwd total CS
212 // if this has to be done in the model.
213 // For the case of forced interaction this will be done in the PostStepDoIt of
214 // the forced interaction. It is important to set the weight before the
215 // creation of the secondary
216 G4double w_corr = fOutsideWeightFactor;
217 if(fInModelWeightCorr)
218 {
219 w_corr = fCSManager->GetPostStepWeightCorrection();
220 }
221
222 // Then another correction is needed due to the fact that a biaised
223 // differential CS has been used rather than the one consistent with the
224 // direct model Here we consider the true diffCS as the one obtained by the
225 // numerical differentiation over Tcut of the direct CS, corrected by the
226 // Migdal term. Basically any other differential CS could be used here
227 // (example Penelope).
228 G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(
229 fCurrentMaterial, projectileKinEnergy, gammaEnergy);
230 w_corr *= diffCS / diffCSUsed;
231
232 G4double new_weight = aTrack.GetWeight() * w_corr;
233 fParticleChange->SetParentWeightByProcess(false);
234 fParticleChange->SetSecondaryWeightByProcess(false);
235 fParticleChange->ProposeParentWeight(new_weight);
236
237 // Kinematic
238 G4double projectileM0 = fAdjEquivDirectPrimPart->GetPDGMass();
239 G4double projectileTotalEnergy = projectileM0 + projectileKinEnergy;
240 G4double projectileP2 =
241 projectileTotalEnergy * projectileTotalEnergy - projectileM0 * projectileM0;
242 G4double projectileP = std::sqrt(projectileP2);
243
244 // Use the angular model of the forward model to generate the gamma direction
245 // Dummy dynamic particle to use the model
246 G4DynamicParticle* aDynPart =
247 new G4DynamicParticle(fElectron, G4ThreeVector(0., 0., 1.) * projectileP);
248
249 // Get the element from the direct model
250 const G4Element* elm = fDirectModel->SelectRandomAtom(
251 fCurrentCouple, fElectron, projectileKinEnergy, fTcutSecond);
252 G4int Z = elm->GetZasInt();
253 G4double energy = aDynPart->GetTotalEnergy() - gammaEnergy;
254 G4ThreeVector projectileMomentum =
255 fDirectModel->GetAngularDistribution()->SampleDirection(aDynPart, energy, Z,
256 fCurrentMaterial) * projectileP;
257 G4double phi = projectileMomentum.getPhi();
258
259 if(isScatProjToProj)
260 { // the adjoint primary is the scattered e-
261 G4ThreeVector gammaMomentum =
262 (projectileTotalEnergy - adjointPrimTotalEnergy) *
263 G4ThreeVector(0., 0., 1.);
264 G4ThreeVector dirProd = projectileMomentum - gammaMomentum;
265 G4double cost1 = std::cos(dirProd.angle(projectileMomentum));
266 G4double sint1 = std::sqrt(1. - cost1 * cost1);
267 projectileMomentum =
268 G4ThreeVector(std::cos(phi) * sint1, std::sin(phi) * sint1, cost1) *
269 projectileP;
270 }
271
272 projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection());
273
274 if(!isScatProjToProj)
275 { // kill the primary and add a secondary
276 fParticleChange->ProposeTrackStatus(fStopAndKill);
277 fParticleChange->AddSecondary(
278 new G4DynamicParticle(fAdjEquivDirectPrimPart, projectileMomentum));
279 }
280 else
281 {
282 fParticleChange->ProposeEnergy(projectileKinEnergy);
283 fParticleChange->ProposeMomentumDirection(projectileMomentum.unit());
284 }
285 }
286
287 ////////////////////////////////////////////////////////////////////////////////
288 G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(
289 const G4Material* aMaterial,
290 G4double kinEnergyProj, // kin energy of primary before interaction
291 G4double kinEnergyProd // kinetic energy of the secondary particle
292 )
293 {
294 if(!fIsDirectModelInitialised)
295 {
296 fEmModelManagerForFwdModels->Initialise(fElectron, fGamma, 0);
297 fIsDirectModelInitialised = true;
298 }
299 return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
300 aMaterial, kinEnergyProj, kinEnergyProd);
301 }
302
303 ////////////////////////////////////////////////////////////////////////////////
304 G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(
305 const G4MaterialCutsCouple* aCouple, G4double primEnergy,
306 G4bool isScatProjToProj)
307 {
308 static constexpr G4double maxEnergy = 100. * MeV / 2.718281828459045;
309 // 2.78.. == std::exp(1.)
310 if(!fIsDirectModelInitialised)
311 {
312 fEmModelManagerForFwdModels->Initialise(fElectron, fGamma, 0);
313 fIsDirectModelInitialised = true;
314 }
315 if(fUseMatrix)
316 return G4VEmAdjointModel::AdjointCrossSection(aCouple, primEnergy,
317 isScatProjToProj);
318 DefineCurrentMaterial(aCouple);
319 G4double Cross = 0.;
320 // this gives the constant above
321 fLastCZ = fDirectModel->CrossSectionPerVolume(
322 aCouple->GetMaterial(), fDirectPrimaryPart, 100. * MeV, maxEnergy);
323
324 if(!isScatProjToProj)
325 {
326 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(primEnergy);
327 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(primEnergy);
328 if(Emax_proj > Emin_proj && primEnergy > fTcutSecond)
329 Cross = fCsBiasingFactor * fLastCZ * std::log(Emax_proj / Emin_proj);
330 }
331 else
332 {
333 G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProj(primEnergy);
334 G4double Emin_proj =
335 GetSecondAdjEnergyMinForScatProjToProj(primEnergy, fTcutSecond);
336 if(Emax_proj > Emin_proj)
337 Cross = fLastCZ * std::log((Emax_proj - primEnergy) * Emin_proj /
338 Emax_proj / (Emin_proj - primEnergy));
339 }
340 return Cross;
341 }
342