61 G4bool IsScatProjToProjCase)
86 G4bool IsScatProjToProjCase)
89 IsScatProjToProjCase);
124 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
135 dSigmadEprod=(sigma1-sigma2)/dE;
150 {
G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
152 if (kinEnergyProd <=0) dSigmadEprod=0;
171 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
179 dSigmadEprod=(sigma1-sigma2)/dE;
193 {
G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
195 if (kinEnergyProd <=0) dSigmadEprod=0;
243 G4int nbin_pro_decade)
256 std::vector< double>* log_ESec_vector =
new std::vector< double>();
257 std::vector< double>* log_Prob_vector =
new std::vector< double>();
258 log_ESec_vector->clear();
259 log_Prob_vector->clear();
260 log_ESec_vector->push_back(std::log(E1));
261 log_Prob_vector->push_back(-50.);
263 G4double E2=std::pow(10.,
double(
int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
264 G4double fE=std::pow(10.,1./nbin_pro_decade);
267 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
270 while (E1 <maxEProj*0.9999999){
273 int_cross_section +=integral.
Simpson(
this,
275 log_ESec_vector->push_back(std::log(
std::min(E2,maxEProj)));
276 log_Prob_vector->push_back(std::log(int_cross_section));
281 std::vector< std::vector<G4double>* > res_mat;
283 if (int_cross_section >0.) {
284 res_mat.push_back(log_ESec_vector);
285 res_mat.push_back(log_Prob_vector);
297 G4int nbin_pro_decade)
308 G4double dEmax=maxEProj-kinEnergyScatProj;
314 std::vector< double>* log_ESec_vector =
new std::vector< double>();
315 std::vector< double>* log_Prob_vector =
new std::vector< double>();
316 log_ESec_vector->push_back(std::log(dEmin));
317 log_Prob_vector->push_back(-50.);
318 G4int nbins=
std::max(
int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
319 G4double fE=std::pow(dEmax/dEmin,1./nbins);
328 while (dE1 <dEmax*0.9999999999999){
330 int_cross_section +=integral.
Simpson(
this,
333 log_ESec_vector->push_back(std::log(
std::min(dE2,maxEProj-minEProj)));
334 log_Prob_vector->push_back(std::log(int_cross_section));
340 std::vector< std::vector<G4double> *> res_mat;
342 if (int_cross_section >0.) {
343 res_mat.push_back(log_ESec_vector);
344 res_mat.push_back(log_Prob_vector);
354 G4int nbin_pro_decade)
364 std::vector< double>* log_ESec_vector =
new std::vector< double>();
365 std::vector< double>* log_Prob_vector =
new std::vector< double>();
366 log_ESec_vector->clear();
367 log_Prob_vector->clear();
368 log_ESec_vector->push_back(std::log(E1));
369 log_Prob_vector->push_back(-50.);
371 G4double E2=std::pow(10.,
double(
int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
372 G4double fE=std::pow(10.,1./nbin_pro_decade);
375 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
378 while (E1 <maxEProj*0.9999999){
380 int_cross_section +=integral.
Simpson(
this,
382 log_ESec_vector->push_back(std::log(
std::min(E2,maxEProj)));
383 log_Prob_vector->push_back(std::log(int_cross_section));
388 std::vector< std::vector<G4double>* > res_mat;
391 if (int_cross_section >0.) {
392 res_mat.push_back(log_ESec_vector);
393 res_mat.push_back(log_Prob_vector);
406 G4int nbin_pro_decade)
418 G4double dEmax=maxEProj-kinEnergyScatProj;
424 std::vector< double>* log_ESec_vector =
new std::vector< double>();
425 std::vector< double>* log_Prob_vector =
new std::vector< double>();
426 log_ESec_vector->push_back(std::log(dEmin));
427 log_Prob_vector->push_back(-50.);
428 G4int nbins=
std::max(
int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
429 G4double fE=std::pow(dEmax/dEmin,1./nbins);
434 while (dE1 <dEmax*0.9999999999999){
436 int_cross_section +=integral.
Simpson(
this,
438 log_ESec_vector->push_back(std::log(
std::min(dE2,maxEProj-minEProj)));
439 log_Prob_vector->push_back(std::log(int_cross_section));
448 std::vector< std::vector<G4double> *> res_mat;
450 if (int_cross_section >0.) {
451 res_mat.push_back(log_ESec_vector);
452 res_mat.push_back(log_Prob_vector);
463 G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
464 if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
467 if (theLogPrimEnergyVector->size() ==0){
468 G4cout<<
"No data are contained in the given AdjointCSMatrix!"<<
G4endl;
469 G4cout<<
"The sampling procedure will be stopped."<<
G4endl;
475 G4double aLogPrimEnergy = std::log(aPrimEnergy);
479 G4double aLogPrimEnergy1,aLogPrimEnergy2;
482 std::vector< double>* aLogSecondEnergyVector1 =0;
483 std::vector< double>* aLogSecondEnergyVector2 =0;
484 std::vector< double>* aLogProbVector1=0;
485 std::vector< double>* aLogProbVector2=0;
486 std::vector< size_t>* aLogProbVectorIndex1=0;
487 std::vector< size_t>* aLogProbVectorIndex2=0;
489 theMatrix->
GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
490 theMatrix->
GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
493 G4double log_rand_var= std::log(rand_var);
497 G4double log_rand_var1,log_rand_var2;
499 log_rand_var1=log_rand_var;
500 log_rand_var2=log_rand_var;
512 log_rand_var1=log_rand_var+theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
513 log_rand_var2=log_rand_var+theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
516 log_dE1 = theInterpolator->
Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,
"Lin");
517 log_dE2 = theInterpolator->
Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,
"Lin");
518 dE=std::exp(theInterpolator->
LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
521 Esec = aPrimEnergy +
dE;
528 log_E1 = theInterpolator->
Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,
"Lin");
529 log_E2 = theInterpolator->
Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,
"Lin");
531 Esec = std::exp(theInterpolator->
LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
562 if ( !IsScatProjToProjCase) {
569 for (
size_t i=0;i<CS_Vs_Element->size();i++){
570 SumCS+=(*CS_Vs_Element)[i];
571 if (rand_var<=SumCS/
lastCS){
586 const G4int iimax = 1000;
589 if ( IsScatProjToProjCase){
598 x = 1./(q*(1./xmin -1.) +1.);
602 if(ii >= iimax) {
break; }
616 x = std::pow(xmin, q);
620 if(ii >= iimax) {
break; }
638 G4bool IsScatProjToProjCase)
649 ,IsScatProjToProjCase );
650 if (post_stepCS>0 &&
lastCS>0) w_corr*=post_stepCS/
lastCS;
656 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)
G4double DiffCrossSectionFunction1(G4double kinEnergyProj)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4VEmAdjointModel(const G4String &nam)
G4double lastAdjointCSForScatProjToProjCase
static const G4double Emin
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4bool second_part_of_same_type
G4double GetLowEnergyLimit()
size_t currentCoupleIndex
G4double DiffCrossSectionFunction2(G4double kinEnergyProj)
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
void SetHighEnergyLimit(G4double)
G4double CS_biasing_factor
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
G4Material * currentMaterial
const G4String & GetParticleName() const
G4double Interpolate(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec, G4String InterPolMethod="Log")
static G4AdjointPositron * AdjointPositron()
static G4AdjointElectron * AdjointElectron()
void SetSecondaryWeightByProcess(G4bool)
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
G4bool UseOnlyOneMatrixForAllElements
G4double lastAdjointCSForProdToProjCase
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double additional_weight_correction_factor_for_post_step_outside_model
static G4AdjointGamma * AdjointGamma()
G4double currentTcutForDirectSecond
G4MaterialCutsCouple * currentCouple
G4double kinEnergyProjForIntegration
G4double kinEnergyProdForIntegration
G4ParticleDefinition * theAdjEquivOfDirectSecondPartDef
G4ParticleDefinition * theAdjEquivOfDirectPrimPartDef
virtual G4double DiffCrossSectionPerAtomPrimToSecond(G4double kinEnergyProj, G4double kinEnergyProd, G4double Z, G4double A=0.)
virtual G4double DiffCrossSectionPerVolumePrimToSecond(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyProd)
G4double GetPostStepWeightCorrection()
void SelectCSMatrix(G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
double A(double temperature)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
static G4AdjointInterpolator * GetInstance()
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
static G4Electron * Electron()
size_t currentMaterialIndex
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 Simpson(T &typeT, F f, G4double a, G4double b, G4int n)
static G4AdjointCSManager * GetAdjointCSManager()
size_t indexOfUsedCrossSectionMatrix
virtual void CorrectPostStepWeight(G4ParticleChange *fParticleChange, G4double old_weight, G4double adjointPrimKinEnergy, G4double projectileKinEnergy, G4bool IsScatProjToProjCase)
static G4ProductionCutsTable * GetProductionCutsTable()
G4VEmModel * theDirectEMModel
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual G4double DiffCrossSectionPerVolumePrimToScatPrim(const G4Material *aMaterial, G4double kinEnergyProj, G4double kinEnergyScatProj)
void SetParentWeightByProcess(G4bool)
void ProposeParentWeight(G4double finalWeight)
void SetLowEnergyLimit(G4double)
std::vector< G4double > CS_Vs_ElementForProdToProjCase
G4double DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double EkinProd)
static const G4double Emax
virtual G4double SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy, G4bool IsScatProjToProjCase)
void SetLowEnergyLimit(G4double aVal)
G4GLOB_DLL std::ostream G4cout
G4double kinEnergyScatProjForIntegration
virtual G4double DiffCrossSectionPerAtomPrimToScatPrim(G4double kinEnergyProj, G4double kinEnergyScatProj, G4double Z, G4double A=0.)
G4double SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex, G4double prim_energy, G4bool IsScatProjToProjCase)
G4bool IsScatProjToProjCase()
void SetHighEnergyLimit(G4double aVal)
const G4Material * GetMaterial() const
G4ParticleDefinition * theDirectPrimaryPartDef
G4Material * SelectedMaterial
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
const G4Element * GetElement(G4int iel) const
std::vector< double > * GetLogPrimEnergyVector()
G4double mass_ratio_projectile
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
virtual G4double GetAdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
std::vector< G4double > CS_Vs_ElementForScatProjToProjCase
void SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition *aPart)
G4bool UseMatrixPerElement
G4double mass_ratio_product
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
virtual ~G4VEmAdjointModel()