60 G4bool IsScatProjToProjCase)
85 G4bool IsScatProjToProjCase)
88 IsScatProjToProjCase);
123 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
134 dSigmadEprod=(sigma1-sigma2)/dE;
149 {
G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
151 if (kinEnergyProd <=0) dSigmadEprod=0;
170 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
178 dSigmadEprod=(sigma1-sigma2)/dE;
192 {
G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
194 if (kinEnergyProd <=0) dSigmadEprod=0;
242 G4int nbin_pro_decade)
255 std::vector< double>* log_ESec_vector =
new std::vector< double>();
256 std::vector< double>* log_Prob_vector =
new std::vector< double>();
257 log_ESec_vector->clear();
258 log_Prob_vector->clear();
259 log_ESec_vector->push_back(std::log(E1));
260 log_Prob_vector->push_back(-50.);
262 G4double E2=std::pow(10.,
double(
int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
263 G4double fE=std::pow(10.,1./nbin_pro_decade);
266 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
268 while (E1 <maxEProj*0.9999999){
271 int_cross_section +=integral.
Simpson(
this,
273 log_ESec_vector->push_back(std::log(
std::min(E2,maxEProj)));
274 log_Prob_vector->push_back(std::log(int_cross_section));
279 std::vector< std::vector<G4double>* > res_mat;
281 if (int_cross_section >0.) {
282 res_mat.push_back(log_ESec_vector);
283 res_mat.push_back(log_Prob_vector);
295 G4int nbin_pro_decade)
306 G4double dEmax=maxEProj-kinEnergyScatProj;
312 std::vector< double>* log_ESec_vector =
new std::vector< double>();
313 std::vector< double>* log_Prob_vector =
new std::vector< double>();
314 log_ESec_vector->push_back(std::log(dEmin));
315 log_Prob_vector->push_back(-50.);
316 G4int nbins=
std::max(
int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
317 G4double fE=std::pow(dEmax/dEmin,1./nbins);
325 while (dE1 <dEmax*0.9999999999999){
327 int_cross_section +=integral.
Simpson(
this,
330 log_ESec_vector->push_back(std::log(
std::min(dE2,maxEProj-minEProj)));
331 log_Prob_vector->push_back(std::log(int_cross_section));
337 std::vector< std::vector<G4double> *> res_mat;
339 if (int_cross_section >0.) {
340 res_mat.push_back(log_ESec_vector);
341 res_mat.push_back(log_Prob_vector);
351 G4int nbin_pro_decade)
361 std::vector< double>* log_ESec_vector =
new std::vector< double>();
362 std::vector< double>* log_Prob_vector =
new std::vector< double>();
363 log_ESec_vector->clear();
364 log_Prob_vector->clear();
365 log_ESec_vector->push_back(std::log(E1));
366 log_Prob_vector->push_back(-50.);
368 G4double E2=std::pow(10.,
double(
int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
369 G4double fE=std::pow(10.,1./nbin_pro_decade);
372 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
374 while (E1 <maxEProj*0.9999999){
376 int_cross_section +=integral.
Simpson(
this,
378 log_ESec_vector->push_back(std::log(
std::min(E2,maxEProj)));
379 log_Prob_vector->push_back(std::log(int_cross_section));
384 std::vector< std::vector<G4double>* > res_mat;
387 if (int_cross_section >0.) {
388 res_mat.push_back(log_ESec_vector);
389 res_mat.push_back(log_Prob_vector);
402 G4int nbin_pro_decade)
414 G4double dEmax=maxEProj-kinEnergyScatProj;
420 std::vector< double>* log_ESec_vector =
new std::vector< double>();
421 std::vector< double>* log_Prob_vector =
new std::vector< double>();
422 log_ESec_vector->push_back(std::log(dEmin));
423 log_Prob_vector->push_back(-50.);
424 G4int nbins=
std::max(
int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
425 G4double fE=std::pow(dEmax/dEmin,1./nbins);
429 while (dE1 <dEmax*0.9999999999999){
431 int_cross_section +=integral.
Simpson(
this,
433 log_ESec_vector->push_back(std::log(
std::min(dE2,maxEProj-minEProj)));
434 log_Prob_vector->push_back(std::log(int_cross_section));
443 std::vector< std::vector<G4double> *> res_mat;
445 if (int_cross_section >0.) {
446 res_mat.push_back(log_ESec_vector);
447 res_mat.push_back(log_Prob_vector);
458 G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
459 if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
462 if (theLogPrimEnergyVector->size() ==0){
463 G4cout<<
"No data are contained in the given AdjointCSMatrix!"<<
G4endl;
464 G4cout<<
"The sampling procedure will be stopped."<<
G4endl;
470 G4double aLogPrimEnergy = std::log(aPrimEnergy);
474 G4double aLogPrimEnergy1,aLogPrimEnergy2;
477 std::vector< double>* aLogSecondEnergyVector1 =0;
478 std::vector< double>* aLogSecondEnergyVector2 =0;
479 std::vector< double>* aLogProbVector1=0;
480 std::vector< double>* aLogProbVector2=0;
481 std::vector< size_t>* aLogProbVectorIndex1=0;
482 std::vector< size_t>* aLogProbVectorIndex2=0;
484 theMatrix->
GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
485 theMatrix->
GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
488 G4double log_rand_var= std::log(rand_var);
492 G4double log_rand_var1,log_rand_var2;
494 log_rand_var1=log_rand_var;
495 log_rand_var2=log_rand_var;
507 log_rand_var1=log_rand_var+theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
508 log_rand_var2=log_rand_var+theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
511 log_dE1 = theInterpolator->
Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,
"Lin");
512 log_dE2 = theInterpolator->
Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,
"Lin");
513 dE=std::exp(theInterpolator->
LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
516 Esec = aPrimEnergy +dE;
523 log_E1 = theInterpolator->
Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,
"Lin");
524 log_E2 = theInterpolator->
Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,
"Lin");
526 Esec = std::exp(theInterpolator->
LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
557 if ( !IsScatProjToProjCase) {
564 for (
size_t i=0;i<CS_Vs_Element->size();i++){
565 SumCS+=(*CS_Vs_Element)[i];
566 if (rand_var<=SumCS/
lastCS){
583 if ( IsScatProjToProjCase){
591 x = 1./(q*(1./xmin -1.) +1.);
607 x = std::pow(xmin, q);
628 G4bool IsScatProjToProjCase)
639 ,IsScatProjToProjCase );
640 if (post_stepCS>0 &&
lastCS>0) w_corr*=post_stepCS/
lastCS;
646 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
static G4AdjointGamma * AdjointGamma()
G4double kinEnergyScatProjForIntegration
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual ~G4VEmAdjointModel()
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4bool GetData(unsigned int i, G4double &aPrimEnergy, G4double &aCS, G4double &log0, std::vector< double > *&aLogSecondEnergyVector, std::vector< double > *&aLogProbVector, std::vector< size_t > *&aLogProbVectorIndex)
G4double lastAdjointCSForScatProjToProjCase
G4double GetPostStepWeightCorrection()
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4ParticleDefinition * theDirectPrimaryPartDef
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
void ProposeParentWeight(G4double finalWeight)
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy, G4bool IsScatProjToProjCase)
G4double DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double EkinProd)
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
static G4AdjointElectron * AdjointElectron()
G4double CS_biasing_factor
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
G4Material * currentMaterial
const G4Element * GetElement(G4int iel) const
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
void SetLowEnergyLimit(G4double aVal)
G4Material * SelectedMaterial
size_t indexOfUsedCrossSectionMatrix
const G4String & GetParticleName() const
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void SetHighEnergyLimit(G4double)
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
void SetHighEnergyLimit(G4double aVal)
G4double DiffCrossSectionFunction1(G4double kinEnergyProj)
G4double Interpolate(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec, G4String InterPolMethod="Log")
void SetSecondaryWeightByProcess(G4bool)
void SetParentWeightByProcess(G4bool)
G4VEmModel * theDirectEMModel
G4GLOB_DLL std::ostream G4cout
G4double mass_ratio_product
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
std::vector< G4double > CS_Vs_ElementForProdToProjCase
G4VEmAdjointModel(const G4String &nam)
std::vector< G4double > CS_Vs_ElementForScatProjToProjCase
std::vector< double > * GetLogPrimEnergyVector()
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
G4double GetLowEnergyLimit()
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition *aPart)
G4double DiffCrossSectionFunction2(G4double kinEnergyProj)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
G4MaterialCutsCouple * currentCouple
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
G4bool UseOnlyOneMatrixForAllElements
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4bool IsScatProjToProjCase()
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
static G4ProductionCutsTable * GetProductionCutsTable()
G4double mass_ratio_projectile
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
G4double kinEnergyProjForIntegration
void SelectCSMatrix(G4bool IsScatProjToProjCase)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
static G4AdjointInterpolator * GetInstance()
static G4Electron * Electron()
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
size_t currentCoupleIndex
G4double currentTcutForDirectSecond
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
G4double lastAdjointCSForProdToProjCase
G4double Simpson(T &typeT, F f, G4double a, G4double b, G4int n)
G4bool UseMatrixPerElement
void SetLowEnergyLimit(G4double)
static G4AdjointCSManager * GetAdjointCSManager()
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
G4double kinEnergyProdForIntegration
size_t currentMaterialIndex
static G4AdjointPositron * AdjointPositron()
const G4Material * GetMaterial() const
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)
G4bool second_part_of_same_type