134 if (anAdjPartDef && aProcess){
149 if (anAdjPartDef && aProcess){
208 G4cout<<
"========== Computation of cross section matrices for adjoint models =========="<<
G4endl;
213 std::vector<G4AdjointCSMatrix*>* aListOfMat1 =
new std::vector<G4AdjointCSMatrix*>();
214 std::vector<G4AdjointCSMatrix*>* aListOfMat2 =
new std::vector<G4AdjointCSMatrix*>();
215 aListOfMat1->clear();
216 aListOfMat2->clear();
219 std::vector<G4AdjointCSMatrix*>
221 aListOfMat1->push_back(two_matrices[0]);
222 aListOfMat2->push_back(two_matrices[1]);
225 for (
size_t j=0; j<theElementTable->size();j++){
226 G4Element* anElement=(*theElementTable)[j];
229 std::vector<G4AdjointCSMatrix*>
231 aListOfMat1->push_back(two_matrices[0]);
232 aListOfMat2->push_back(two_matrices[1]);
237 for (
size_t j=0; j<theMaterialTable->size();j++){
239 std::vector<G4AdjointCSMatrix*>
241 aListOfMat1->push_back(two_matrices[0]);
242 aListOfMat2->push_back(two_matrices[1]);
250 else {
G4cout<<
"The model "<<aModel->
GetName()<<
" does not use cross section matrices"<<
G4endl;
251 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
257 G4cout<<
" All adjoint cross section matrices are computed!"<<
G4endl;
258 G4cout<<
"======================================================================"<<
G4endl;
330 size_t mat_index = couple->
GetIndex();
339 if (totCS>sigma_max){
346 if (totCS>0 && !Emin_found) {
371 if (totCS>sigma_max){
377 if (totCS>0 && !Emin_found) {
513 corr_fac *=std::exp((pre_adjCS-pre_fwdCS)*step_length);
537 G4bool IsScatProjToProjCase,
538 std::vector<G4double>& CS_Vs_Element)
544 if (IsScatProjToProjCase){
552 if (EminSec >= EmaxSec)
return 0.;
555 G4bool need_to_compute=
false;
563 need_to_compute=
true;
567 if (!need_to_compute){
568 need_to_compute=
true;
572 need_to_compute=
false;
579 if (need_to_compute){
591 CS_Vs_Element.clear();
601 if (IsScatProjToProjCase){
606 if (PrimEnergy > Tlow)
609 for (
size_t i=0;i<n_el;i++){
614 CS_Vs_Element.push_back(CS);
618 for (
size_t i=0;i<n_el;i++){
622 if (IsScatProjToProjCase){
627 if (PrimEnergy > Tlow)
636 size_t ind_mat = aMaterial->
GetIndex();
638 if (IsScatProjToProjCase){
643 if (PrimEnergy > Tlow)
645 CS_Vs_Element.push_back(CS);
655 for (
size_t i=0;i<CS_Vs_Element.size();i++){
656 CS+=CS_Vs_Element[i];
667 G4bool IsScatProjToProjCase)
668 { std::vector<G4double> CS_Vs_Element;
673 for (
size_t i=0;i<CS_Vs_Element.size();i++){
674 SumCS+=CS_Vs_Element[i];
675 if (rand_var<=SumCS/CS){
697 std::vector<G4double> CS_Vs_Element;
718 if (aPartDef ==
listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectPrimaryParticleDefinition()){
721 Ekin, Tlow,
true,CS_Vs_Element);
725 if (aPartDef ==
listOfAdjointEMModel[i]->GetAdjointEquivalentOfDirectSecondaryParticleDefinition()){
728 Ekin, Tlow,
false, CS_Vs_Element);
746 std::vector<G4AdjointCSMatrix*>
748 G4int nbin_pro_decade)
764 G4double fE=std::pow(10.,1./nbin_pro_decade);
765 G4double E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
768 while (E1 <EkinMaxForProd){
772 if (aMat.size()>=2) {
773 std::vector< double>* log_ESecVec=aMat[0];
774 std::vector< double>* log_CSVec=aMat[1];
775 G4double log_adjointCS=log_CSVec->back();
777 for (
size_t j=0;j<log_CSVec->size();j++) {
778 if (j==0) (*log_CSVec)[j] = 0.;
779 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS) +1
e-50);
781 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
782 theCSMatForProdToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
791 E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
794 while (E1 <EkinMaxForScat){
798 if (aMat.size()>=2) {
799 std::vector< double>* log_ESecVec=aMat[0];
800 std::vector< double>* log_CSVec=aMat[1];
801 G4double log_adjointCS=log_CSVec->back();
803 for (
size_t j=0;j<log_CSVec->size();j++) {
804 if (j==0) (*log_CSVec)[j] = 0.;
805 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS)+1
e-50);
807 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
808 theCSMatForScatProjToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
815 std::vector<G4AdjointCSMatrix*> res;
817 res.push_back(theCSMatForProdToProjBackwardScattering);
818 res.push_back(theCSMatForScatProjToProjBackwardScattering);
839 std::vector<G4AdjointCSMatrix*>
842 G4int nbin_pro_decade)
863 G4double fE=std::pow(10.,1./nbin_pro_decade);
864 G4double E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
867 while (E1 <EkinMaxForProd){
871 if (aMat.size()>=2) {
872 std::vector< double>* log_ESecVec=aMat[0];
873 std::vector< double>* log_CSVec=aMat[1];
874 G4double log_adjointCS=log_CSVec->back();
877 for (
size_t j=0;j<log_CSVec->size();j++) {
879 if (j==0) (*log_CSVec)[j] = 0.;
880 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
883 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
884 theCSMatForProdToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
896 E2=std::pow(10.,
double(
int(std::log10(EkinMin)*nbin_pro_decade)+1)/nbin_pro_decade)/fE;
899 while (E1 <EkinMaxForScat){
903 if (aMat.size()>=2) {
904 std::vector< double>* log_ESecVec=aMat[0];
905 std::vector< double>* log_CSVec=aMat[1];
906 G4double log_adjointCS=log_CSVec->back();
908 for (
size_t j=0;j<log_CSVec->size();j++) {
910 if (j==0) (*log_CSVec)[j] = 0.;
911 else (*log_CSVec)[j]=std::log(1.-std::exp((*log_CSVec)[j]-log_adjointCS));
915 (*log_CSVec)[log_CSVec->size()-1]=(*log_CSVec)[log_CSVec->size()-2]-std::log(1000.);
917 theCSMatForScatProjToProjBackwardScattering->
AddData(std::log(E1),log_adjointCS,log_ESecVec,log_CSVec,0);
929 std::vector<G4AdjointCSMatrix*> res;
932 res.push_back(theCSMatForProdToProjBackwardScattering);
933 res.push_back(theCSMatForScatProjToProjBackwardScattering);
1006 if (theLogPrimEnergyVector->size() ==0){
1007 G4cout<<
"No data are contained in the given AdjointCSMatrix!"<<
G4endl;
1012 G4double log_Tcut = std::log(Tcut);
1013 G4double log_E =std::log(aPrimEnergy);
1015 if (aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back())
return 0.;
1022 G4double aLogPrimEnergy1,aLogPrimEnergy2;
1025 std::vector< double>* aLogSecondEnergyVector1 =0;
1026 std::vector< double>* aLogSecondEnergyVector2 =0;
1027 std::vector< double>* aLogProbVector1=0;
1028 std::vector< double>* aLogProbVector2=0;
1029 std::vector< size_t>* aLogProbVectorIndex1=0;
1030 std::vector< size_t>* aLogProbVectorIndex2=0;
1033 anAdjointCSMatrix->
GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
1034 anAdjointCSMatrix->
GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
1036 G4double log_minimum_prob1, log_minimum_prob2;
1037 log_minimum_prob1=theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
1038 log_minimum_prob2=theInterpolator->
InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
1039 aLogCS1+= log_minimum_prob1;
1040 aLogCS2+= log_minimum_prob2;
1044 return std::exp(log_adjointCS);