Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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