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