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