90 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
93 G4double E2 = kinEnergyProd * 1.000001;
99 dSigmadEprod = (sigma1 - sigma2) / (E2 - E1);
108 G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
110 if(kinEnergyProd <= 0.)
126 if(kinEnergyProj > Emin_proj && kinEnergyProj <= Emax_proj)
129 G4double E2 = kinEnergyProd * 1.0001;
134 dSigmadEprod = (sigma1 - sigma2) / (E2 - E1);
144 G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
146 if(kinEnergyProd <= 0.)
150 aMaterial, kinEnergyProj, kinEnergyProd);
197std::vector<std::vector<G4double>*>
200 G4int nbin_pro_decade)
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.);
218 std::pow(10.,
G4double(
G4int(std::log10(minEProj) * nbin_pro_decade) + 1) /
220 G4double fE = std::pow(10., 1. / nbin_pro_decade);
222 if(std::pow(fE, 5.) > (maxEProj / minEProj))
223 fE = std::pow(maxEProj / minEProj, 0.2);
227 while(E1 < maxEProj * 0.9999999)
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));
237 std::vector<std::vector<G4double>*> res_mat;
238 if(int_cross_section > 0.)
240 res_mat.push_back(log_ESec_vector);
241 res_mat.push_back(log_Prob_vector);
244 delete log_ESec_vector;
245 delete log_Prob_vector;
251std::vector<std::vector<G4double>*>
254 G4int nbin_pro_decade)
265 G4double dEmax = maxEProj - kinEnergyScatProj;
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.);
275 G4double fE = std::pow(dEmax / dEmin, 1. / nbins);
279 while(dE1 < dEmax * 0.9999999999999)
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));
290 std::vector<std::vector<G4double>*> res_mat;
291 if(int_cross_section > 0.)
293 res_mat.push_back(log_ESec_vector);
294 res_mat.push_back(log_Prob_vector);
297 delete log_ESec_vector;
298 delete log_Prob_vector;
305std::vector<std::vector<G4double>*>
308 G4int nbin_pro_decade)
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.);
325 std::pow(10.,
G4double(
G4int(std::log10(minEProj) * nbin_pro_decade) + 1) /
327 G4double fE = std::pow(10., 1. / nbin_pro_decade);
329 if(std::pow(fE, 5.) > (maxEProj / minEProj))
330 fE = std::pow(maxEProj / minEProj, 0.2);
334 while(E1 < maxEProj * 0.9999999)
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));
344 std::vector<std::vector<G4double>*> res_mat;
346 if(int_cross_section > 0.)
348 res_mat.push_back(log_ESec_vector);
349 res_mat.push_back(log_Prob_vector);
352 delete log_ESec_vector;
353 delete log_Prob_vector;
359std::vector<std::vector<G4double>*>
362 G4int nbin_pro_decade)
373 G4double dEmax = maxEProj - kinEnergyScatProj;
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);
387 while(dE1 < dEmax * 0.9999999999999)
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));
398 std::vector<std::vector<G4double>*> res_mat;
399 if(int_cross_section > 0.)
401 res_mat.push_back(log_ESec_vector);
402 res_mat.push_back(log_Prob_vector);
405 delete log_ESec_vector;
406 delete log_Prob_vector;
414 size_t MatrixIndex,
G4double aPrimEnergy,
G4bool isScatProjToProj)
418 theMatrix = (*fCSMatrixProjToProjBackScat)[MatrixIndex];
419 std::vector<G4double>* theLogPrimEnergyVector =
422 if(theLogPrimEnergyVector->empty())
424 G4cout <<
"No data are contained in the given AdjointCSMatrix!" <<
G4endl;
425 G4cout <<
"The sampling procedure will be stopped." <<
G4endl;
430 G4double aLogPrimEnergy = std::log(aPrimEnergy);
432 aLogPrimEnergy, *theLogPrimEnergyVector);
434 G4double aLogPrimEnergy1, aLogPrimEnergy2;
437 std::vector<double>* aLogSecondEnergyVector1 =
nullptr;
438 std::vector<double>* aLogSecondEnergyVector2 =
nullptr;
439 std::vector<double>* aLogProbVector1 =
nullptr;
440 std::vector<double>* aLogProbVector2 =
nullptr;
441 std::vector<size_t>* aLogProbVectorIndex1 =
nullptr;
442 std::vector<size_t>* aLogProbVectorIndex2 =
nullptr;
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);
451 if (! (aLogProbVector1 && aLogProbVector2 &&
452 aLogSecondEnergyVector1 && aLogSecondEnergyVector2)){
457 G4double log_rand_var = std::log(rand_var);
461 G4double log_rand_var1, log_rand_var2;
463 log_rand_var1 = log_rand_var;
464 log_rand_var2 = log_rand_var;
480 log_rand_var1 = log_rand_var +
482 log_Tcut, *aLogSecondEnergyVector1, *aLogProbVector1);
483 log_rand_var2 = log_rand_var +
485 log_Tcut, *aLogSecondEnergyVector2, *aLogProbVector2);
487 log_dE1 = theInterpolator->
Interpolate(log_rand_var1, *aLogProbVector1,
488 *aLogSecondEnergyVector1,
"Lin");
489 log_dE2 = theInterpolator->
Interpolate(log_rand_var2, *aLogProbVector2,
490 *aLogSecondEnergyVector2,
"Lin");
492 aLogPrimEnergy, aLogPrimEnergy1, aLogPrimEnergy2, log_dE1, log_dE2));
495 Esec = aPrimEnergy +
dE;
502 log_E1 = theInterpolator->
Interpolate(log_rand_var, *aLogProbVector1,
503 *aLogSecondEnergyVector1,
"Lin");
504 log_E2 = theInterpolator->
Interpolate(log_rand_var, *aLogProbVector2,
505 *aLogSecondEnergyVector2,
"Lin");
508 aLogPrimEnergy, aLogPrimEnergy1, aLogPrimEnergy2, log_E1, log_E2));
536 if(!isScatProjToProj)
543 for(
size_t i = 0; i < CS_Vs_Element->size(); ++i)
545 SumCS += (*CS_Vs_Element)[i];
561 constexpr G4int iimax = 1000;
627 if(!isScatProjToProj)
633 if(post_stepCS > 0. &&
fLastCS > 0.)
634 w_corr *= post_stepCS /
fLastCS;
637 new_weight *= w_corr;
638 new_weight *= projectileKinEnergy / adjointPrimKinEnergy;
653 maxEProj =
std::min(kinEnergyScatProj * 2., maxEProj);
679 minEProj = primAdjEnergy * 2.;
703 const std::vector<G4double>* aVec =
static const G4double Emax
static const G4double Emin
G4GLOB_DLL std::ostream G4cout
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool isScatProjToProj, std::vector< G4double > &AdjointCS_for_each_element)
G4double GetPostStepWeightCorrection()
static G4AdjointCSManager * GetAdjointCSManager()
G4bool IsScatProjToProj()
G4bool GetData(unsigned int i, G4double &aPrimEnergy, G4double &aCS, G4double &log0, std::vector< double > *&aLogSecondEnergyVector, std::vector< double > *&aLogProbVector, std::vector< size_t > *&aLogProbVectorIndex)
std::vector< double > * GetLogPrimEnergyVector()
static G4AdjointElectron * AdjointElectron()
static G4AdjointGamma * AdjointGamma()
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
G4double Interpolate(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec, G4String InterPolMethod="Log")
static G4AdjointInterpolator * GetInstance()
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
static G4AdjointPositron * AdjointPositron()
static G4Electron * Electron()
const G4Material * GetMaterial() const
const G4Element * GetElement(G4int iel) const
const G4String & GetParticleName() const
const std::vector< G4double > * GetEnergyCutsVector(std::size_t pcIdx) const
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double GetSecondAdjEnergyMaxForProdToProj(G4double primAdjEnergy)
G4double fLastAdjointCSForScatProjToProj
std::vector< G4AdjointCSMatrix * > * fCSMatrixProdToProjBackScat
G4bool fSecondPartSameType
G4double fCsBiasingFactor
G4ParticleDefinition * fAdjEquivDirectSecondPart
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
virtual ~G4VEmAdjointModel()
G4double fKinEnergyScatProjForIntegration
G4double DiffCrossSectionFunction1(G4double kinEnergyProj)
virtual G4double GetSecondAdjEnergyMinForScatProjToProj(G4double primAdjEnergy, G4double tcut=0.)
G4double fKinEnergyProdForIntegration
G4ParticleDefinition * fAdjEquivDirectPrimPart
G4MaterialCutsCouple * fCurrentCouple
virtual G4double GetSecondAdjEnergyMaxForScatProjToProj(G4double primAdjEnergy)
G4Material * fCurrentMaterial
std::vector< G4AdjointCSMatrix * > * fCSMatrixProjToProjBackScat
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool isScatProjToProj)
std::vector< G4double > fElementCSScatProjToProj
G4double GetLowEnergyLimit()
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool isScatProjToProj)
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
G4AdjointCSManager * fCSManager
void SetLowEnergyLimit(G4double aVal)
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool fUseMatrixPerElement
G4double fHighEnergyLimit
void SelectCSMatrix(G4bool isScatProjToProj)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool isScatProjToProj)
G4VEmModel * fDirectModel
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
G4Material * fSelectedMaterial
virtual G4double GetSecondAdjEnergyMinForProdToProj(G4double primAdjEnergy)
G4double DiffCrossSectionFunction2(G4double kinEnergyProj)
G4ParticleDefinition * fDirectPrimaryPart
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition *aPart)
G4VEmAdjointModel(const G4String &nam)
G4bool fOneMatrixForAllElements
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
std::vector< G4double > fElementCSProdToProj
void SetHighEnergyLimit(G4double aVal)
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4double GetHighEnergyLimit()
G4double fLastAdjointCSForProdToProj
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy, G4bool isScatProjToProj)
void SetHighEnergyLimit(G4double)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
void SetLowEnergyLimit(G4double)
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
void ProposeParentWeight(G4double finalWeight)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments