135 if (anAdjPartDef && aProcess){
150 if (anAdjPartDef && aProcess){
209 G4cout<<
"========== Computation of cross section matrices for adjoint models =========="<<
G4endl;
214 std::vector<G4AdjointCSMatrix*>* aListOfMat1 =
new std::vector<G4AdjointCSMatrix*>();
215 std::vector<G4AdjointCSMatrix*>* aListOfMat2 =
new std::vector<G4AdjointCSMatrix*>();
216 aListOfMat1->clear();
217 aListOfMat2->clear();
220 std::vector<G4AdjointCSMatrix*>
222 aListOfMat1->push_back(two_matrices[0]);
223 aListOfMat2->push_back(two_matrices[1]);
226 for (
size_t j=0; j<theElementTable->size();j++){
227 G4Element* anElement=(*theElementTable)[j];
230 std::vector<G4AdjointCSMatrix*>
232 aListOfMat1->push_back(two_matrices[0]);
233 aListOfMat2->push_back(two_matrices[1]);
238 for (
size_t j=0; j<theMaterialTable->size();j++){
240 std::vector<G4AdjointCSMatrix*>
242 aListOfMat1->push_back(two_matrices[0]);
243 aListOfMat2->push_back(two_matrices[1]);
251 else {
G4cout<<
"The model "<<aModel->
GetName()<<
" does not use cross section matrices"<<
G4endl;
252 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
258 G4cout<<
" All adjoint cross section matrices are computed!"<<
G4endl;
259 G4cout<<
"======================================================================"<<
G4endl;
331 size_t mat_index = couple->
GetIndex();
340 if (totCS>sigma_max){
347 if (totCS>0 && !Emin_found) {
372 if (totCS>sigma_max){
378 if (totCS>0 && !Emin_found) {
514 corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
538 G4bool IsScatProjToProjCase,
539 std::vector<G4double>& CS_Vs_Element)
545 if (IsScatProjToProjCase){
553 if (EminSec >= EmaxSec)
return 0.;
556 G4bool need_to_compute=
false;
564 need_to_compute=
true;
568 if (!need_to_compute){
569 need_to_compute=
true;
573 need_to_compute=
false;
580 if (need_to_compute){
592 CS_Vs_Element.clear();
602 if (IsScatProjToProjCase){
607 if (PrimEnergy > Tlow)
610 for (
size_t i=0;i<n_el;i++){
615 CS_Vs_Element.push_back(CS);
619 for (
size_t i=0;i<n_el;i++){
623 if (IsScatProjToProjCase){
628 if (PrimEnergy > Tlow)
637 size_t ind_mat = aMaterial->
GetIndex();
639 if (IsScatProjToProjCase){
644 if (PrimEnergy > Tlow)
646 CS_Vs_Element.push_back(CS);
656 for (
size_t i=0;i<CS_Vs_Element.size();i++){
657 CS+=CS_Vs_Element[i];
668 G4bool IsScatProjToProjCase)
669 { std::vector<G4double> CS_Vs_Element;
674 for (
size_t i=0;i<CS_Vs_Element.size();i++){
675 SumCS+=CS_Vs_Element[i];
676 if (rand_var<=SumCS/CS){
698 std::vector<G4double> CS_Vs_Element;
719 if (aPartDef ==
listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
722 Ekin, Tlow,
true,CS_Vs_Element);
726 if (aPartDef ==
listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
729 Ekin, Tlow,
false, CS_Vs_Element);
747 std::vector<G4AdjointCSMatrix*>
749 G4int nbin_pro_decade)
765 G4double fE=std::pow(10.,1./nbin_pro_decade);
766 G4double E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
769 while (E1 <EkinMaxForProd){
773 if (aMat.size()>=2) {
774 std::vector< double>* log_ESecVec=aMat[0];
775 std::vector< double>* log_CSVec=aMat[1];
776 G4double log_adjointCS=log_CSVec->back();
778 for (
size_t j=0;j<log_CSVec->size();j++) {
779 if (j==0) (*log_CSVec)[j] = 0.;
780 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1
e-50);
782 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
783 theCSMatForProdToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
792 E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
795 while (E1 <EkinMaxForScat){
799 if (aMat.size()>=2) {
800 std::vector< double>* log_ESecVec=aMat[0];
801 std::vector< double>* log_CSVec=aMat[1];
802 G4double log_adjointCS=log_CSVec->back();
804 for (
size_t j=0;j<log_CSVec->size();j++) {
805 if (j==0) (*log_CSVec)[j] = 0.;
806 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1
e-50);
808 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
809 theCSMatForScatProjToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
816 std::vector<G4AdjointCSMatrix*> res;
818 res.push_back(theCSMatForProdToProjBackwardScattering);
819 res.push_back(theCSMatForScatProjToProjBackwardScattering);
840 std::vector<G4AdjointCSMatrix*>
843 G4int nbin_pro_decade)
864 G4double fE=std::pow(10.,1./nbin_pro_decade);
865 G4double E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
868 while (E1 <EkinMaxForProd){
872 if (aMat.size()>=2) {
873 std::vector< double>* log_ESecVec=aMat[0];
874 std::vector< double>* log_CSVec=aMat[1];
875 G4double log_adjointCS=log_CSVec->back();
878 for (
size_t j=0;j<log_CSVec->size();j++) {
880 if (j==0) (*log_CSVec)[j] = 0.;
881 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
884 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
885 theCSMatForProdToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
897 E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
900 while (E1 <EkinMaxForScat){
904 if (aMat.size()>=2) {
905 std::vector< double>* log_ESecVec=aMat[0];
906 std::vector< double>* log_CSVec=aMat[1];
907 G4double log_adjointCS=log_CSVec->back();
909 for (
size_t j=0;j<log_CSVec->size();j++) {
911 if (j==0) (*log_CSVec)[j] = 0.;
912 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
916 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
918 theCSMatForScatProjToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
930 std::vector<G4AdjointCSMatrix*> res;
933 res.push_back(theCSMatForProdToProjBackwardScattering);
934 res.push_back(theCSMatForScatProjToProjBackwardScattering);
1007 if (theLogPrimEnergyVector->size() ==0){
1008 G4cout<<
"No data are contained in the given AdjointCSMatrix!"<<
G4endl;
1013 G4double log_Tcut = std::log(Tcut);
1014 G4double log_E =std::log(aPrimEnergy);
1016 if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back())
return 0.;
1023 G4double aLogPrimEnergy1,aLogPrimEnergy2;
1026 std::vector< double>* aLogSecondEnergyVector1 =0;
1027 std::vector< double>* aLogSecondEnergyVector2 =0;
1028 std::vector< double>* aLogProbVector1=0;
1029 std::vector< double>* aLogProbVector2=0;
1030 std::vector< size_t>* aLogProbVectorIndex1=0;
1031 std::vector< size_t>* aLogProbVectorIndex2=0;
1034 anAdjointCSMatrix->
GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
1035 anAdjointCSMatrix->
GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
1037 G4double log_minimum_prob1, log_minimum_prob2;
1038 log_minimum_prob1=theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
1039 log_minimum_prob2=theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
1040 aLogCS1+= log_minimum_prob1;
1041 aLogCS2+= log_minimum_prob2;
1045 return std::exp(log_adjointCS);
G4bool CrossSectionMatrixesAreBuilt
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4PhysicsTable * > listSigmaTableForAdjointModelScatProjToProj
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForSecond(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4bool GetSecondPartOfSameType()
void BuildTotalSigmaTables()
G4double GetLowEnergyLimit()
G4double GetAdjointSigma(G4double Ekin_nuc, size_t index_model, G4bool is_scat_proj_to_proj, const G4MaterialCutsCouple *aCouple)
std::vector< G4PhysicsTable * > theTotalForwardSigmaTableVector
virtual G4double GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy, G4double Tcut=0)
G4double LastCSCorrectionFactor
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4double GetLowEdgeEnergy(size_t binNumber) const
G4bool TotalSigmaTableAreBuilt
void SetCSMatrices(std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
static G4MaterialTable * GetMaterialTable()
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
G4bool GetApplyCutInRange()
std::vector< std::vector< G4VEnergyLossProcess * > * > listOfForwardEnergyLossProcess
static constexpr double keV
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
std::vector< std::vector< G4double > > EkinofFwdSigmaMax
const G4String & GetParticleName() const
G4bool forward_CS_is_used
static G4AdjointElectron * AdjointElectron()
G4double lastPrimaryEnergy
std::vector< std::vector< G4double > > lastAdjointCSVsModelsAndElements
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
std::vector< std::vector< G4VEmProcess * > * > listOfForwardEmProcess
G4ParticleDefinition * GetForwardParticleEquivalent(G4ParticleDefinition *theAdjPartDef)
G4double GetHighEnergyLimit()
const G4double * GetVecNbOfAtomsPerVolume() const
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
void RegisterEnergyLossProcess(G4VEnergyLossProcess *aProcess, G4ParticleDefinition *aPartDef)
static G4Proton * Proton()
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
G4double GetCrossSectionCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, const G4MaterialCutsCouple *aCouple, G4bool &fwd_is_used, G4double &fwd_TotCS)
G4double GetPDGMass() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4bool GetUseMatrixPerElement()
G4ParticleDefinition * GetAdjointParticleEquivalent(G4ParticleDefinition *theFwdPartDef)
static G4AdjointGamma * AdjointGamma()
G4MaterialCutsCouple * currentCouple
static constexpr double proton_mass_c2
G4double GetTotalAdjointCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
G4Material * currentMaterial
std::vector< G4PhysicsTable * > listSigmaTableForAdjointModelProdToProj
void AddData(G4double aPrimEnergy, G4double aCS, std::vector< double > *aLogSecondEnergyVector, std::vector< double > *aLogProbVector, size_t n_pro_decade=0)
std::vector< std::vector< G4AdjointCSMatrix * > > theAdjointCSMatricesForScatProjToProj
static constexpr double TeV
void GetMaxFwdTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
G4Material * lastMaterial
std::vector< G4AdjointCSMatrix * > BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel *aModel, G4Material *aMaterial, G4int nbin_pro_decade)
std::vector< G4PhysicsTable * > theTotalAdjointSigmaTableVector
G4double GetPostStepWeightCorrection()
void DefineCurrentParticle(const G4ParticleDefinition *aPartDef)
std::vector< G4Element * > G4ElementTable
virtual G4double GetSecondAdjEnergyMaxForProdToProjCase(G4double PrimAdjEnergy)
G4bool GetUseOnlyOneMatrixForAllElements()
std::vector< G4ParticleDefinition * > theListOfAdjointParticlesInAction
G4ParticleDefinition * currentParticleDef
double A(double temperature)
static G4AdjointProton * AdjointProton()
static G4ThreadLocal G4AdjointCSManager * theInstance
static G4AdjointInterpolator * GetInstance()
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
std::vector< std::vector< G4AdjointCSMatrix * > > theAdjointCSMatricesForProdToProj
G4ParticleDefinition * theFwdIon
void RegisterAdjointParticle(G4ParticleDefinition *aPartDef)
static G4Electron * Electron()
size_t GetTableSize() const
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< std::vector< G4double > > EminForAdjSigmaTables
static G4AdjointCSManager * GetAdjointCSManager()
std::vector< G4Material * > G4MaterialTable
std::vector< size_t > listOfIndexOfAdjointEMModelInAction
std::vector< G4bool > listOfIsScatProjToProjCase
G4double GetTotalForwardCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
static G4ProductionCutsTable * GetProductionCutsTable()
G4ParticleDefinition * lastPartDefForCS
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
size_t currentParticleIndex
G4Element * SampleElementFromCSMatrices(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase)
std::vector< std::vector< G4double > > EkinofAdjSigmaMax
G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple *aMatCutCouple, G4ParticleDefinition *aPart, G4double PrimEnergy)
G4ParticleDefinition * theAdjIon
G4GLOB_DLL std::ostream G4cout
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
void BuildCrossSectionMatrices()
G4bool IsScatProjToProjCase()
std::vector< G4AdjointCSMatrix * > BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel *aModel, G4int Z, G4int A, G4int nbin_pro_decade)
void RegisterEmProcess(G4VEmProcess *aProcess, G4ParticleDefinition *aPartDef)
const G4Material * GetMaterial() const
void GetEminForTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &emin_adj, G4double &emin_fwd)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
void GetMaxAdjTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
const G4Element * GetElement(G4int iel) const
G4double GetContinuousWeightCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, G4double AfterStepEkin, const G4MaterialCutsCouple *aCouple, G4double step_length)
std::vector< double > * GetLogPrimEnergyVector()
void PutValue(size_t index, G4double theValue)
std::vector< G4VEmAdjointModel * > listOfAdjointEMModel
size_t GetNumberOfElements() const
size_t GetVectorLength() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
std::vector< std::vector< G4double > > EminForFwdSigmaTables
static G4ElementTable * GetElementTable()