133 if (anAdjPartDef && aProcess){
148 if (anAdjPartDef && aProcess){
207 G4cout<<
"========== Computation of cross section matrices for adjoint models =========="<<
G4endl;
212 std::vector<G4AdjointCSMatrix*>* aListOfMat1 =
new std::vector<G4AdjointCSMatrix*>();
213 std::vector<G4AdjointCSMatrix*>* aListOfMat2 =
new std::vector<G4AdjointCSMatrix*>();
214 aListOfMat1->clear();
215 aListOfMat2->clear();
218 std::vector<G4AdjointCSMatrix*>
220 aListOfMat1->push_back(two_matrices[0]);
221 aListOfMat2->push_back(two_matrices[1]);
224 for (
size_t j=0; j<theElementTable->size();j++){
225 G4Element* anElement=(*theElementTable)[j];
228 std::vector<G4AdjointCSMatrix*>
230 aListOfMat1->push_back(two_matrices[0]);
231 aListOfMat2->push_back(two_matrices[1]);
236 for (
size_t j=0; j<theMaterialTable->size();j++){
238 std::vector<G4AdjointCSMatrix*>
240 aListOfMat1->push_back(two_matrices[0]);
241 aListOfMat2->push_back(two_matrices[1]);
249 else {
G4cout<<
"The model "<<aModel->
GetName()<<
" does not use cross section matrices"<<
G4endl;
250 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
256 G4cout<<
" All adjoint cross section matrices are computed!"<<
G4endl;
257 G4cout<<
"======================================================================"<<
G4endl;
329 size_t mat_index = couple->
GetIndex();
338 if (totCS>sigma_max){
345 if (totCS>0 && !Emin_found) {
370 if (totCS>sigma_max){
376 if (totCS>0 && !Emin_found) {
512 corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
536 G4bool IsScatProjToProjCase,
537 std::vector<G4double>& CS_Vs_Element)
543 if (IsScatProjToProjCase){
551 if (EminSec >= EmaxSec)
return 0.;
554 G4bool need_to_compute=
false;
562 need_to_compute=
true;
566 if (!need_to_compute){
567 need_to_compute=
true;
571 need_to_compute=
false;
578 if (need_to_compute){
590 CS_Vs_Element.clear();
600 if (IsScatProjToProjCase){
605 if (PrimEnergy > Tlow)
608 for (
size_t i=0;i<n_el;i++){
613 CS_Vs_Element.push_back(CS);
617 for (
size_t i=0;i<n_el;i++){
621 if (IsScatProjToProjCase){
626 if (PrimEnergy > Tlow)
635 size_t ind_mat = aMaterial->
GetIndex();
637 if (IsScatProjToProjCase){
642 if (PrimEnergy > Tlow)
644 CS_Vs_Element.push_back(CS);
654 for (
size_t i=0;i<CS_Vs_Element.size();i++){
655 CS+=CS_Vs_Element[i];
666 G4bool IsScatProjToProjCase)
667 { std::vector<G4double> CS_Vs_Element;
672 for (
size_t i=0;i<CS_Vs_Element.size();i++){
673 SumCS+=CS_Vs_Element[i];
674 if (rand_var<=SumCS/CS){
696 std::vector<G4double> CS_Vs_Element;
717 if (aPartDef ==
listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
720 Ekin, Tlow,
true,CS_Vs_Element);
724 if (aPartDef ==
listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
727 Ekin, Tlow,
false, CS_Vs_Element);
745 std::vector<G4AdjointCSMatrix*>
747 G4int nbin_pro_decade)
763 G4double fE=std::pow(10.,1./nbin_pro_decade);
764 G4double E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
766 while (E1 <EkinMaxForProd){
770 if (aMat.size()>=2) {
771 std::vector< double>* log_ESecVec=aMat[0];
772 std::vector< double>* log_CSVec=aMat[1];
773 G4double log_adjointCS=log_CSVec->back();
775 for (
size_t j=0;j<log_CSVec->size();j++) {
776 if (j==0) (*log_CSVec)[j] = 0.;
777 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1e-50);
779 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
780 theCSMatForProdToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
789 E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
791 while (E1 <EkinMaxForScat){
795 if (aMat.size()>=2) {
796 std::vector< double>* log_ESecVec=aMat[0];
797 std::vector< double>* log_CSVec=aMat[1];
798 G4double log_adjointCS=log_CSVec->back();
800 for (
size_t j=0;j<log_CSVec->size();j++) {
801 if (j==0) (*log_CSVec)[j] = 0.;
802 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1e-50);
804 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
805 theCSMatForScatProjToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
812 std::vector<G4AdjointCSMatrix*> res;
814 res.push_back(theCSMatForProdToProjBackwardScattering);
815 res.push_back(theCSMatForScatProjToProjBackwardScattering);
836 std::vector<G4AdjointCSMatrix*>
839 G4int nbin_pro_decade)
860 G4double fE=std::pow(10.,1./nbin_pro_decade);
861 G4double E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
863 while (E1 <EkinMaxForProd){
867 if (aMat.size()>=2) {
868 std::vector< double>* log_ESecVec=aMat[0];
869 std::vector< double>* log_CSVec=aMat[1];
870 G4double log_adjointCS=log_CSVec->back();
873 for (
size_t j=0;j<log_CSVec->size();j++) {
875 if (j==0) (*log_CSVec)[j] = 0.;
876 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
879 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
880 theCSMatForProdToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
892 E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
894 while (E1 <EkinMaxForScat){
898 if (aMat.size()>=2) {
899 std::vector< double>* log_ESecVec=aMat[0];
900 std::vector< double>* log_CSVec=aMat[1];
901 G4double log_adjointCS=log_CSVec->back();
903 for (
size_t j=0;j<log_CSVec->size();j++) {
905 if (j==0) (*log_CSVec)[j] = 0.;
906 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
910 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
912 theCSMatForScatProjToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
924 std::vector<G4AdjointCSMatrix*> res;
927 res.push_back(theCSMatForProdToProjBackwardScattering);
928 res.push_back(theCSMatForScatProjToProjBackwardScattering);
1001 if (theLogPrimEnergyVector->size() ==0){
1002 G4cout<<
"No data are contained in the given AdjointCSMatrix!"<<
G4endl;
1007 G4double log_Tcut = std::log(Tcut);
1008 G4double log_E =std::log(aPrimEnergy);
1010 if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back())
return 0.;
1017 G4double aLogPrimEnergy1,aLogPrimEnergy2;
1020 std::vector< double>* aLogSecondEnergyVector1 =0;
1021 std::vector< double>* aLogSecondEnergyVector2 =0;
1022 std::vector< double>* aLogProbVector1=0;
1023 std::vector< double>* aLogProbVector2=0;
1024 std::vector< size_t>* aLogProbVectorIndex1=0;
1025 std::vector< size_t>* aLogProbVectorIndex2=0;
1028 anAdjointCSMatrix->
GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
1029 anAdjointCSMatrix->
GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
1031 G4double log_minimum_prob1, log_minimum_prob2;
1032 log_minimum_prob1=theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
1033 log_minimum_prob2=theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
1034 aLogCS1+= log_minimum_prob1;
1035 aLogCS2+= log_minimum_prob2;
1039 return std::exp(log_adjointCS);
static G4AdjointGamma * AdjointGamma()
G4double GetTotalAdjointCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
G4double GetCrossSectionCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, const G4MaterialCutsCouple *aCouple, G4bool &fwd_is_used, G4double &fwd_TotCS)
G4double GetAdjointSigma(G4double Ekin_nuc, size_t index_model, G4bool is_scat_proj_to_proj, const G4MaterialCutsCouple *aCouple)
G4ParticleDefinition * theFwdIon
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4double ComputeTotalAdjointCS(const G4MaterialCutsCouple *aMatCutCouple, G4ParticleDefinition *aPart, G4double PrimEnergy)
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 GetPostStepWeightCorrection()
G4double LastCSCorrectionFactor
std::vector< G4PhysicsTable * > theTotalAdjointSigmaTableVector
G4MaterialCutsCouple * currentCouple
std::vector< G4PhysicsTable * > listSigmaTableForAdjointModelProdToProj
void GetMaxFwdTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
size_t RegisterEmAdjointModel(G4VEmAdjointModel *)
void DefineCurrentMaterial(const G4MaterialCutsCouple *couple)
static G4MaterialTable * GetMaterialTable()
std::vector< G4Material * > G4MaterialTable
static G4AdjointElectron * AdjointElectron()
G4bool GetUseOnlyOneMatrixForAllElements()
G4bool forward_CS_is_used
void RegisterEnergyLossProcess(G4VEnergyLossProcess *aProcess, G4ParticleDefinition *aPartDef)
size_t GetVectorLength() const
G4double GetLowEdgeEnergy(size_t binNumber) const
std::vector< size_t > listOfIndexOfAdjointEMModelInAction
const G4Element * GetElement(G4int iel) const
void AddData(G4double aPrimEnergy, G4double aCS, std::vector< double > *aLogSecondEnergyVector, std::vector< double > *aLogProbVector, size_t n_pro_decade=0)
const G4String & GetParticleName() const
std::vector< G4PhysicsTable * > listSigmaTableForAdjointModelScatProjToProj
virtual G4double GetSecondAdjEnergyMaxForScatProjToProjCase(G4double PrimAdjEnergy)
std::vector< G4PhysicsTable * > theTotalForwardSigmaTableVector
G4double InterpolateForLogVector(G4double &x, std::vector< G4double > &x_vec, std::vector< G4double > &y_vec)
void RegisterAdjointParticle(G4ParticleDefinition *aPartDef)
std::vector< std::vector< G4AdjointCSMatrix * > > theAdjointCSMatricesForProdToProj
std::vector< std::vector< G4double > > EminForFwdSigmaTables
const G4double * GetVecNbOfAtomsPerVolume() const
static G4ThreadLocal G4AdjointCSManager * theInstance
G4GLOB_DLL std::ostream G4cout
G4bool GetUseMatrixPerElement()
virtual G4double GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
size_t GetTableSize() const
G4ParticleDefinition * GetAdjointParticleEquivalent(G4ParticleDefinition *theFwdPartDef)
std::vector< std::vector< G4double > > EkinofAdjSigmaMax
std::vector< std::vector< G4AdjointCSMatrix * > > theAdjointCSMatricesForScatProjToProj
std::vector< G4ParticleDefinition * > theListOfAdjointParticlesInAction
G4ParticleDefinition * theAdjIon
std::vector< double > * GetLogPrimEnergyVector()
std::vector< std::vector< G4VEnergyLossProcess * > * > listOfForwardEnergyLossProcess
void PutValue(size_t index, G4double theValue)
size_t FindPositionForLogVector(G4double &x, std::vector< G4double > &x_vec)
std::vector< G4AdjointCSMatrix * > BuildCrossSectionsMatricesForAGivenModelAndMaterial(G4VEmAdjointModel *aModel, G4Material *aMaterial, G4int nbin_pro_decade)
static G4Proton * Proton()
G4double lastPrimaryEnergy
G4double GetLowEnergyLimit()
void SetCSMatrices(std::vector< G4AdjointCSMatrix * > *Vec1CSMatrix, std::vector< G4AdjointCSMatrix * > *Vec2CSMatrix)
G4double GetTotalForwardCS(G4ParticleDefinition *aPartDef, G4double Ekin, const G4MaterialCutsCouple *aCouple)
static const G4double A[nN]
std::vector< std::vector< G4double > > EkinofFwdSigmaMax
std::vector< std::vector< G4double > > EminForAdjSigmaTables
void GetEminForTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &emin_adj, G4double &emin_fwd)
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerAtomForScatProj(G4double kinEnergyProd, G4double Z, G4double A=0., G4int nbin_pro_decade=10)
G4bool IsScatProjToProjCase()
G4ParticleDefinition * GetForwardParticleEquivalent(G4ParticleDefinition *theAdjPartDef)
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
virtual G4double AdjointCrossSection(const G4MaterialCutsCouple *aCouple, G4double primEnergy, G4bool IsScatProjToProjCase)
static G4ProductionCutsTable * GetProductionCutsTable()
G4double GetContinuousWeightCorrection(G4ParticleDefinition *aPartDef, G4double PreStepEkin, G4double AfterStepEkin, const G4MaterialCutsCouple *aCouple, G4double step_length)
G4ParticleDefinition * lastPartDefForCS
void RegisterEmProcess(G4VEmProcess *aProcess, G4ParticleDefinition *aPartDef)
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4bool > listOfIsScatProjToProjCase
static G4AdjointProton * AdjointProton()
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 GetHighEnergyLimit()
G4double LinearInterpolation(G4double &x, G4double &x1, G4double &x2, G4double &y1, G4double &y2)
void BuildTotalSigmaTables()
std::vector< std::vector< G4VEmProcess * > * > listOfForwardEmProcess
void GetMaxAdjTotalCS(G4ParticleDefinition *aPartDef, const G4MaterialCutsCouple *aCouple, G4double &e_sigma_max, G4double &sigma_max)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4Material * lastMaterial
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForSecond(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
G4Material * currentMaterial
G4bool GetApplyCutInRange()
G4ParticleDefinition * currentParticleDef
static G4AdjointInterpolator * GetInstance()
G4bool GetSecondPartOfSameType()
std::vector< G4VEmAdjointModel * > listOfAdjointEMModel
static G4Electron * Electron()
std::vector< std::vector< G4double > > lastAdjointCSVsModelsAndElements
size_t GetNumberOfElements() const
std::vector< std::vector< double > * > ComputeAdjointCrossSectionVectorPerVolumeForScatProj(G4Material *aMaterial, G4double kinEnergyProd, G4int nbin_pro_decade=10)
void BuildCrossSectionMatrices()
std::vector< G4Element * > G4ElementTable
G4Element * SampleElementFromCSMatrices(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase)
static G4ElementTable * GetElementTable()
size_t currentParticleIndex
void DefineCurrentParticle(const G4ParticleDefinition *aPartDef)
static G4AdjointCSManager * GetAdjointCSManager()
std::vector< G4AdjointCSMatrix * > BuildCrossSectionsMatricesForAGivenModelAndElement(G4VEmAdjointModel *aModel, G4int Z, G4int A, G4int nbin_pro_decade)
const G4Material * GetMaterial() const
G4bool CrossSectionMatrixesAreBuilt
G4bool TotalSigmaTableAreBuilt
G4double ComputeAdjointCS(G4Material *aMaterial, G4VEmAdjointModel *aModel, G4double PrimEnergy, G4double Tcut, G4bool IsScatProjToProjCase, std::vector< G4double > &AdjointCS_for_each_element)