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);
269 while (E1 <maxEProj*0.9999999){
272 int_cross_section +=integral.
Simpson(
this,
274 log_ESec_vector->push_back(std::log(
std::min(E2,maxEProj)));
275 log_Prob_vector->push_back(std::log(int_cross_section));
280 std::vector< std::vector<G4double>* > res_mat;
282 if (int_cross_section >0.) {
283 res_mat.push_back(log_ESec_vector);
284 res_mat.push_back(log_Prob_vector);
296 G4int nbin_pro_decade)
307 G4double dEmax=maxEProj-kinEnergyScatProj;
313 std::vector< double>* log_ESec_vector =
new std::vector< double>();
314 std::vector< double>* log_Prob_vector =
new std::vector< double>();
315 log_ESec_vector->push_back(std::log(dEmin));
316 log_Prob_vector->push_back(-50.);
317 G4int nbins=
std::max(
int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
318 G4double fE=std::pow(dEmax/dEmin,1./nbins);
327 while (dE1 <dEmax*0.9999999999999){
329 int_cross_section +=integral.
Simpson(
this,
332 log_ESec_vector->push_back(std::log(
std::min(dE2,maxEProj-minEProj)));
333 log_Prob_vector->push_back(std::log(int_cross_section));
339 std::vector< std::vector<G4double> *> res_mat;
341 if (int_cross_section >0.) {
342 res_mat.push_back(log_ESec_vector);
343 res_mat.push_back(log_Prob_vector);
353 G4int nbin_pro_decade)
363 std::vector< double>* log_ESec_vector =
new std::vector< double>();
364 std::vector< double>* log_Prob_vector =
new std::vector< double>();
365 log_ESec_vector->clear();
366 log_Prob_vector->clear();
367 log_ESec_vector->push_back(std::log(E1));
368 log_Prob_vector->push_back(-50.);
370 G4double E2=std::pow(10.,
double(
int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
371 G4double fE=std::pow(10.,1./nbin_pro_decade);
374 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
377 while (E1 <maxEProj*0.9999999){
379 int_cross_section +=integral.
Simpson(
this,
381 log_ESec_vector->push_back(std::log(
std::min(E2,maxEProj)));
382 log_Prob_vector->push_back(std::log(int_cross_section));
387 std::vector< std::vector<G4double>* > res_mat;
390 if (int_cross_section >0.) {
391 res_mat.push_back(log_ESec_vector);
392 res_mat.push_back(log_Prob_vector);
405 G4int nbin_pro_decade)
417 G4double dEmax=maxEProj-kinEnergyScatProj;
423 std::vector< double>* log_ESec_vector =
new std::vector< double>();
424 std::vector< double>* log_Prob_vector =
new std::vector< double>();
425 log_ESec_vector->push_back(std::log(dEmin));
426 log_Prob_vector->push_back(-50.);
427 G4int nbins=
std::max(
int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
428 G4double fE=std::pow(dEmax/dEmin,1./nbins);
433 while (dE1 <dEmax*0.9999999999999){
435 int_cross_section +=integral.
Simpson(
this,
437 log_ESec_vector->push_back(std::log(
std::min(dE2,maxEProj-minEProj)));
438 log_Prob_vector->push_back(std::log(int_cross_section));
447 std::vector< std::vector<G4double> *> res_mat;
449 if (int_cross_section >0.) {
450 res_mat.push_back(log_ESec_vector);
451 res_mat.push_back(log_Prob_vector);
462 G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
463 if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
466 if (theLogPrimEnergyVector->size() ==0){
467 G4cout<<
"No data are contained in the given AdjointCSMatrix!"<<
G4endl;
468 G4cout<<
"The sampling procedure will be stopped."<<
G4endl;
474 G4double aLogPrimEnergy = std::log(aPrimEnergy);
478 G4double aLogPrimEnergy1,aLogPrimEnergy2;
481 std::vector< double>* aLogSecondEnergyVector1 =0;
482 std::vector< double>* aLogSecondEnergyVector2 =0;
483 std::vector< double>* aLogProbVector1=0;
484 std::vector< double>* aLogProbVector2=0;
485 std::vector< size_t>* aLogProbVectorIndex1=0;
486 std::vector< size_t>* aLogProbVectorIndex2=0;
488 theMatrix->
GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
489 theMatrix->
GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
492 G4double log_rand_var= std::log(rand_var);
496 G4double log_rand_var1,log_rand_var2;
498 log_rand_var1=log_rand_var;
499 log_rand_var2=log_rand_var;
511 log_rand_var1=log_rand_var+theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
512 log_rand_var2=log_rand_var+theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
515 log_dE1 = theInterpolator->
Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,
"Lin");
516 log_dE2 = theInterpolator->
Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,
"Lin");
517 dE=std::exp(theInterpolator->
LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
520 Esec = aPrimEnergy +
dE;
527 log_E1 = theInterpolator->
Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,
"Lin");
528 log_E2 = theInterpolator->
Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,
"Lin");
530 Esec = std::exp(theInterpolator->
LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
561 if ( !IsScatProjToProjCase) {
568 for (
size_t i=0;i<CS_Vs_Element->size();i++){
569 SumCS+=(*CS_Vs_Element)[i];
570 if (rand_var<=SumCS/
lastCS){
585 const G4int iimax = 1000;
588 if ( IsScatProjToProjCase){
597 x = 1./(q*(1./xmin -1.) +1.);
601 if(ii >= iimax) {
break; }
615 x = std::pow(xmin, q);
619 if(ii >= iimax) {
break; }
637 G4bool IsScatProjToProjCase)
648 ,IsScatProjToProjCase );
649 if (post_stepCS>0 &&
lastCS>0) w_corr*=post_stepCS/
lastCS;
655 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
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
double A(double temperature)
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)
const G4double x[NPOINTSGL]
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 const G4double Emin
static G4AdjointInterpolator * GetInstance()
static G4Electron * Electron()
static const G4double Emax
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