Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4VEmAdjointModel.cc
Go to the documentation of this file.
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////////////////////////////////////////////////////////////////////////////////
42G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam)
43 : fName(nam)
44{
45 fCSManager = G4AdjointCSManager::GetAdjointCSManager();
46 fCSManager->RegisterEmAdjointModel(this);
47}
48
49///////////////////////////////////e////////////////////////////////////////////
50G4VEmAdjointModel::~G4VEmAdjointModel()
51{
52 delete fCSMatrixProdToProjBackScat;
53 fCSMatrixProdToProjBackScat = nullptr;
54
55 delete fCSMatrixProjToProjBackScat;
56 fCSMatrixProjToProjBackScat = nullptr;
57}
58
59////////////////////////////////////////////////////////////////////////////////
60G4double 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////////////////////////////////////////////////////////////////////////////////
82G4double 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////////////////////////////////////////////////////////////////////////////////
105G4double 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////////////////////////////////////////////////////////////////////////////////
119G4double 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////////////////////////////////////////////////////////////////////////////////
140G4double 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////////////////////////////////////////////////////////////////////////////////
155G4double 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////////////////////////////////////////////////////////////////////////////////
176G4double 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////////////////////////////////////////////////////////////////////////////////
197std::vector<std::vector<G4double>*>
198G4VEmAdjointModel::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/////////////////////////////////////////////////////////////////////////////////////
251std::vector<std::vector<G4double>*>
252G4VEmAdjointModel::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////////////////////////////////////////////////////////////////////////////////
305std::vector<std::vector<G4double>*>
306G4VEmAdjointModel::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////////////////////////////////////////////////////////////////////////////////
359std::vector<std::vector<G4double>*>
360G4VEmAdjointModel::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//////////////////////////////////////////////////////////////////////////////
413G4double 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//////////////////////////////////////////////////////////////////////////////
518G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(
519 G4double aPrimEnergy, G4bool isScatProjToProj)
520{
521 SelectCSMatrix(isScatProjToProj);
522 return SampleAdjSecEnergyFromCSMatrix(fCSMatrixUsed, aPrimEnergy,
523 isScatProjToProj);
524}
525
526//////////////////////////////////////////////////////////////////////////////
527void 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//////////////////////////////////////////////////////////////////////////////
557G4double 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////////////////////////////////////////////////////////////////////////////////
616void 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//////////////////////////////////////////////////////////////////////////////
648G4double 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//////////////////////////////////////////////////////////////////////////////
658G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProj(
659 G4double primAdjEnergy, G4double tcut)
660{
661 G4double Emin = primAdjEnergy;
662 if(fApplyCutInRange)
663 Emin += tcut;
664 return Emin;
665}
666
667//////////////////////////////////////////////////////////////////////////////
668G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProj(G4double)
669{
670 return fHighEnergyLimit;
671}
672
673//////////////////////////////////////////////////////////////////////////////
674G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProj(
675 G4double primAdjEnergy)
676{
677 G4double minEProj = primAdjEnergy;
678 if(fSecondPartSameType)
679 minEProj = primAdjEnergy * 2.;
680 return minEProj;
681}
682
683////////////////////////////////////////////////////////////////////////////////////////////
684void 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////////////////////////////////////////////////////////////////////////////////////////////
713void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
714{
715 fHighEnergyLimit = aVal;
716 if(fDirectModel)
717 fDirectModel->SetHighEnergyLimit(aVal);
718}
719
720////////////////////////////////////////////////////////////////////////////////////////////
721void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
722{
723 fLowEnergyLimit = aVal;
724 if(fDirectModel)
725 fDirectModel->SetLowEnergyLimit(aVal);
726}
727
728////////////////////////////////////////////////////////////////////////////////////////////
729void 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}
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
Z
const G4int Z[17]
Definition: G4WaterStopping.cc:51
A
const G4double A[17]
Definition: G4WaterStopping.cc:53
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:52
G4AdjointCSManager::ComputeAdjointCS
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool isScatProjToProj, std::vector< G4double > &AdjointCS_for_each_element)
Definition: G4AdjointCSManager.cc:584
G4AdjointCSManager::GetPostStepWeightCorrection
G4double GetPostStepWeightCorrection()
Definition: G4AdjointCSManager.cc:578
G4AdjointCSManager::GetAdjointCSManager
static G4AdjointCSManager * GetAdjointCSManager()
Definition: G4AdjointCSManager.cc:57
G4AdjointCSManager::RegisterEmAdjointModel
std::size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
Definition: G4AdjointCSManager.cc:150
G4AdjointCSMatrix::GetData
G4bool GetData(unsigned int i, G4double &aPrimEnergy, G4double &aCS, G4double &log0, std::vector< double > *&aLogSecondEnergyVector, std::vector< double > *&aLogProbVector, std::vector< size_t > *&aLogProbVectorIndex)
Definition: G4AdjointCSMatrix.cc:122
G4AdjointCSMatrix::GetLogPrimEnergyVector
std::vector< double > * GetLogPrimEnergyVector()
Definition: G4AdjointCSMatrix.hh:66
G4AdjointElectron::AdjointElectron
static G4AdjointElectron * AdjointElectron()
Definition: G4AdjointElectron.cc:99
G4AdjointGamma::AdjointGamma
static G4AdjointGamma * AdjointGamma()
Definition: G4AdjointGamma.cc:83
G4AdjointInterpolator::LinearInterpolation
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
Definition: G4AdjointInterpolator.cc:54
G4AdjointInterpolator::Interpolate
G4double Interpolate(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec, G4String InterPolMethod="Log")
Definition: G4AdjointInterpolator.cc:163
G4AdjointInterpolator::GetInstance
static G4AdjointInterpolator * GetInstance()
Definition: G4AdjointInterpolator.cc:38
G4AdjointInterpolator::FindPositionForLogVector
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
Definition: G4AdjointInterpolator.cc:155
G4AdjointInterpolator::InterpolateForLogVector
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
Definition: G4AdjointInterpolator.cc:198
G4AdjointPositron::AdjointPositron
static G4AdjointPositron * AdjointPositron()
Definition: G4AdjointPositron.cc:98
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4Element::GetIndex
size_t GetIndex() const
Definition: G4Element.hh:182
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
G4MaterialCutsCouple::GetMaterial
const G4Material * GetMaterial() const
Definition: G4MaterialCutsCouple.hh:166
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:197
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:255
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:109
G4ProductionCutsTable::GetEnergyCutsVector
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
Definition: G4ProductionCutsTable.hh:233
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:58
G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProj
virtual G4double GetSecondAdjEnergyMaxForProdToProj(G4double primAdjEnergy)
Definition: G4VEmAdjointModel.cc:668
G4VEmAdjointModel::fLastAdjointCSForScatProjToProj
G4double fLastAdjointCSForScatProjToProj
Definition: G4VEmAdjointModel.hh:274
G4VEmAdjointModel::fCSMatrixProdToProjBackScat
std::vector< G4AdjointCSMatrix * > * fCSMatrixProdToProjBackScat
Definition: G4VEmAdjointModel.hh:264
G4VEmAdjointModel::fAdjEquivDirectSecondPart
G4ParticleDefinition * fAdjEquivDirectSecondPart
Definition: G4VEmAdjointModel.hh:260
G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
Definition: G4VEmAdjointModel.cc:140
G4VEmAdjointModel::~G4VEmAdjointModel
virtual ~G4VEmAdjointModel()
Definition: G4VEmAdjointModel.cc:50
G4VEmAdjointModel::fKinEnergyScatProjForIntegration
G4double fKinEnergyScatProjForIntegration
Definition: G4VEmAdjointModel.hh:271
G4VEmAdjointModel::DiffCrossSectionFunction1
G4double DiffCrossSectionFunction1(G4double kinEnergyProj)
Definition: G4VEmAdjointModel.cc:155
G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProj
virtual G4double GetSecondAdjEnergyMinForScatProjToProj(G4double primAdjEnergy, G4double tcut=0.)
Definition: G4VEmAdjointModel.cc:658
G4VEmAdjointModel::fKinEnergyProdForIntegration
G4double fKinEnergyProdForIntegration
Definition: G4VEmAdjointModel.hh:270
G4VEmAdjointModel::fAdjEquivDirectPrimPart
G4ParticleDefinition * fAdjEquivDirectPrimPart
Definition: G4VEmAdjointModel.hh:259
G4VEmAdjointModel::fCurrentCouple
G4MaterialCutsCouple * fCurrentCouple
Definition: G4VEmAdjointModel.hh:256
G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProj
virtual G4double GetSecondAdjEnergyMaxForScatProjToProj(G4double primAdjEnergy)
Definition: G4VEmAdjointModel.cc:648
G4VEmAdjointModel::fCurrentMaterial
G4Material * fCurrentMaterial
Definition: G4VEmAdjointModel.hh:255
G4VEmAdjointModel::fCSMatrixProjToProjBackScat
std::vector< G4AdjointCSMatrix * > * fCSMatrixProjToProjBackScat
Definition: G4VEmAdjointModel.hh:265
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:198
G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
Definition: G4VEmAdjointModel.cc:119
G4VEmAdjointModel::CorrectPostStepWeight
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool isScatProjToProj)
Definition: G4VEmAdjointModel.cc:616
G4VEmAdjointModel::fElementCSScatProjToProj
std::vector< G4double > fElementCSScatProjToProj
Definition: G4VEmAdjointModel.hh:267
G4VEmAdjointModel::GetLowEnergyLimit
G4double GetLowEnergyLimit()
Definition: G4VEmAdjointModel.hh:154
G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
Definition: G4VEmAdjointModel.cc:105
G4VEmAdjointModel::fCSManager
G4AdjointCSManager * fCSManager
Definition: G4VEmAdjointModel.hh:249
G4VEmAdjointModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double aVal)
Definition: G4VEmAdjointModel.cc:721
G4VEmAdjointModel::DefineCurrentMaterial
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
Definition: G4VEmAdjointModel.cc:684
G4VEmAdjointModel::SelectCSMatrix
void SelectCSMatrix(G4bool isScatProjToProj)
Definition: G4VEmAdjointModel.cc:527
G4VEmAdjointModel::AdjointCrossSection
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool isScatProjToProj)
Definition: G4VEmAdjointModel.cc:60
G4VEmAdjointModel::fDirectModel
G4VEmModel * fDirectModel
Definition: G4VEmAdjointModel.hh:250
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:306
G4VEmAdjointModel::fSelectedMaterial
G4Material * fSelectedMaterial
Definition: G4VEmAdjointModel.hh:254
G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProj
virtual G4double GetSecondAdjEnergyMinForProdToProj(G4double primAdjEnergy)
Definition: G4VEmAdjointModel.cc:674
G4VEmAdjointModel::DiffCrossSectionFunction2
G4double DiffCrossSectionFunction2(G4double kinEnergyProj)
Definition: G4VEmAdjointModel.cc:176
G4VEmAdjointModel::fDirectPrimaryPart
G4ParticleDefinition * fDirectPrimaryPart
Definition: G4VEmAdjointModel.hh:261
G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition *aPart)
Definition: G4VEmAdjointModel.cc:729
G4VEmAdjointModel::G4VEmAdjointModel
G4VEmAdjointModel(const G4String &nam)
Definition: G4VEmAdjointModel.cc:42
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:360
G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
Definition: G4VEmAdjointModel.cc:82
G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix
G4double SampleAdjSecEnergyFromCSMatrix(std::size_t MatrixIndex, G4double prim_energy, G4bool isScatProjToProj)
Definition: G4VEmAdjointModel.cc:413
G4VEmAdjointModel::fElementCSProdToProj
std::vector< G4double > fElementCSProdToProj
Definition: G4VEmAdjointModel.hh:268
G4VEmAdjointModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double aVal)
Definition: G4VEmAdjointModel.cc:713
G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
Definition: G4VEmAdjointModel.cc:252
G4VEmAdjointModel::GetHighEnergyLimit
G4double GetHighEnergyLimit()
Definition: G4VEmAdjointModel.hh:152
G4VEmAdjointModel::fLastAdjointCSForProdToProj
G4double fLastAdjointCSForProdToProj
Definition: G4VEmAdjointModel.hh:275
G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy, G4bool isScatProjToProj)
Definition: G4VEmAdjointModel.cc:557
G4VEmModel::SetHighEnergyLimit
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:746
G4VEmModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:284
G4VEmModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:181
G4VEmModel::SetLowEnergyLimit
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:753
G4VParticleChange::SetSecondaryWeightByProcess
void SetSecondaryWeightByProcess(G4bool)
G4VParticleChange::SetParentWeightByProcess
void SetParentWeightByProcess(G4bool)
Definition: G4VParticleChange.cc:375
G4VParticleChange::ProposeParentWeight
void ProposeParentWeight(G4double finalWeight)
G4lrint
int G4lrint(double ad)
Definition: templates.hh:134