| 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 "G4VEmAdjointModel.hh" 27 #include "G4VEmAdjointModel.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 "G4AdjointInterpolator.hh" 32 #include "G4AdjointInterpolator.hh"
33 #include "G4AdjointPositron.hh" 33 #include "G4AdjointPositron.hh"
34 #include "G4Electron.hh" 34 #include "G4Electron.hh"
35 #include "G4Gamma.hh" 35 #include "G4Gamma.hh"
36 #include "G4Integrator.hh" 36 #include "G4Integrator.hh"
37 #include "G4ParticleChange.hh" 37 #include "G4ParticleChange.hh"
38 #include "G4ProductionCutsTable.hh" 38 #include "G4ProductionCutsTable.hh"
39 #include "G4TrackStatus.hh" 39 #include "G4TrackStatus.hh"
40 40
41 ////////////////////////////////////////////// 41 ////////////////////////////////////////////////////////////////////////////////
42 G4VEmAdjointModel::G4VEmAdjointModel(const G4S 42 G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam)
43 : fName(nam) 43 : fName(nam)
44 { 44 {
45 fCSManager = G4AdjointCSManager::GetAdjointC 45 fCSManager = G4AdjointCSManager::GetAdjointCSManager();
46 fCSManager->RegisterEmAdjointModel(this); 46 fCSManager->RegisterEmAdjointModel(this);
47 } 47 }
48 48
49 ///////////////////////////////////e////////// 49 ///////////////////////////////////e////////////////////////////////////////////
50 G4VEmAdjointModel::~G4VEmAdjointModel() 50 G4VEmAdjointModel::~G4VEmAdjointModel()
51 { 51 {
52 delete fCSMatrixProdToProjBackScat; 52 delete fCSMatrixProdToProjBackScat;
53 fCSMatrixProdToProjBackScat = nullptr; 53 fCSMatrixProdToProjBackScat = nullptr;
54 54
55 delete fCSMatrixProjToProjBackScat; 55 delete fCSMatrixProjToProjBackScat;
56 fCSMatrixProjToProjBackScat = nullptr; 56 fCSMatrixProjToProjBackScat = nullptr;
57 } 57 }
58 58
59 ////////////////////////////////////////////// 59 ////////////////////////////////////////////////////////////////////////////////
60 G4double G4VEmAdjointModel::AdjointCrossSectio 60 G4double G4VEmAdjointModel::AdjointCrossSection(
61 const G4MaterialCutsCouple* aCouple, G4doubl 61 const G4MaterialCutsCouple* aCouple, G4double primEnergy,
62 G4bool isScatProjToProj) 62 G4bool isScatProjToProj)
63 { 63 {
64 DefineCurrentMaterial(aCouple); 64 DefineCurrentMaterial(aCouple);
65 fPreStepEnergy = primEnergy; 65 fPreStepEnergy = primEnergy;
66 66
67 std::vector<G4double>* CS_Vs_Element = &fEle 67 std::vector<G4double>* CS_Vs_Element = &fElementCSProdToProj;
68 if(isScatProjToProj) 68 if(isScatProjToProj)
69 CS_Vs_Element = &fElementCSScatProjToProj; 69 CS_Vs_Element = &fElementCSScatProjToProj;
70 fLastCS = 70 fLastCS =
71 fCSManager->ComputeAdjointCS(fCurrentMater 71 fCSManager->ComputeAdjointCS(fCurrentMaterial, this, primEnergy,
72 fTcutSecond, 72 fTcutSecond, isScatProjToProj, *CS_Vs_Element);
73 if(isScatProjToProj) 73 if(isScatProjToProj)
74 fLastAdjointCSForScatProjToProj = fLastCS; 74 fLastAdjointCSForScatProjToProj = fLastCS;
75 else 75 else
76 fLastAdjointCSForProdToProj = fLastCS; 76 fLastAdjointCSForProdToProj = fLastCS;
77 77
78 return fLastCS; 78 return fLastCS;
79 } 79 }
80 80
81 ////////////////////////////////////////////// 81 ////////////////////////////////////////////////////////////////////////////////
82 G4double G4VEmAdjointModel::DiffCrossSectionPe 82 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
83 G4double kinEnergyProj, G4double kinEnergyPr 83 G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A)
84 { 84 {
85 G4double dSigmadEprod = 0.; 85 G4double dSigmadEprod = 0.;
86 G4double Emax_proj = GetSecondAdjEnergyMa 86 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
87 G4double Emin_proj = GetSecondAdjEnergyMi 87 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
88 88
89 // the produced particle should have a kinet 89 // the produced particle should have a kinetic energy less than the projectile
90 if(kinEnergyProj > Emin_proj && kinEnergyPro 90 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
91 { 91 {
92 G4double E1 = kinEnergyProd; 92 G4double E1 = kinEnergyProd;
93 G4double E2 = kinEnergyProd * 1.000001 93 G4double E2 = kinEnergyProd * 1.000001;
94 G4double sigma1 = fDirectModel->ComputeCro 94 G4double sigma1 = fDirectModel->ComputeCrossSectionPerAtom(
95 fDirectPrimaryPart, kinEnergyProj, Z, A, 95 fDirectPrimaryPart, kinEnergyProj, Z, A, E1, 1.e20);
96 G4double sigma2 = fDirectModel->ComputeCro 96 G4double sigma2 = fDirectModel->ComputeCrossSectionPerAtom(
97 fDirectPrimaryPart, kinEnergyProj, Z, A, 97 fDirectPrimaryPart, kinEnergyProj, Z, A, E2, 1.e20);
98 98
99 dSigmadEprod = (sigma1 - sigma2) / (E2 - E 99 dSigmadEprod = (sigma1 - sigma2) / (E2 - E1);
100 } 100 }
101 return dSigmadEprod; 101 return dSigmadEprod;
102 } 102 }
103 103
104 ////////////////////////////////////////////// 104 ////////////////////////////////////////////////////////////////////////////////
105 G4double G4VEmAdjointModel::DiffCrossSectionPe 105 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
106 G4double kinEnergyProj, G4double kinEnergySc 106 G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A)
107 { 107 {
108 G4double kinEnergyProd = kinEnergyProj - kin 108 G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
109 G4double dSigmadEprod; 109 G4double dSigmadEprod;
110 if(kinEnergyProd <= 0.) 110 if(kinEnergyProd <= 0.)
111 dSigmadEprod = 0.; 111 dSigmadEprod = 0.;
112 else 112 else
113 dSigmadEprod = 113 dSigmadEprod =
114 DiffCrossSectionPerAtomPrimToSecond(kinE 114 DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj, kinEnergyProd, Z, A);
115 return dSigmadEprod; 115 return dSigmadEprod;
116 } 116 }
117 117
118 ////////////////////////////////////////////// 118 ////////////////////////////////////////////////////////////////////////////////
119 G4double G4VEmAdjointModel::DiffCrossSectionPe 119 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
120 const G4Material* aMaterial, G4double kinEne 120 const G4Material* aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
121 { 121 {
122 G4double dSigmadEprod = 0.; 122 G4double dSigmadEprod = 0.;
123 G4double Emax_proj = GetSecondAdjEnergyMa 123 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
124 G4double Emin_proj = GetSecondAdjEnergyMi 124 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
125 125
126 if(kinEnergyProj > Emin_proj && kinEnergyPro 126 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
127 { 127 {
128 G4double E1 = kinEnergyProd; 128 G4double E1 = kinEnergyProd;
129 G4double E2 = kinEnergyProd * 1.0001; 129 G4double E2 = kinEnergyProd * 1.0001;
130 G4double sigma1 = fDirectModel->CrossSecti 130 G4double sigma1 = fDirectModel->CrossSectionPerVolume(
131 aMaterial, fDirectPrimaryPart, kinEnergy 131 aMaterial, fDirectPrimaryPart, kinEnergyProj, E1, 1.e20);
132 G4double sigma2 = fDirectModel->CrossSecti 132 G4double sigma2 = fDirectModel->CrossSectionPerVolume(
133 aMaterial, fDirectPrimaryPart, kinEnergy 133 aMaterial, fDirectPrimaryPart, kinEnergyProj, E2, 1.e20);
134 dSigmadEprod = (sigma1 - sigma2) / (E2 - E 134 dSigmadEprod = (sigma1 - sigma2) / (E2 - E1);
135 } 135 }
136 return dSigmadEprod; 136 return dSigmadEprod;
137 } 137 }
138 138
139 ////////////////////////////////////////////// 139 ////////////////////////////////////////////////////////////////////////////////
140 G4double G4VEmAdjointModel::DiffCrossSectionPe 140 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
141 const G4Material* aMaterial, G4double kinEne 141 const G4Material* aMaterial, G4double kinEnergyProj,
142 G4double kinEnergyScatProj) 142 G4double kinEnergyScatProj)
143 { 143 {
144 G4double kinEnergyProd = kinEnergyProj - kin 144 G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
145 G4double dSigmadEprod; 145 G4double dSigmadEprod;
146 if(kinEnergyProd <= 0.) 146 if(kinEnergyProd <= 0.)
147 dSigmadEprod = 0.; 147 dSigmadEprod = 0.;
148 else 148 else
149 dSigmadEprod = DiffCrossSectionPerVolumePr 149 dSigmadEprod = DiffCrossSectionPerVolumePrimToSecond(
150 aMaterial, kinEnergyProj, kinEnergyProd) 150 aMaterial, kinEnergyProj, kinEnergyProd);
151 return dSigmadEprod; 151 return dSigmadEprod;
152 } 152 }
153 153
154 ////////////////////////////////////////////// 154 ////////////////////////////////////////////////////////////////////////////////
155 G4double G4VEmAdjointModel::DiffCrossSectionFu 155 G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj)
156 { 156 {
157 G4double bias_factor = 157 G4double bias_factor =
158 fCsBiasingFactor * fKinEnergyProdForIntegr 158 fCsBiasingFactor * fKinEnergyProdForIntegration / kinEnergyProj;
159 159
160 if(fUseMatrixPerElement) 160 if(fUseMatrixPerElement)
161 { 161 {
162 return DiffCrossSectionPerAtomPrimToSecond 162 return DiffCrossSectionPerAtomPrimToSecond(
163 kinEnergyProj, fKinEnergyProdForI 163 kinEnergyProj, fKinEnergyProdForIntegration, fZSelectedNucleus,
164 fASelectedNucleus) * 164 fASelectedNucleus) *
165 bias_factor; 165 bias_factor;
166 } 166 }
167 else 167 else
168 { 168 {
169 return DiffCrossSectionPerVolumePrimToSeco 169 return DiffCrossSectionPerVolumePrimToSecond(
170 fSelectedMaterial, kinEnergyProj, 170 fSelectedMaterial, kinEnergyProj, fKinEnergyProdForIntegration) *
171 bias_factor; 171 bias_factor;
172 } 172 }
173 } 173 }
174 174
175 ////////////////////////////////////////////// 175 ////////////////////////////////////////////////////////////////////////////////
176 G4double G4VEmAdjointModel::DiffCrossSectionFu 176 G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj)
177 { 177 {
178 G4double bias_factor = 178 G4double bias_factor =
179 fCsBiasingFactor * fKinEnergyScatProjForIn 179 fCsBiasingFactor * fKinEnergyScatProjForIntegration / kinEnergyProj;
180 if(fUseMatrixPerElement) 180 if(fUseMatrixPerElement)
181 { 181 {
182 return DiffCrossSectionPerAtomPrimToScatPr 182 return DiffCrossSectionPerAtomPrimToScatPrim(
183 kinEnergyProj, fKinEnergyScatProj 183 kinEnergyProj, fKinEnergyScatProjForIntegration, fZSelectedNucleus,
184 fASelectedNucleus) * 184 fASelectedNucleus) *
185 bias_factor; 185 bias_factor;
186 } 186 }
187 else 187 else
188 { 188 {
189 return DiffCrossSectionPerVolumePrimToScat 189 return DiffCrossSectionPerVolumePrimToScatPrim(
190 fSelectedMaterial, kinEnergyProj, 190 fSelectedMaterial, kinEnergyProj,
191 fKinEnergyScatProjForIntegration) 191 fKinEnergyScatProjForIntegration) *
192 bias_factor; 192 bias_factor;
193 } 193 }
194 } 194 }
195 195
196 ////////////////////////////////////////////// 196 ////////////////////////////////////////////////////////////////////////////////
197 std::vector<std::vector<G4double>*> 197 std::vector<std::vector<G4double>*>
198 G4VEmAdjointModel::ComputeAdjointCrossSectionV 198 G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(
199 G4double kinEnergyProd, G4double Z, G4double 199 G4double kinEnergyProd, G4double Z, G4double A,
200 G4int nbin_pro_decade) // nb bins per order 200 G4int nbin_pro_decade) // nb bins per order of magnitude of energy
201 { 201 {
202 G4Integrator<G4VEmAdjointModel, G4double (G4 202 G4Integrator<G4VEmAdjointModel, G4double (G4VEmAdjointModel::*)(G4double)>
203 integral; 203 integral;
204 fASelectedNucleus = G4lrint(A); 204 fASelectedNucleus = G4lrint(A);
205 fZSelectedNucleus = G4lrint(Z); 205 fZSelectedNucleus = G4lrint(Z);
206 fKinEnergyProdForIntegration = kinEnergyProd 206 fKinEnergyProdForIntegration = kinEnergyProd;
207 207
208 // compute the vector of integrated cross se 208 // compute the vector of integrated cross sections
209 G4double minEProj = GetSecondAdjEnergyMinFor 209 G4double minEProj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
210 G4double maxEProj = GetSecondAdjEnergyMaxFor 210 G4double maxEProj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
211 G4double E1 = minEProj; 211 G4double E1 = minEProj;
212 std::vector<G4double>* log_ESec_vector = new 212 std::vector<G4double>* log_ESec_vector = new std::vector<G4double>();
213 std::vector<G4double>* log_Prob_vector = new 213 std::vector<G4double>* log_Prob_vector = new std::vector<G4double>();
214 log_ESec_vector->push_back(std::log(E1)); 214 log_ESec_vector->push_back(std::log(E1));
215 log_Prob_vector->push_back(-50.); 215 log_Prob_vector->push_back(-50.);
216 216
217 G4double E2 = 217 G4double E2 =
218 std::pow(10., G4double(G4int(std::log10(mi 218 std::pow(10., G4double(G4int(std::log10(minEProj) * nbin_pro_decade) + 1) /
219 nbin_pro_decade); 219 nbin_pro_decade);
220 G4double fE = std::pow(10., 1. / nbin_pro_de 220 G4double fE = std::pow(10., 1. / nbin_pro_decade);
221 221
222 if(std::pow(fE, 5.) > (maxEProj / minEProj)) 222 if(std::pow(fE, 5.) > (maxEProj / minEProj))
223 fE = std::pow(maxEProj / minEProj, 0.2); 223 fE = std::pow(maxEProj / minEProj, 0.2);
224 224
225 G4double int_cross_section = 0.; 225 G4double int_cross_section = 0.;
226 // Loop checking, 07-Aug-2015, Vladimir Ivan 226 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
227 while(E1 < maxEProj * 0.9999999) 227 while(E1 < maxEProj * 0.9999999)
228 { 228 {
229 int_cross_section += 229 int_cross_section +=
230 integral.Simpson(this, &G4VEmAdjointMode 230 integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1, E1,
231 std::min(E2, maxEProj * 231 std::min(E2, maxEProj * 0.99999999), 5);
232 log_ESec_vector->push_back(std::log(std::m 232 log_ESec_vector->push_back(std::log(std::min(E2, maxEProj)));
233 log_Prob_vector->push_back(std::log(int_cr 233 log_Prob_vector->push_back(std::log(int_cross_section));
234 E1 = E2; 234 E1 = E2;
235 E2 *= fE; 235 E2 *= fE;
236 } 236 }
237 std::vector<std::vector<G4double>*> res_mat; 237 std::vector<std::vector<G4double>*> res_mat;
238 if(int_cross_section > 0.) 238 if(int_cross_section > 0.)
239 { 239 {
240 res_mat.push_back(log_ESec_vector); 240 res_mat.push_back(log_ESec_vector);
241 res_mat.push_back(log_Prob_vector); 241 res_mat.push_back(log_Prob_vector);
242 } 242 }
243 else { 243 else {
244 delete log_ESec_vector; 244 delete log_ESec_vector;
245 delete log_Prob_vector; 245 delete log_Prob_vector;
246 } 246 }
247 return res_mat; 247 return res_mat;
248 } 248 }
249 249
250 ////////////////////////////////////////////// 250 /////////////////////////////////////////////////////////////////////////////////////
251 std::vector<std::vector<G4double>*> 251 std::vector<std::vector<G4double>*>
252 G4VEmAdjointModel::ComputeAdjointCrossSectionV 252 G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
253 G4double kinEnergyScatProj, G4double Z, G4do 253 G4double kinEnergyScatProj, G4double Z, G4double A,
254 G4int nbin_pro_decade) // nb bins pro order 254 G4int nbin_pro_decade) // nb bins pro order of magnitude of energy
255 { 255 {
256 G4Integrator<G4VEmAdjointModel, G4double (G4 256 G4Integrator<G4VEmAdjointModel, G4double (G4VEmAdjointModel::*)(G4double)>
257 integral; 257 integral;
258 fASelectedNucleus = G4lrint(A 258 fASelectedNucleus = G4lrint(A);
259 fZSelectedNucleus = G4lrint(Z 259 fZSelectedNucleus = G4lrint(Z);
260 fKinEnergyScatProjForIntegration = kinEnergy 260 fKinEnergyScatProjForIntegration = kinEnergyScatProj;
261 261
262 // compute the vector of integrated cross se 262 // compute the vector of integrated cross sections
263 G4double minEProj = GetSecondAdjEnergyMinFor 263 G4double minEProj = GetSecondAdjEnergyMinForScatProjToProj(kinEnergyScatProj);
264 G4double maxEProj = GetSecondAdjEnergyMaxFor 264 G4double maxEProj = GetSecondAdjEnergyMaxForScatProjToProj(kinEnergyScatProj);
265 G4double dEmax = maxEProj - kinEnergyScat 265 G4double dEmax = maxEProj - kinEnergyScatProj;
266 G4double dEmin = GetLowEnergyLimit(); 266 G4double dEmin = GetLowEnergyLimit();
267 G4double dE1 = dEmin; 267 G4double dE1 = dEmin;
268 G4double dE2 = dEmin; 268 G4double dE2 = dEmin;
269 269
270 std::vector<G4double>* log_ESec_vector = new 270 std::vector<G4double>* log_ESec_vector = new std::vector<G4double>();
271 std::vector<G4double>* log_Prob_vector = new 271 std::vector<G4double>* log_Prob_vector = new std::vector<G4double>();
272 log_ESec_vector->push_back(std::log(dEmin)); 272 log_ESec_vector->push_back(std::log(dEmin));
273 log_Prob_vector->push_back(-50.); 273 log_Prob_vector->push_back(-50.);
274 G4int nbins = std::max(G4int(std::log10(dEma 274 G4int nbins = std::max(G4int(std::log10(dEmax / dEmin)) * nbin_pro_decade, 5);
275 G4double fE = std::pow(dEmax / dEmin, 1. / n 275 G4double fE = std::pow(dEmax / dEmin, 1. / nbins);
276 276
277 G4double int_cross_section = 0.; 277 G4double int_cross_section = 0.;
278 // Loop checking, 07-Aug-2015, Vladimir Ivan 278 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
279 while(dE1 < dEmax * 0.9999999999999) 279 while(dE1 < dEmax * 0.9999999999999)
280 { 280 {
281 dE2 = dE1 * fE; 281 dE2 = dE1 * fE;
282 int_cross_section += 282 int_cross_section +=
283 integral.Simpson(this, &G4VEmAdjointMode 283 integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction2,
284 minEProj + dE1, std::mi 284 minEProj + dE1, std::min(minEProj + dE2, maxEProj), 5);
285 log_ESec_vector->push_back(std::log(std::m 285 log_ESec_vector->push_back(std::log(std::min(dE2, maxEProj - minEProj)));
286 log_Prob_vector->push_back(std::log(int_cr 286 log_Prob_vector->push_back(std::log(int_cross_section));
287 dE1 = dE2; 287 dE1 = dE2;
288 } 288 }
289 289
290 std::vector<std::vector<G4double>*> res_mat; 290 std::vector<std::vector<G4double>*> res_mat;
291 if(int_cross_section > 0.) 291 if(int_cross_section > 0.)
292 { 292 {
293 res_mat.push_back(log_ESec_vector); 293 res_mat.push_back(log_ESec_vector);
294 res_mat.push_back(log_Prob_vector); 294 res_mat.push_back(log_Prob_vector);
295 } 295 }
296 else { 296 else {
297 delete log_ESec_vector; 297 delete log_ESec_vector;
298 delete log_Prob_vector; 298 delete log_Prob_vector;
299 } 299 }
300 300
301 return res_mat; 301 return res_mat;
302 } 302 }
303 303
304 ////////////////////////////////////////////// 304 ////////////////////////////////////////////////////////////////////////////////
305 std::vector<std::vector<G4double>*> 305 std::vector<std::vector<G4double>*>
306 G4VEmAdjointModel::ComputeAdjointCrossSectionV 306 G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(
307 G4Material* aMaterial, G4double kinEnergyPro 307 G4Material* aMaterial, G4double kinEnergyProd,
308 G4int nbin_pro_decade) // nb bins pro order 308 G4int nbin_pro_decade) // nb bins pro order of magnitude of energy
309 { 309 {
310 G4Integrator<G4VEmAdjointModel, G4double (G4 310 G4Integrator<G4VEmAdjointModel, G4double (G4VEmAdjointModel::*)(G4double)>
311 integral; 311 integral;
312 fSelectedMaterial = aMaterial; 312 fSelectedMaterial = aMaterial;
313 fKinEnergyProdForIntegration = kinEnergyProd 313 fKinEnergyProdForIntegration = kinEnergyProd;
314 314
315 // compute the vector of integrated cross se 315 // compute the vector of integrated cross sections
316 G4double minEProj = GetSecondAdjEnergyMinFor 316 G4double minEProj = GetSecondAdjEnergyMinForProdToProj(kinEnergyProd);
317 G4double maxEProj = GetSecondAdjEnergyMaxFor 317 G4double maxEProj = GetSecondAdjEnergyMaxForProdToProj(kinEnergyProd);
318 G4double E1 = minEProj; 318 G4double E1 = minEProj;
319 std::vector<G4double>* log_ESec_vector = new 319 std::vector<G4double>* log_ESec_vector = new std::vector<G4double>();
320 std::vector<G4double>* log_Prob_vector = new 320 std::vector<G4double>* log_Prob_vector = new std::vector<G4double>();
321 log_ESec_vector->push_back(std::log(E1)); 321 log_ESec_vector->push_back(std::log(E1));
322 log_Prob_vector->push_back(-50.); 322 log_Prob_vector->push_back(-50.);
323 323
324 G4double E2 = 324 G4double E2 =
325 std::pow(10., G4double(G4int(std::log10(mi 325 std::pow(10., G4double(G4int(std::log10(minEProj) * nbin_pro_decade) + 1) /
326 nbin_pro_decade); 326 nbin_pro_decade);
327 G4double fE = std::pow(10., 1. / nbin_pro_de 327 G4double fE = std::pow(10., 1. / nbin_pro_decade);
328 328
329 if(std::pow(fE, 5.) > (maxEProj / minEProj)) 329 if(std::pow(fE, 5.) > (maxEProj / minEProj))
330 fE = std::pow(maxEProj / minEProj, 0.2); 330 fE = std::pow(maxEProj / minEProj, 0.2);
331 331
332 G4double int_cross_section = 0.; 332 G4double int_cross_section = 0.;
333 // Loop checking, 07-Aug-2015, Vladimir Ivan 333 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
334 while(E1 < maxEProj * 0.9999999) 334 while(E1 < maxEProj * 0.9999999)
335 { 335 {
336 int_cross_section += 336 int_cross_section +=
337 integral.Simpson(this, &G4VEmAdjointMode 337 integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction1, E1,
338 std::min(E2, maxEProj * 338 std::min(E2, maxEProj * 0.99999999), 5);
339 log_ESec_vector->push_back(std::log(std::m 339 log_ESec_vector->push_back(std::log(std::min(E2, maxEProj)));
340 log_Prob_vector->push_back(std::log(int_cr 340 log_Prob_vector->push_back(std::log(int_cross_section));
341 E1 = E2; 341 E1 = E2;
342 E2 *= fE; 342 E2 *= fE;
343 } 343 }
344 std::vector<std::vector<G4double>*> res_mat; 344 std::vector<std::vector<G4double>*> res_mat;
345 345
346 if(int_cross_section > 0.) 346 if(int_cross_section > 0.)
347 { 347 {
348 res_mat.push_back(log_ESec_vector); 348 res_mat.push_back(log_ESec_vector);
349 res_mat.push_back(log_Prob_vector); 349 res_mat.push_back(log_Prob_vector);
350 } 350 }
351 else { 351 else {
352 delete log_ESec_vector; 352 delete log_ESec_vector;
353 delete log_Prob_vector; 353 delete log_Prob_vector;
354 } 354 }
355 return res_mat; 355 return res_mat;
356 } 356 }
357 357
358 ////////////////////////////////////////////// 358 ////////////////////////////////////////////////////////////////////////////////
359 std::vector<std::vector<G4double>*> 359 std::vector<std::vector<G4double>*>
360 G4VEmAdjointModel::ComputeAdjointCrossSectionV 360 G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
361 G4Material* aMaterial, G4double kinEnergySca 361 G4Material* aMaterial, G4double kinEnergyScatProj,
362 G4int nbin_pro_decade) // nb bins pro order 362 G4int nbin_pro_decade) // nb bins pro order of magnitude of energy
363 { 363 {
364 G4Integrator<G4VEmAdjointModel, G4double (G4 364 G4Integrator<G4VEmAdjointModel, G4double (G4VEmAdjointModel::*)(G4double)>
365 integral; 365 integral;
366 fSelectedMaterial = aMaterial 366 fSelectedMaterial = aMaterial;
367 fKinEnergyScatProjForIntegration = kinEnergy 367 fKinEnergyScatProjForIntegration = kinEnergyScatProj;
368 368
369 // compute the vector of integrated cross se 369 // compute the vector of integrated cross sections
370 G4double minEProj = GetSecondAdjEnergyMinFor 370 G4double minEProj = GetSecondAdjEnergyMinForScatProjToProj(kinEnergyScatProj);
371 G4double maxEProj = GetSecondAdjEnergyMaxFor 371 G4double maxEProj = GetSecondAdjEnergyMaxForScatProjToProj(kinEnergyScatProj);
372 372
373 G4double dEmax = maxEProj - kinEnergyScatPro 373 G4double dEmax = maxEProj - kinEnergyScatProj;
374 G4double dEmin = GetLowEnergyLimit(); 374 G4double dEmin = GetLowEnergyLimit();
375 G4double dE1 = dEmin; 375 G4double dE1 = dEmin;
376 G4double dE2 = dEmin; 376 G4double dE2 = dEmin;
377 377
378 std::vector<G4double>* log_ESec_vector = new 378 std::vector<G4double>* log_ESec_vector = new std::vector<G4double>();
379 std::vector<G4double>* log_Prob_vector = new 379 std::vector<G4double>* log_Prob_vector = new std::vector<G4double>();
380 log_ESec_vector->push_back(std::log(dEmin)); 380 log_ESec_vector->push_back(std::log(dEmin));
381 log_Prob_vector->push_back(-50.); 381 log_Prob_vector->push_back(-50.);
382 G4int nbins = std::max(int(std::log10(dEmax 382 G4int nbins = std::max(int(std::log10(dEmax / dEmin)) * nbin_pro_decade, 5);
383 G4double fE = std::pow(dEmax / dEmin, 1. / n 383 G4double fE = std::pow(dEmax / dEmin, 1. / nbins);
384 384
385 G4double int_cross_section = 0.; 385 G4double int_cross_section = 0.;
386 // Loop checking, 07-Aug-2015, Vladimir Ivan 386 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
387 while(dE1 < dEmax * 0.9999999999999) 387 while(dE1 < dEmax * 0.9999999999999)
388 { 388 {
389 dE2 = dE1 * fE; 389 dE2 = dE1 * fE;
390 int_cross_section += 390 int_cross_section +=
391 integral.Simpson(this, &G4VEmAdjointMode 391 integral.Simpson(this, &G4VEmAdjointModel::DiffCrossSectionFunction2,
392 minEProj + dE1, std::mi 392 minEProj + dE1, std::min(minEProj + dE2, maxEProj), 5);
393 log_ESec_vector->push_back(std::log(std::m 393 log_ESec_vector->push_back(std::log(std::min(dE2, maxEProj - minEProj)));
394 log_Prob_vector->push_back(std::log(int_cr 394 log_Prob_vector->push_back(std::log(int_cross_section));
395 dE1 = dE2; 395 dE1 = dE2;
396 } 396 }
397 397
398 std::vector<std::vector<G4double>*> res_mat; 398 std::vector<std::vector<G4double>*> res_mat;
399 if(int_cross_section > 0.) 399 if(int_cross_section > 0.)
400 { 400 {
401 res_mat.push_back(log_ESec_vector); 401 res_mat.push_back(log_ESec_vector);
402 res_mat.push_back(log_Prob_vector); 402 res_mat.push_back(log_Prob_vector);
403 } 403 }
404 else { 404 else {
405 delete log_ESec_vector; 405 delete log_ESec_vector;
406 delete log_Prob_vector; 406 delete log_Prob_vector;
407 } 407 }
408 408
409 return res_mat; 409 return res_mat;
410 } 410 }
411 411
412 ////////////////////////////////////////////// 412 //////////////////////////////////////////////////////////////////////////////
413 G4double G4VEmAdjointModel::SampleAdjSecEnergy 413 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(
414 std::size_t MatrixIndex, G4double aPrimEnerg 414 std::size_t MatrixIndex, G4double aPrimEnergy, G4bool isScatProjToProj)
415 { 415 {
416 G4AdjointCSMatrix* theMatrix = (*fCSMatrixPr 416 G4AdjointCSMatrix* theMatrix = (*fCSMatrixProdToProjBackScat)[MatrixIndex];
417 if(isScatProjToProj) 417 if(isScatProjToProj)
418 theMatrix = (*fCSMatrixProjToProjBackScat) 418 theMatrix = (*fCSMatrixProjToProjBackScat)[MatrixIndex];
419 std::vector<G4double>* theLogPrimEnergyVecto 419 std::vector<G4double>* theLogPrimEnergyVector =
420 theMatrix->GetLogPrimEnergyVector(); 420 theMatrix->GetLogPrimEnergyVector();
421 421
422 if(theLogPrimEnergyVector->empty()) 422 if(theLogPrimEnergyVector->empty())
423 { 423 {
424 G4cout << "No data are contained in the gi 424 G4cout << "No data are contained in the given AdjointCSMatrix!" << G4endl;
425 G4cout << "The sampling procedure will be 425 G4cout << "The sampling procedure will be stopped." << G4endl;
426 return 0.; 426 return 0.;
427 } 427 }
428 428
429 G4AdjointInterpolator* theInterpolator = G4A 429 G4AdjointInterpolator* theInterpolator = G4AdjointInterpolator::GetInstance();
430 G4double aLogPrimEnergy = std 430 G4double aLogPrimEnergy = std::log(aPrimEnergy);
431 G4int ind = (G4int)theInterpolator->FindPosi 431 G4int ind = (G4int)theInterpolator->FindPositionForLogVector(
432 aLogPrimEnergy, *theLogPrimEnergyVector); 432 aLogPrimEnergy, *theLogPrimEnergyVector);
433 433
434 G4double aLogPrimEnergy1, aLogPrimEnergy2; 434 G4double aLogPrimEnergy1, aLogPrimEnergy2;
435 G4double aLogCS1, aLogCS2; 435 G4double aLogCS1, aLogCS2;
436 G4double log01, log02; 436 G4double log01, log02;
437 std::vector<G4double>* aLogSecondEnergyVecto 437 std::vector<G4double>* aLogSecondEnergyVector1 = nullptr;
438 std::vector<G4double>* aLogSecondEnergyVecto 438 std::vector<G4double>* aLogSecondEnergyVector2 = nullptr;
439 std::vector<G4double>* aLogProbVector1 439 std::vector<G4double>* aLogProbVector1 = nullptr;
440 std::vector<G4double>* aLogProbVector2 440 std::vector<G4double>* aLogProbVector2 = nullptr;
441 std::vector<std::size_t>* aLogProbVectorInde 441 std::vector<std::size_t>* aLogProbVectorIndex1 = nullptr;
442 std::vector<std::size_t>* aLogProbVectorInde 442 std::vector<std::size_t>* aLogProbVectorIndex2 = nullptr;
443 443
444 theMatrix->GetData(ind, aLogPrimEnergy1, aLo 444 theMatrix->GetData(ind, aLogPrimEnergy1, aLogCS1, log01,
445 aLogSecondEnergyVector1, 445 aLogSecondEnergyVector1, aLogProbVector1,
446 aLogProbVectorIndex1 ); 446 aLogProbVectorIndex1 );
447 theMatrix->GetData(ind + 1, aLogPrimEnergy2, 447 theMatrix->GetData(ind + 1, aLogPrimEnergy2, aLogCS2, log02,
448 aLogSecondEnergyVector2, 448 aLogSecondEnergyVector2, aLogProbVector2,
449 aLogProbVectorIndex2); 449 aLogProbVectorIndex2);
450 450
451 if (! (aLogProbVector1 && aLogProbVector2 && 451 if (! (aLogProbVector1 && aLogProbVector2 &&
452 aLogSecondEnergyVector1 && aLogSeco 452 aLogSecondEnergyVector1 && aLogSecondEnergyVector2)){
453 return 0.; 453 return 0.;
454 } 454 }
455 455
456 G4double rand_var = G4UniformRand(); 456 G4double rand_var = G4UniformRand();
457 G4double log_rand_var = std::log(rand_var); 457 G4double log_rand_var = std::log(rand_var);
458 G4double log_Tcut = std::log(fTcutSecond 458 G4double log_Tcut = std::log(fTcutSecond);
459 G4double Esec = 0.; 459 G4double Esec = 0.;
460 G4double log_dE1, log_dE2; 460 G4double log_dE1, log_dE2;
461 G4double log_rand_var1, log_rand_var2; 461 G4double log_rand_var1, log_rand_var2;
462 G4double log_E1, log_E2; 462 G4double log_E1, log_E2;
463 log_rand_var1 = log_rand_var; 463 log_rand_var1 = log_rand_var;
464 log_rand_var2 = log_rand_var; 464 log_rand_var2 = log_rand_var;
465 465
466 G4double Emin = 0.; 466 G4double Emin = 0.;
467 G4double Emax = 0.; 467 G4double Emax = 0.;
468 if(theMatrix->IsScatProjToProj()) 468 if(theMatrix->IsScatProjToProj())
469 { // case where Tcut plays a role 469 { // case where Tcut plays a role
470 Emin = GetSecondAdjEnergyMinForScatProjToP 470 Emin = GetSecondAdjEnergyMinForScatProjToProj(aPrimEnergy, fTcutSecond);
471 Emax = GetSecondAdjEnergyMaxForScatProjToP 471 Emax = GetSecondAdjEnergyMaxForScatProjToProj(aPrimEnergy);
472 G4double dE = 0.; 472 G4double dE = 0.;
473 if(Emin < Emax) 473 if(Emin < Emax)
474 { 474 {
475 if(fApplyCutInRange) 475 if(fApplyCutInRange)
476 { 476 {
477 if(fSecondPartSameType && fTcutSecond 477 if(fSecondPartSameType && fTcutSecond > aPrimEnergy)
478 return aPrimEnergy; 478 return aPrimEnergy;
479 479
480 log_rand_var1 = log_rand_var + 480 log_rand_var1 = log_rand_var +
481 theInterpolator->Inter 481 theInterpolator->InterpolateForLogVector(
482 log_Tcut, *aLogSecon 482 log_Tcut, *aLogSecondEnergyVector1, *aLogProbVector1);
483 log_rand_var2 = log_rand_var + 483 log_rand_var2 = log_rand_var +
484 theInterpolator->Inter 484 theInterpolator->InterpolateForLogVector(
485 log_Tcut, *aLogSecon 485 log_Tcut, *aLogSecondEnergyVector2, *aLogProbVector2);
486 } 486 }
487 log_dE1 = theInterpolator->Interpolate(l 487 log_dE1 = theInterpolator->Interpolate(log_rand_var1, *aLogProbVector1,
488 * 488 *aLogSecondEnergyVector1, "Lin");
489 log_dE2 = theInterpolator->Interpolate(l 489 log_dE2 = theInterpolator->Interpolate(log_rand_var2, *aLogProbVector2,
490 * 490 *aLogSecondEnergyVector2, "Lin");
491 dE = std::exp(theInterpolator->Line 491 dE = std::exp(theInterpolator->LinearInterpolation(
492 aLogPrimEnergy, aLogPrimEnergy1, aLogP 492 aLogPrimEnergy, aLogPrimEnergy1, aLogPrimEnergy2, log_dE1, log_dE2));
493 } 493 }
494 494
495 Esec = aPrimEnergy + dE; 495 Esec = aPrimEnergy + dE;
496 Esec = std::max(Esec, Emin); 496 Esec = std::max(Esec, Emin);
497 Esec = std::min(Esec, Emax); 497 Esec = std::min(Esec, Emax);
498 } 498 }
499 else 499 else
500 { // Tcut condition is already full-filled 500 { // Tcut condition is already full-filled
501 501
502 log_E1 = theInterpolator->Interpolate(log_ 502 log_E1 = theInterpolator->Interpolate(log_rand_var, *aLogProbVector1,
503 *aLo 503 *aLogSecondEnergyVector1, "Lin");
504 log_E2 = theInterpolator->Interpolate(log_ 504 log_E2 = theInterpolator->Interpolate(log_rand_var, *aLogProbVector2,
505 *aLo 505 *aLogSecondEnergyVector2, "Lin");
506 506
507 Esec = std::exp(theInterpolator->LinearInt 507 Esec = std::exp(theInterpolator->LinearInterpolation(
508 aLogPrimEnergy, aLogPrimEnergy1, aLogPri 508 aLogPrimEnergy, aLogPrimEnergy1, aLogPrimEnergy2, log_E1, log_E2));
509 Emin = GetSecondAdjEnergyMinForProdToProj( 509 Emin = GetSecondAdjEnergyMinForProdToProj(aPrimEnergy);
510 Emax = GetSecondAdjEnergyMaxForProdToProj( 510 Emax = GetSecondAdjEnergyMaxForProdToProj(aPrimEnergy);
511 Esec = std::max(Esec, Emin); 511 Esec = std::max(Esec, Emin);
512 Esec = std::min(Esec, Emax); 512 Esec = std::min(Esec, Emax);
513 } 513 }
514 return Esec; 514 return Esec;
515 } 515 }
516 516
517 ////////////////////////////////////////////// 517 //////////////////////////////////////////////////////////////////////////////
518 G4double G4VEmAdjointModel::SampleAdjSecEnergy 518 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(
519 G4double aPrimEnergy, G4bool isScatProjToPro 519 G4double aPrimEnergy, G4bool isScatProjToProj)
520 { 520 {
521 SelectCSMatrix(isScatProjToProj); 521 SelectCSMatrix(isScatProjToProj);
522 return SampleAdjSecEnergyFromCSMatrix(fCSMat 522 return SampleAdjSecEnergyFromCSMatrix(fCSMatrixUsed, aPrimEnergy,
523 isScat 523 isScatProjToProj);
524 } 524 }
525 525
526 ////////////////////////////////////////////// 526 //////////////////////////////////////////////////////////////////////////////
527 void G4VEmAdjointModel::SelectCSMatrix(G4bool 527 void G4VEmAdjointModel::SelectCSMatrix(G4bool isScatProjToProj)
528 { 528 {
529 fCSMatrixUsed = 0; 529 fCSMatrixUsed = 0;
530 if(!fUseMatrixPerElement) 530 if(!fUseMatrixPerElement)
531 fCSMatrixUsed = fCurrentMaterial->GetIndex 531 fCSMatrixUsed = fCurrentMaterial->GetIndex();
532 else if(!fOneMatrixForAllElements) 532 else if(!fOneMatrixForAllElements)
533 { // Select Material 533 { // Select Material
534 std::vector<G4double>* CS_Vs_Element = &fE 534 std::vector<G4double>* CS_Vs_Element = &fElementCSScatProjToProj;
535 fLastCS = fLa 535 fLastCS = fLastAdjointCSForScatProjToProj;
536 if(!isScatProjToProj) 536 if(!isScatProjToProj)
537 { 537 {
538 CS_Vs_Element = &fElementCSProdToProj; 538 CS_Vs_Element = &fElementCSProdToProj;
539 fLastCS = fLastAdjointCSForProdToP 539 fLastCS = fLastAdjointCSForProdToProj;
540 } 540 }
541 G4double SumCS = 0.; 541 G4double SumCS = 0.;
542 std::size_t ind = 0; 542 std::size_t ind = 0;
543 for(std::size_t i = 0; i < CS_Vs_Element-> 543 for(std::size_t i = 0; i < CS_Vs_Element->size(); ++i)
544 { 544 {
545 SumCS += (*CS_Vs_Element)[i]; 545 SumCS += (*CS_Vs_Element)[i];
546 if(G4UniformRand() <= SumCS / fLastCS) 546 if(G4UniformRand() <= SumCS / fLastCS)
547 { 547 {
548 ind = i; 548 ind = i;
549 break; 549 break;
550 } 550 }
551 } 551 }
552 fCSMatrixUsed = fCurrentMaterial->GetEleme 552 fCSMatrixUsed = fCurrentMaterial->GetElement((G4int)ind)->GetIndex();
553 } 553 }
554 } 554 }
555 555
556 ////////////////////////////////////////////// 556 //////////////////////////////////////////////////////////////////////////////
557 G4double G4VEmAdjointModel::SampleAdjSecEnergy 557 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(
558 G4double prim_energy, G4bool isScatProjToPro 558 G4double prim_energy, G4bool isScatProjToProj)
559 { 559 {
560 // here we try to use the rejection method 560 // here we try to use the rejection method
561 constexpr G4int iimax = 1000; 561 constexpr G4int iimax = 1000;
562 G4double E = 0.; 562 G4double E = 0.;
563 G4double x, xmin, greject; 563 G4double x, xmin, greject;
564 if(isScatProjToProj) 564 if(isScatProjToProj)
565 { 565 {
566 G4double Emax = GetSecondAdjEnergyMaxForSc 566 G4double Emax = GetSecondAdjEnergyMaxForScatProjToProj(prim_energy);
567 G4double Emin = prim_energy + fTcutSecond; 567 G4double Emin = prim_energy + fTcutSecond;
568 xmin = Emin / Emax; 568 xmin = Emin / Emax;
569 G4double grejmax = 569 G4double grejmax =
570 DiffCrossSectionPerAtomPrimToScatPrim(Em 570 DiffCrossSectionPerAtomPrimToScatPrim(Emin, prim_energy, 1) * prim_energy;
571 571
572 G4int ii = 0; 572 G4int ii = 0;
573 do 573 do
574 { 574 {
575 // q = G4UniformRand(); 575 // q = G4UniformRand();
576 x = 1. / (G4UniformRand() * (1. / xmin - 576 x = 1. / (G4UniformRand() * (1. / xmin - 1.) + 1.);
577 E = x * Emax; 577 E = x * Emax;
578 greject = 578 greject =
579 DiffCrossSectionPerAtomPrimToScatPrim( 579 DiffCrossSectionPerAtomPrimToScatPrim(E, prim_energy, 1) * prim_energy;
580 ++ii; 580 ++ii;
581 if(ii >= iimax) 581 if(ii >= iimax)
582 { 582 {
583 break; 583 break;
584 } 584 }
585 } 585 }
586 // Loop checking, 07-Aug-2015, Vladimir Iv 586 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
587 while(greject < G4UniformRand() * grejmax) 587 while(greject < G4UniformRand() * grejmax);
588 } 588 }
589 else 589 else
590 { 590 {
591 G4double Emax = GetSecondAdjEnergyMaxForPr 591 G4double Emax = GetSecondAdjEnergyMaxForProdToProj(prim_energy);
592 G4double Emin = GetSecondAdjEnergyMinForPr 592 G4double Emin = GetSecondAdjEnergyMinForProdToProj(prim_energy);
593 xmin = Emin / Emax; 593 xmin = Emin / Emax;
594 G4double grejmax = 594 G4double grejmax =
595 DiffCrossSectionPerAtomPrimToSecond(Emin 595 DiffCrossSectionPerAtomPrimToSecond(Emin, prim_energy, 1);
596 G4int ii = 0; 596 G4int ii = 0;
597 do 597 do
598 { 598 {
599 x = std::pow(xmin, G4UniformRand() 599 x = std::pow(xmin, G4UniformRand());
600 E = x * Emax; 600 E = x * Emax;
601 greject = DiffCrossSectionPerAtomPrimToS 601 greject = DiffCrossSectionPerAtomPrimToSecond(E, prim_energy, 1);
602 ++ii; 602 ++ii;
603 if(ii >= iimax) 603 if(ii >= iimax)
604 { 604 {
605 break; 605 break;
606 } 606 }
607 } 607 }
608 // Loop checking, 07-Aug-2015, Vladimir Iv 608 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
609 while(greject < G4UniformRand() * grejmax) 609 while(greject < G4UniformRand() * grejmax);
610 } 610 }
611 611
612 return E; 612 return E;
613 } 613 }
614 614
615 ////////////////////////////////////////////// 615 ////////////////////////////////////////////////////////////////////////////////
616 void G4VEmAdjointModel::CorrectPostStepWeight( 616 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
617 617 G4double old_weight,
618 618 G4double adjointPrimKinEnergy,
619 619 G4double projectileKinEnergy,
620 620 G4bool isScatProjToProj)
621 { 621 {
622 G4double new_weight = old_weight; 622 G4double new_weight = old_weight;
623 G4double w_corr = 623 G4double w_corr =
624 fCSManager->GetPostStepWeightCorrection() 624 fCSManager->GetPostStepWeightCorrection() / fCsBiasingFactor;
625 625
626 fLastCS = fLastAdjointCSForScatProjToProj; 626 fLastCS = fLastAdjointCSForScatProjToProj;
627 if(!isScatProjToProj) 627 if(!isScatProjToProj)
628 fLastCS = fLastAdjointCSForProdToProj; 628 fLastCS = fLastAdjointCSForProdToProj;
629 if((adjointPrimKinEnergy - fPreStepEnergy) / 629 if((adjointPrimKinEnergy - fPreStepEnergy) / fPreStepEnergy > 0.001)
630 { 630 {
631 G4double post_stepCS = AdjointCrossSection 631 G4double post_stepCS = AdjointCrossSection(
632 fCurrentCouple, adjointPrimKinEnergy, is 632 fCurrentCouple, adjointPrimKinEnergy, isScatProjToProj);
633 if(post_stepCS > 0. && fLastCS > 0.) 633 if(post_stepCS > 0. && fLastCS > 0.)
634 w_corr *= post_stepCS / fLastCS; 634 w_corr *= post_stepCS / fLastCS;
635 } 635 }
636 636
637 new_weight *= w_corr; 637 new_weight *= w_corr;
638 new_weight *= projectileKinEnergy / adjointP 638 new_weight *= projectileKinEnergy / adjointPrimKinEnergy;
639 // This is needed due to the biasing of diff 639 // This is needed due to the biasing of diff CS
640 // by the factor adjointPrimKinEnergy/projec 640 // by the factor adjointPrimKinEnergy/projectileKinEnergy
641 641
642 fParticleChange->SetParentWeightByProcess(fa 642 fParticleChange->SetParentWeightByProcess(false);
643 fParticleChange->SetSecondaryWeightByProcess 643 fParticleChange->SetSecondaryWeightByProcess(false);
644 fParticleChange->ProposeParentWeight(new_wei 644 fParticleChange->ProposeParentWeight(new_weight);
645 } 645 }
646 646
647 ////////////////////////////////////////////// 647 //////////////////////////////////////////////////////////////////////////////
648 G4double G4VEmAdjointModel::GetSecondAdjEnergy 648 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProj(
649 G4double kinEnergyScatProj) 649 G4double kinEnergyScatProj)
650 { 650 {
651 G4double maxEProj = GetHighEnergyLimit(); 651 G4double maxEProj = GetHighEnergyLimit();
652 if(fSecondPartSameType) 652 if(fSecondPartSameType)
653 maxEProj = std::min(kinEnergyScatProj * 2. 653 maxEProj = std::min(kinEnergyScatProj * 2., maxEProj);
654 return maxEProj; 654 return maxEProj;
655 } 655 }
656 656
657 ////////////////////////////////////////////// 657 //////////////////////////////////////////////////////////////////////////////
658 G4double G4VEmAdjointModel::GetSecondAdjEnergy 658 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProj(
659 G4double primAdjEnergy, G4double tcut) 659 G4double primAdjEnergy, G4double tcut)
660 { 660 {
661 G4double Emin = primAdjEnergy; 661 G4double Emin = primAdjEnergy;
662 if(fApplyCutInRange) 662 if(fApplyCutInRange)
663 Emin += tcut; 663 Emin += tcut;
664 return Emin; 664 return Emin;
665 } 665 }
666 666
667 ////////////////////////////////////////////// 667 //////////////////////////////////////////////////////////////////////////////
668 G4double G4VEmAdjointModel::GetSecondAdjEnergy 668 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProj(G4double)
669 { 669 {
670 return fHighEnergyLimit; 670 return fHighEnergyLimit;
671 } 671 }
672 672
673 ////////////////////////////////////////////// 673 //////////////////////////////////////////////////////////////////////////////
674 G4double G4VEmAdjointModel::GetSecondAdjEnergy 674 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProj(
675 G4double primAdjEnergy) 675 G4double primAdjEnergy)
676 { 676 {
677 G4double minEProj = primAdjEnergy; 677 G4double minEProj = primAdjEnergy;
678 if(fSecondPartSameType) 678 if(fSecondPartSameType)
679 minEProj = primAdjEnergy * 2.; 679 minEProj = primAdjEnergy * 2.;
680 return minEProj; 680 return minEProj;
681 } 681 }
682 682
683 ////////////////////////////////////////////// 683 ////////////////////////////////////////////////////////////////////////////////////////////
684 void G4VEmAdjointModel::DefineCurrentMaterial( 684 void G4VEmAdjointModel::DefineCurrentMaterial(
685 const G4MaterialCutsCouple* couple) 685 const G4MaterialCutsCouple* couple)
686 { 686 {
687 if(couple != fCurrentCouple) 687 if(couple != fCurrentCouple)
688 { 688 {
689 fCurrentCouple = const_cast<G4MaterialCu 689 fCurrentCouple = const_cast<G4MaterialCutsCouple*>(couple);
690 fCurrentMaterial = const_cast<G4Material*> 690 fCurrentMaterial = const_cast<G4Material*>(couple->GetMaterial());
691 std::size_t idx = 56; 691 std::size_t idx = 56;
692 fTcutSecond = 1.e-11; 692 fTcutSecond = 1.e-11;
693 if(fAdjEquivDirectSecondPart) 693 if(fAdjEquivDirectSecondPart)
694 { 694 {
695 if(fAdjEquivDirectSecondPart == G4Adjoin 695 if(fAdjEquivDirectSecondPart == G4AdjointGamma::AdjointGamma())
696 idx = 0; 696 idx = 0;
697 else if(fAdjEquivDirectSecondPart == G4A 697 else if(fAdjEquivDirectSecondPart == G4AdjointElectron::AdjointElectron())
698 idx = 1; 698 idx = 1;
699 else if(fAdjEquivDirectSecondPart == G4A 699 else if(fAdjEquivDirectSecondPart == G4AdjointPositron::AdjointPositron())
700 idx = 2; 700 idx = 2;
701 if(idx < 56) 701 if(idx < 56)
702 { 702 {
703 const std::vector<G4double>* aVec = 703 const std::vector<G4double>* aVec =
704 G4ProductionCutsTable::GetProduction 704 G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(
705 idx); 705 idx);
706 fTcutSecond = (*aVec)[couple->GetIndex 706 fTcutSecond = (*aVec)[couple->GetIndex()];
707 } 707 }
708 } 708 }
709 } 709 }
710 } 710 }
711 711
712 ////////////////////////////////////////////// 712 ////////////////////////////////////////////////////////////////////////////////////////////
713 void G4VEmAdjointModel::SetHighEnergyLimit(G4d 713 void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
714 { 714 {
715 fHighEnergyLimit = aVal; 715 fHighEnergyLimit = aVal;
716 if(fDirectModel) 716 if(fDirectModel)
717 fDirectModel->SetHighEnergyLimit(aVal); 717 fDirectModel->SetHighEnergyLimit(aVal);
718 } 718 }
719 719
720 ////////////////////////////////////////////// 720 ////////////////////////////////////////////////////////////////////////////////////////////
721 void G4VEmAdjointModel::SetLowEnergyLimit(G4do 721 void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
722 { 722 {
723 fLowEnergyLimit = aVal; 723 fLowEnergyLimit = aVal;
724 if(fDirectModel) 724 if(fDirectModel)
725 fDirectModel->SetLowEnergyLimit(aVal); 725 fDirectModel->SetLowEnergyLimit(aVal);
726 } 726 }
727 727
728 ////////////////////////////////////////////// 728 ////////////////////////////////////////////////////////////////////////////////////////////
729 void G4VEmAdjointModel::SetAdjointEquivalentOf 729 void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(
730 G4ParticleDefinition* aPart) 730 G4ParticleDefinition* aPart)
731 { 731 {
732 fAdjEquivDirectPrimPart = aPart; 732 fAdjEquivDirectPrimPart = aPart;
733 if(fAdjEquivDirectPrimPart->GetParticleName( 733 if(fAdjEquivDirectPrimPart->GetParticleName() == "adj_e-")
734 fDirectPrimaryPart = G4Electron::Electron( 734 fDirectPrimaryPart = G4Electron::Electron();
735 else if(fAdjEquivDirectPrimPart->GetParticle 735 else if(fAdjEquivDirectPrimPart->GetParticleName() == "adj_gamma")
736 fDirectPrimaryPart = G4Gamma::Gamma(); 736 fDirectPrimaryPart = G4Gamma::Gamma();
737 } 737 }
738 738