| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 << 26 // $Id: G4VEmAdjointModel.cc 66892 2013-01-17 10:57:59Z gunter $
>> 27 //
27 #include "G4VEmAdjointModel.hh" 28 #include "G4VEmAdjointModel.hh"
28 <<
29 #include "G4AdjointCSManager.hh" 29 #include "G4AdjointCSManager.hh"
>> 30 #include "G4Integrator.hh"
>> 31 #include "G4TrackStatus.hh"
>> 32 #include "G4ParticleChange.hh"
30 #include "G4AdjointElectron.hh" 33 #include "G4AdjointElectron.hh"
31 #include "G4AdjointGamma.hh" 34 #include "G4AdjointGamma.hh"
32 #include "G4AdjointInterpolator.hh" <<
33 #include "G4AdjointPositron.hh" 35 #include "G4AdjointPositron.hh"
34 #include "G4Electron.hh" << 36 #include "G4AdjointInterpolator.hh"
35 #include "G4Gamma.hh" << 37 #include "G4PhysicsTable.hh"
36 #include "G4Integrator.hh" <<
37 #include "G4ParticleChange.hh" <<
38 #include "G4ProductionCutsTable.hh" <<
39 #include "G4TrackStatus.hh" <<
40 38
41 ////////////////////////////////////////////// 39 ////////////////////////////////////////////////////////////////////////////////
42 G4VEmAdjointModel::G4VEmAdjointModel(const G4S << 40 //
43 : fName(nam) << 41 G4VEmAdjointModel::G4VEmAdjointModel(const G4String& nam):
44 { << 42 name(nam)
45 fCSManager = G4AdjointCSManager::GetAdjointC << 43 // lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
46 fCSManager->RegisterEmAdjointModel(this); << 44 {
>> 45 model_index = G4AdjointCSManager::GetAdjointCSManager()->RegisterEmAdjointModel(this);
>> 46 second_part_of_same_type =false;
>> 47 theDirectEMModel=0;
>> 48 mass_ratio_product=1.;
>> 49 mass_ratio_projectile=1.;
>> 50 currentCouple=0;
47 } 51 }
48 << 52 ////////////////////////////////////////////////////////////////////////////////
49 ///////////////////////////////////e////////// << 53 //
50 G4VEmAdjointModel::~G4VEmAdjointModel() 54 G4VEmAdjointModel::~G4VEmAdjointModel()
51 { << 55 {;}
52 delete fCSMatrixProdToProjBackScat; <<
53 fCSMatrixProdToProjBackScat = nullptr; <<
54 <<
55 delete fCSMatrixProjToProjBackScat; <<
56 fCSMatrixProjToProjBackScat = nullptr; <<
57 } <<
58 <<
59 ////////////////////////////////////////////// 56 ////////////////////////////////////////////////////////////////////////////////
60 G4double G4VEmAdjointModel::AdjointCrossSectio << 57 //
61 const G4MaterialCutsCouple* aCouple, G4doubl << 58 G4double G4VEmAdjointModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple,
62 G4bool isScatProjToProj) << 59 G4double primEnergy,
63 { << 60 G4bool IsScatProjToProjCase)
>> 61 {
64 DefineCurrentMaterial(aCouple); 62 DefineCurrentMaterial(aCouple);
65 fPreStepEnergy = primEnergy; << 63 preStepEnergy=primEnergy;
66 << 64
67 std::vector<G4double>* CS_Vs_Element = &fEle << 65 std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
68 if(isScatProjToProj) << 66 if (IsScatProjToProjCase) CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
69 CS_Vs_Element = &fElementCSScatProjToProj; << 67 lastCS = G4AdjointCSManager::GetAdjointCSManager()->ComputeAdjointCS(currentMaterial,
70 fLastCS = << 68 this,
71 fCSManager->ComputeAdjointCS(fCurrentMater << 69 primEnergy,
72 fTcutSecond, << 70 currentTcutForDirectSecond,
73 if(isScatProjToProj) << 71 IsScatProjToProjCase,
74 fLastAdjointCSForScatProjToProj = fLastCS; << 72 *CS_Vs_Element);
75 else << 73 if (IsScatProjToProjCase) lastAdjointCSForScatProjToProjCase = lastCS;
76 fLastAdjointCSForProdToProj = fLastCS; << 74 else lastAdjointCSForProdToProjCase =lastCS;
77 << 75
78 return fLastCS; << 76
79 } << 77
80 << 78 return lastCS;
>> 79
>> 80 }
>> 81 ////////////////////////////////////////////////////////////////////////////////
>> 82 //
>> 83 G4double G4VEmAdjointModel::GetAdjointCrossSection(const G4MaterialCutsCouple* aCouple,
>> 84 G4double primEnergy,
>> 85 G4bool IsScatProjToProjCase)
>> 86 {
>> 87 return AdjointCrossSection(aCouple, primEnergy,
>> 88 IsScatProjToProjCase);
>> 89
>> 90 /*
>> 91 //To continue
>> 92 DefineCurrentMaterial(aCouple);
>> 93 preStepEnergy=primEnergy;
>> 94 if (IsScatProjToProjCase){
>> 95 G4double ekin=primEnergy*mass_ratio_projectile;
>> 96 lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,true, aCouple);
>> 97 lastAdjointCSForScatProjToProjCase = lastCS;
>> 98 //G4cout<<ekin<<std::endl;
>> 99 }
>> 100 else {
>> 101 G4double ekin=primEnergy*mass_ratio_product;
>> 102 lastCS = G4AdjointCSManager::GetAdjointCSManager()->GetAdjointSigma(ekin, model_index,false, aCouple);
>> 103 lastAdjointCSForProdToProjCase = lastCS;
>> 104 //G4cout<<ekin<<std::endl;
>> 105 }
>> 106 return lastCS;
>> 107 */
>> 108 }
81 ////////////////////////////////////////////// 109 ////////////////////////////////////////////////////////////////////////////////
>> 110 //
>> 111 //General implementation correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
82 G4double G4VEmAdjointModel::DiffCrossSectionPe 112 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToSecond(
83 G4double kinEnergyProj, G4double kinEnergyPr << 113 G4double kinEnergyProj,
84 { << 114 G4double kinEnergyProd,
85 G4double dSigmadEprod = 0.; << 115 G4double Z,
86 G4double Emax_proj = GetSecondAdjEnergyMa << 116 G4double A)
87 G4double Emin_proj = GetSecondAdjEnergyMi << 117 {
88 << 118 G4double dSigmadEprod=0;
89 // the produced particle should have a kinet << 119 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
90 if(kinEnergyProj > Emin_proj && kinEnergyPro << 120 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
91 { << 121
92 G4double E1 = kinEnergyProd; << 122
93 G4double E2 = kinEnergyProd * 1.000001 << 123 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ //the produced particle should have a kinetic energy smaller than the projectile
94 G4double sigma1 = fDirectModel->ComputeCro << 124
95 fDirectPrimaryPart, kinEnergyProj, Z, A, << 125 /*G4double Tmax=kinEnergyProj;
96 G4double sigma2 = fDirectModel->ComputeCro << 126 if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
97 fDirectPrimaryPart, kinEnergyProj, Z, A, << 127
98 << 128 G4double E1=kinEnergyProd;
99 dSigmadEprod = (sigma1 - sigma2) / (E2 - E << 129 G4double E2=kinEnergyProd*1.000001;
100 } << 130 G4double dE=(E2-E1);
101 return dSigmadEprod; << 131 G4double sigma1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E1,1.e20);
>> 132 G4double sigma2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,Z,A ,E2,1.e20);
>> 133
>> 134 dSigmadEprod=(sigma1-sigma2)/dE;
>> 135 }
>> 136 return dSigmadEprod;
>> 137
>> 138
>> 139
102 } 140 }
103 << 141 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
104 ////////////////////////////////////////////// 142 ////////////////////////////////////////////////////////////////////////////////
>> 143 //
105 G4double G4VEmAdjointModel::DiffCrossSectionPe 144 G4double G4VEmAdjointModel::DiffCrossSectionPerAtomPrimToScatPrim(
106 G4double kinEnergyProj, G4double kinEnergySc << 145 G4double kinEnergyProj,
107 { << 146 G4double kinEnergyScatProj,
108 G4double kinEnergyProd = kinEnergyProj - kin << 147 G4double Z,
>> 148 G4double A)
>> 149 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
109 G4double dSigmadEprod; 150 G4double dSigmadEprod;
110 if(kinEnergyProd <= 0.) << 151 if (kinEnergyProd <=0) dSigmadEprod=0;
111 dSigmadEprod = 0.; << 152 else dSigmadEprod=DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProd,Z,A);
112 else << 153 return dSigmadEprod;
113 dSigmadEprod = << 154
114 DiffCrossSectionPerAtomPrimToSecond(kinE <<
115 return dSigmadEprod; <<
116 } 155 }
117 156
118 ////////////////////////////////////////////// 157 ////////////////////////////////////////////////////////////////////////////////
>> 158 //
>> 159 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
119 G4double G4VEmAdjointModel::DiffCrossSectionPe 160 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(
120 const G4Material* aMaterial, G4double kinEne << 161 const G4Material* aMaterial,
121 { << 162 G4double kinEnergyProj,
122 G4double dSigmadEprod = 0.; << 163 G4double kinEnergyProd)
123 G4double Emax_proj = GetSecondAdjEnergyMa << 164 {
124 G4double Emin_proj = GetSecondAdjEnergyMi << 165 G4double dSigmadEprod=0;
125 << 166 G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
126 if(kinEnergyProj > Emin_proj && kinEnergyPro << 167 G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
127 { << 168
128 G4double E1 = kinEnergyProd; << 169
129 G4double E2 = kinEnergyProd * 1.0001; << 170 if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){
130 G4double sigma1 = fDirectModel->CrossSecti << 171 /*G4double Tmax=kinEnergyProj;
131 aMaterial, fDirectPrimaryPart, kinEnergy << 172 if (second_part_of_same_type) Tmax = kinEnergyProj/2.;*/
132 G4double sigma2 = fDirectModel->CrossSecti << 173 G4double E1=kinEnergyProd;
133 aMaterial, fDirectPrimaryPart, kinEnergy << 174 G4double E2=kinEnergyProd*1.0001;
134 dSigmadEprod = (sigma1 - sigma2) / (E2 - E << 175 G4double dE=(E2-E1);
135 } << 176 G4double sigma1=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E1,1.e20);
136 return dSigmadEprod; << 177 G4double sigma2=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,E2,1.e20);
>> 178 dSigmadEprod=(sigma1-sigma2)/dE;
>> 179 }
>> 180 return dSigmadEprod;
>> 181
>> 182
>> 183
137 } 184 }
138 << 185 //The implementation here is correct for energy loss process, for the photoelectric and compton scattering the method should be redefine
139 ////////////////////////////////////////////// 186 ////////////////////////////////////////////////////////////////////////////////
>> 187 //
140 G4double G4VEmAdjointModel::DiffCrossSectionPe 188 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToScatPrim(
141 const G4Material* aMaterial, G4double kinEne << 189 const G4Material* aMaterial,
142 G4double kinEnergyScatProj) << 190 G4double kinEnergyProj,
143 { << 191 G4double kinEnergyScatProj)
144 G4double kinEnergyProd = kinEnergyProj - kin << 192 { G4double kinEnergyProd = kinEnergyProj - kinEnergyScatProj;
145 G4double dSigmadEprod; 193 G4double dSigmadEprod;
146 if(kinEnergyProd <= 0.) << 194 if (kinEnergyProd <=0) dSigmadEprod=0;
147 dSigmadEprod = 0.; << 195 else dSigmadEprod=DiffCrossSectionPerVolumePrimToSecond(aMaterial,kinEnergyProj,kinEnergyProd);
148 else << 196 return dSigmadEprod;
149 dSigmadEprod = DiffCrossSectionPerVolumePr << 197
150 aMaterial, kinEnergyProj, kinEnergyProd) << 198 }
151 return dSigmadEprod; << 199 ///////////////////////////////////////////////////////////////////////////////////////////////////////////
152 } << 200 //
153 << 201 G4double G4VEmAdjointModel::DiffCrossSectionFunction1(G4double kinEnergyProj){
154 ////////////////////////////////////////////// << 202
155 G4double G4VEmAdjointModel::DiffCrossSectionFu << 203
156 { << 204 G4double bias_factor = CS_biasing_factor*kinEnergyProdForIntegration/kinEnergyProj;
157 G4double bias_factor = << 205
158 fCsBiasingFactor * fKinEnergyProdForIntegr << 206
159 << 207 if (UseMatrixPerElement ) {
160 if(fUseMatrixPerElement) << 208 return DiffCrossSectionPerAtomPrimToSecond(kinEnergyProj,kinEnergyProdForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
161 { << 209 }
162 return DiffCrossSectionPerAtomPrimToSecond << 210 else {
163 kinEnergyProj, fKinEnergyProdForI << 211 return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProj,kinEnergyProdForIntegration)*bias_factor;
164 fASelectedNucleus) * << 212 }
165 bias_factor; << 213 }
166 } << 214
167 else << 215 ////////////////////////////////////////////////////////////////////////////////
168 { << 216 //
169 return DiffCrossSectionPerVolumePrimToSeco << 217 G4double G4VEmAdjointModel::DiffCrossSectionFunction2(G4double kinEnergyProj){
170 fSelectedMaterial, kinEnergyProj, << 218
171 bias_factor; << 219 G4double bias_factor = CS_biasing_factor*kinEnergyScatProjForIntegration/kinEnergyProj;
172 } << 220 if (UseMatrixPerElement ) {
173 } << 221 return DiffCrossSectionPerAtomPrimToScatPrim(kinEnergyProj,kinEnergyScatProjForIntegration,ZSelectedNucleus,ASelectedNucleus)*bias_factor;
174 << 222 }
175 ////////////////////////////////////////////// << 223 else {
176 G4double G4VEmAdjointModel::DiffCrossSectionFu << 224 return DiffCrossSectionPerVolumePrimToScatPrim(SelectedMaterial,kinEnergyProj,kinEnergyScatProjForIntegration)*bias_factor;
177 { << 225
178 G4double bias_factor = << 226 }
179 fCsBiasingFactor * fKinEnergyScatProjForIn << 227
180 if(fUseMatrixPerElement) << 228 }
181 { << 229 ////////////////////////////////////////////////////////////////////////////////
182 return DiffCrossSectionPerAtomPrimToScatPr << 230 //
183 kinEnergyProj, fKinEnergyScatProj << 231
184 fASelectedNucleus) * << 232 G4double G4VEmAdjointModel::DiffCrossSectionPerVolumeFunctionForIntegrationOverEkinProj(G4double kinEnergyProd)
185 bias_factor; << 233 {
186 } << 234 return DiffCrossSectionPerVolumePrimToSecond(SelectedMaterial,kinEnergyProjForIntegration,kinEnergyProd);
187 else << 235 }
188 { << 236 ////////////////////////////////////////////////////////////////////////////////
189 return DiffCrossSectionPerVolumePrimToScat << 237 //
190 fSelectedMaterial, kinEnergyProj, << 238 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForSecond(
191 fKinEnergyScatProjForIntegration) << 239 G4double kinEnergyProd,
192 bias_factor; << 240 G4double Z,
193 } << 241 G4double A ,
194 } << 242 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
195 << 243 {
196 ////////////////////////////////////////////// << 244 G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
197 std::vector<std::vector<G4double>*> << 245 ASelectedNucleus= int(A);
198 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 246 ZSelectedNucleus=int(Z);
199 G4double kinEnergyProd, G4double Z, G4double << 247 kinEnergyProdForIntegration = kinEnergyProd;
200 G4int nbin_pro_decade) // nb bins per order << 248
201 { << 249 //compute the vector of integrated cross sections
202 G4Integrator<G4VEmAdjointModel, G4double (G4 << 250 //-------------------
203 integral; << 251
204 fASelectedNucleus = G4lrint(A); << 252 G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
205 fZSelectedNucleus = G4lrint(Z); << 253 G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
206 fKinEnergyProdForIntegration = kinEnergyProd << 254 G4double E1=minEProj;
207 << 255 std::vector< double>* log_ESec_vector = new std::vector< double>();
208 // compute the vector of integrated cross se << 256 std::vector< double>* log_Prob_vector = new std::vector< double>();
209 G4double minEProj = GetSecondAdjEnergyMinFor << 257 log_ESec_vector->clear();
210 G4double maxEProj = GetSecondAdjEnergyMaxFor << 258 log_Prob_vector->clear();
211 G4double E1 = minEProj; <<
212 std::vector<G4double>* log_ESec_vector = new <<
213 std::vector<G4double>* log_Prob_vector = new <<
214 log_ESec_vector->push_back(std::log(E1)); 259 log_ESec_vector->push_back(std::log(E1));
215 log_Prob_vector->push_back(-50.); 260 log_Prob_vector->push_back(-50.);
216 << 261
217 G4double E2 = << 262 G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
218 std::pow(10., G4double(G4int(std::log10(mi << 263 G4double fE=std::pow(10.,1./nbin_pro_decade);
219 nbin_pro_decade); << 264 G4double int_cross_section=0.;
220 G4double fE = std::pow(10., 1. / nbin_pro_de << 265
221 << 266 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
222 if(std::pow(fE, 5.) > (maxEProj / minEProj)) << 267
223 fE = std::pow(maxEProj / minEProj, 0.2); << 268 while (E1 <maxEProj*0.9999999){
224 << 269 //G4cout<<E1<<'\t'<<E2<<G4endl;
225 G4double int_cross_section = 0.; << 270
226 // Loop checking, 07-Aug-2015, Vladimir Ivan << 271 int_cross_section +=integral.Simpson(this,
227 while(E1 < maxEProj * 0.9999999) << 272 &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
228 { << 273 log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
229 int_cross_section += << 274 log_Prob_vector->push_back(std::log(int_cross_section));
230 integral.Simpson(this, &G4VEmAdjointMode << 275 E1=E2;
231 std::min(E2, maxEProj * << 276 E2*=fE;
232 log_ESec_vector->push_back(std::log(std::m << 277
233 log_Prob_vector->push_back(std::log(int_cr << 278 }
234 E1 = E2; << 279 std::vector< std::vector<G4double>* > res_mat;
235 E2 *= fE; << 280 res_mat.clear();
236 } << 281 if (int_cross_section >0.) {
237 std::vector<std::vector<G4double>*> res_mat; << 282 res_mat.push_back(log_ESec_vector);
238 if(int_cross_section > 0.) << 283 res_mat.push_back(log_Prob_vector);
239 { << 284 }
240 res_mat.push_back(log_ESec_vector); << 285
241 res_mat.push_back(log_Prob_vector); <<
242 } <<
243 else { <<
244 delete log_ESec_vector; <<
245 delete log_Prob_vector; <<
246 } <<
247 return res_mat; 286 return res_mat;
248 } 287 }
249 288
250 ////////////////////////////////////////////// 289 /////////////////////////////////////////////////////////////////////////////////////
251 std::vector<std::vector<G4double>*> << 290 //
252 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 291 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
253 G4double kinEnergyScatProj, G4double Z, G4do << 292 G4double kinEnergyScatProj,
254 G4int nbin_pro_decade) // nb bins pro order << 293 G4double Z,
255 { << 294 G4double A ,
256 G4Integrator<G4VEmAdjointModel, G4double (G4 << 295 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
257 integral; << 296 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
258 fASelectedNucleus = G4lrint(A << 297 ASelectedNucleus=int(A);
259 fZSelectedNucleus = G4lrint(Z << 298 ZSelectedNucleus=int(Z);
260 fKinEnergyScatProjForIntegration = kinEnergy << 299 kinEnergyScatProjForIntegration = kinEnergyScatProj;
261 << 300
262 // compute the vector of integrated cross se << 301 //compute the vector of integrated cross sections
263 G4double minEProj = GetSecondAdjEnergyMinFor << 302 //-------------------
264 G4double maxEProj = GetSecondAdjEnergyMaxFor << 303
265 G4double dEmax = maxEProj - kinEnergyScat << 304 G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
266 G4double dEmin = GetLowEnergyLimit(); << 305 G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
267 G4double dE1 = dEmin; << 306 G4double dEmax=maxEProj-kinEnergyScatProj;
268 G4double dE2 = dEmin; << 307 G4double dEmin=GetLowEnergyLimit();
269 << 308 G4double dE1=dEmin;
270 std::vector<G4double>* log_ESec_vector = new << 309 G4double dE2=dEmin;
271 std::vector<G4double>* log_Prob_vector = new << 310
>> 311
>> 312 std::vector< double>* log_ESec_vector = new std::vector< double>();
>> 313 std::vector< double>* log_Prob_vector = new std::vector< double>();
272 log_ESec_vector->push_back(std::log(dEmin)); 314 log_ESec_vector->push_back(std::log(dEmin));
273 log_Prob_vector->push_back(-50.); 315 log_Prob_vector->push_back(-50.);
274 G4int nbins = std::max(G4int(std::log10(dEma << 316 G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
275 G4double fE = std::pow(dEmax / dEmin, 1. / n << 317 G4double fE=std::pow(dEmax/dEmin,1./nbins);
276 << 318
277 G4double int_cross_section = 0.; << 319
278 // Loop checking, 07-Aug-2015, Vladimir Ivan << 320
279 while(dE1 < dEmax * 0.9999999999999) << 321
280 { << 322
281 dE2 = dE1 * fE; << 323 G4double int_cross_section=0.;
282 int_cross_section += << 324
283 integral.Simpson(this, &G4VEmAdjointMode << 325 while (dE1 <dEmax*0.9999999999999){
284 minEProj + dE1, std::mi << 326 dE2=dE1*fE;
285 log_ESec_vector->push_back(std::log(std::m << 327 int_cross_section +=integral.Simpson(this,
286 log_Prob_vector->push_back(std::log(int_cr << 328 &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
287 dE1 = dE2; << 329 //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
288 } << 330 log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
289 << 331 log_Prob_vector->push_back(std::log(int_cross_section));
290 std::vector<std::vector<G4double>*> res_mat; << 332 dE1=dE2;
291 if(int_cross_section > 0.) << 333
292 { << 334 }
293 res_mat.push_back(log_ESec_vector); << 335
294 res_mat.push_back(log_Prob_vector); << 336
295 } << 337 std::vector< std::vector<G4double> *> res_mat;
296 else { << 338 res_mat.clear();
297 delete log_ESec_vector; << 339 if (int_cross_section >0.) {
298 delete log_Prob_vector; << 340 res_mat.push_back(log_ESec_vector);
299 } << 341 res_mat.push_back(log_Prob_vector);
300 << 342 }
>> 343
301 return res_mat; 344 return res_mat;
302 } 345 }
303 <<
304 ////////////////////////////////////////////// 346 ////////////////////////////////////////////////////////////////////////////////
305 std::vector<std::vector<G4double>*> << 347 //
306 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 348 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(
307 G4Material* aMaterial, G4double kinEnergyPro << 349 G4Material* aMaterial,
308 G4int nbin_pro_decade) // nb bins pro order << 350 G4double kinEnergyProd,
309 { << 351 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
310 G4Integrator<G4VEmAdjointModel, G4double (G4 << 352 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
311 integral; << 353 SelectedMaterial= aMaterial;
312 fSelectedMaterial = aMaterial; << 354 kinEnergyProdForIntegration = kinEnergyProd;
313 fKinEnergyProdForIntegration = kinEnergyProd << 355 //compute the vector of integrated cross sections
314 << 356 //-------------------
315 // compute the vector of integrated cross se << 357
316 G4double minEProj = GetSecondAdjEnergyMinFor << 358 G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
317 G4double maxEProj = GetSecondAdjEnergyMaxFor << 359 G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
318 G4double E1 = minEProj; << 360 G4double E1=minEProj;
319 std::vector<G4double>* log_ESec_vector = new << 361 std::vector< double>* log_ESec_vector = new std::vector< double>();
320 std::vector<G4double>* log_Prob_vector = new << 362 std::vector< double>* log_Prob_vector = new std::vector< double>();
>> 363 log_ESec_vector->clear();
>> 364 log_Prob_vector->clear();
321 log_ESec_vector->push_back(std::log(E1)); 365 log_ESec_vector->push_back(std::log(E1));
322 log_Prob_vector->push_back(-50.); 366 log_Prob_vector->push_back(-50.);
323 << 367
324 G4double E2 = << 368 G4double E2=std::pow(10.,double( int(std::log10(minEProj)*nbin_pro_decade)+1)/nbin_pro_decade);
325 std::pow(10., G4double(G4int(std::log10(mi << 369 G4double fE=std::pow(10.,1./nbin_pro_decade);
326 nbin_pro_decade); << 370 G4double int_cross_section=0.;
327 G4double fE = std::pow(10., 1. / nbin_pro_de << 371
328 << 372 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
329 if(std::pow(fE, 5.) > (maxEProj / minEProj)) << 373
330 fE = std::pow(maxEProj / minEProj, 0.2); << 374 while (E1 <maxEProj*0.9999999){
331 << 375
332 G4double int_cross_section = 0.; << 376 int_cross_section +=integral.Simpson(this,
333 // Loop checking, 07-Aug-2015, Vladimir Ivan << 377 &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
334 while(E1 < maxEProj * 0.9999999) << 378 log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
335 { << 379 log_Prob_vector->push_back(std::log(int_cross_section));
336 int_cross_section += << 380 E1=E2;
337 integral.Simpson(this, &G4VEmAdjointMode << 381 E2*=fE;
338 std::min(E2, maxEProj * << 382
339 log_ESec_vector->push_back(std::log(std::m << 383 }
340 log_Prob_vector->push_back(std::log(int_cr << 384 std::vector< std::vector<G4double>* > res_mat;
341 E1 = E2; << 385 res_mat.clear();
342 E2 *= fE; << 386
343 } << 387 if (int_cross_section >0.) {
344 std::vector<std::vector<G4double>*> res_mat; << 388 res_mat.push_back(log_ESec_vector);
345 << 389 res_mat.push_back(log_Prob_vector);
346 if(int_cross_section > 0.) << 390 }
347 { << 391
348 res_mat.push_back(log_ESec_vector); << 392
349 res_mat.push_back(log_Prob_vector); << 393
350 } <<
351 else { <<
352 delete log_ESec_vector; <<
353 delete log_Prob_vector; <<
354 } <<
355 return res_mat; 394 return res_mat;
356 } 395 }
357 396
358 ////////////////////////////////////////////// << 397 /////////////////////////////////////////////////////////////////////////////////////
359 std::vector<std::vector<G4double>*> << 398 //
360 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 399 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
361 G4Material* aMaterial, G4double kinEnergySca << 400 G4Material* aMaterial,
362 G4int nbin_pro_decade) // nb bins pro order << 401 G4double kinEnergyScatProj,
363 { << 402 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
364 G4Integrator<G4VEmAdjointModel, G4double (G4 << 403 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
365 integral; << 404 SelectedMaterial= aMaterial;
366 fSelectedMaterial = aMaterial << 405 kinEnergyScatProjForIntegration = kinEnergyScatProj;
367 fKinEnergyScatProjForIntegration = kinEnergy << 406
368 << 407 //compute the vector of integrated cross sections
369 // compute the vector of integrated cross se << 408 //-------------------
370 G4double minEProj = GetSecondAdjEnergyMinFor << 409
371 G4double maxEProj = GetSecondAdjEnergyMaxFor << 410 G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
372 << 411 G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
373 G4double dEmax = maxEProj - kinEnergyScatPro << 412
374 G4double dEmin = GetLowEnergyLimit(); << 413
375 G4double dE1 = dEmin; << 414 G4double dEmax=maxEProj-kinEnergyScatProj;
376 G4double dE2 = dEmin; << 415 G4double dEmin=GetLowEnergyLimit();
377 << 416 G4double dE1=dEmin;
378 std::vector<G4double>* log_ESec_vector = new << 417 G4double dE2=dEmin;
379 std::vector<G4double>* log_Prob_vector = new << 418
>> 419
>> 420 std::vector< double>* log_ESec_vector = new std::vector< double>();
>> 421 std::vector< double>* log_Prob_vector = new std::vector< double>();
380 log_ESec_vector->push_back(std::log(dEmin)); 422 log_ESec_vector->push_back(std::log(dEmin));
381 log_Prob_vector->push_back(-50.); 423 log_Prob_vector->push_back(-50.);
382 G4int nbins = std::max(int(std::log10(dEmax << 424 G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
383 G4double fE = std::pow(dEmax / dEmin, 1. / n << 425 G4double fE=std::pow(dEmax/dEmin,1./nbins);
384 << 426
385 G4double int_cross_section = 0.; << 427 G4double int_cross_section=0.;
386 // Loop checking, 07-Aug-2015, Vladimir Ivan << 428
387 while(dE1 < dEmax * 0.9999999999999) << 429 while (dE1 <dEmax*0.9999999999999){
388 { << 430 dE2=dE1*fE;
389 dE2 = dE1 * fE; << 431 int_cross_section +=integral.Simpson(this,
390 int_cross_section += << 432 &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
391 integral.Simpson(this, &G4VEmAdjointMode << 433 log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
392 minEProj + dE1, std::mi << 434 log_Prob_vector->push_back(std::log(int_cross_section));
393 log_ESec_vector->push_back(std::log(std::m << 435 dE1=dE2;
394 log_Prob_vector->push_back(std::log(int_cr << 436
395 dE1 = dE2; << 437 }
396 } << 438
397 << 439
398 std::vector<std::vector<G4double>*> res_mat; << 440
399 if(int_cross_section > 0.) << 441
400 { << 442
401 res_mat.push_back(log_ESec_vector); << 443 std::vector< std::vector<G4double> *> res_mat;
402 res_mat.push_back(log_Prob_vector); << 444 res_mat.clear();
403 } << 445 if (int_cross_section >0.) {
404 else { << 446 res_mat.push_back(log_ESec_vector);
405 delete log_ESec_vector; << 447 res_mat.push_back(log_Prob_vector);
406 delete log_Prob_vector; << 448 }
407 } << 449
408 <<
409 return res_mat; 450 return res_mat;
410 } 451 }
411 <<
412 ////////////////////////////////////////////// 452 //////////////////////////////////////////////////////////////////////////////
413 G4double G4VEmAdjointModel::SampleAdjSecEnergy << 453 //
414 std::size_t MatrixIndex, G4double aPrimEnerg << 454 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
415 { << 455 {
416 G4AdjointCSMatrix* theMatrix = (*fCSMatrixPr << 456
417 if(isScatProjToProj) << 457
418 theMatrix = (*fCSMatrixProjToProjBackScat) << 458 G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
419 std::vector<G4double>* theLogPrimEnergyVecto << 459 if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
420 theMatrix->GetLogPrimEnergyVector(); << 460 std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
421 << 461
422 if(theLogPrimEnergyVector->empty()) << 462 if (theLogPrimEnergyVector->size() ==0){
423 { << 463 G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
424 G4cout << "No data are contained in the gi << 464 G4cout<<"The sampling procedure will be stopped."<<G4endl;
425 G4cout << "The sampling procedure will be << 465 return 0.;
426 return 0.; << 466
427 } << 467 }
428 << 468
429 G4AdjointInterpolator* theInterpolator = G4A << 469 G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
430 G4double aLogPrimEnergy = std << 470 G4double aLogPrimEnergy = std::log(aPrimEnergy);
431 G4int ind = (G4int)theInterpolator->FindPosi << 471 size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
432 aLogPrimEnergy, *theLogPrimEnergyVector); << 472
433 << 473
434 G4double aLogPrimEnergy1, aLogPrimEnergy2; << 474 G4double aLogPrimEnergy1,aLogPrimEnergy2;
435 G4double aLogCS1, aLogCS2; << 475 G4double aLogCS1,aLogCS2;
436 G4double log01, log02; << 476 G4double log01,log02;
437 std::vector<G4double>* aLogSecondEnergyVecto << 477 std::vector< double>* aLogSecondEnergyVector1 =0;
438 std::vector<G4double>* aLogSecondEnergyVecto << 478 std::vector< double>* aLogSecondEnergyVector2 =0;
439 std::vector<G4double>* aLogProbVector1 << 479 std::vector< double>* aLogProbVector1=0;
440 std::vector<G4double>* aLogProbVector2 << 480 std::vector< double>* aLogProbVector2=0;
441 std::vector<std::size_t>* aLogProbVectorInde << 481 std::vector< size_t>* aLogProbVectorIndex1=0;
442 std::vector<std::size_t>* aLogProbVectorInde << 482 std::vector< size_t>* aLogProbVectorIndex2=0;
443 << 483
444 theMatrix->GetData(ind, aLogPrimEnergy1, aLo << 484 theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
445 aLogSecondEnergyVector1, << 485 theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
446 aLogProbVectorIndex1 ); << 486
447 theMatrix->GetData(ind + 1, aLogPrimEnergy2, << 487 G4double rand_var = G4UniformRand();
448 aLogSecondEnergyVector2, << 488 G4double log_rand_var= std::log(rand_var);
449 aLogProbVectorIndex2); << 489 G4double log_Tcut =std::log(currentTcutForDirectSecond);
450 << 490 G4double Esec=0;
451 if (! (aLogProbVector1 && aLogProbVector2 && << 491 G4double log_dE1,log_dE2;
452 aLogSecondEnergyVector1 && aLogSeco << 492 G4double log_rand_var1,log_rand_var2;
453 return 0.; << 493 G4double log_E1,log_E2;
454 } << 494 log_rand_var1=log_rand_var;
455 << 495 log_rand_var2=log_rand_var;
456 G4double rand_var = G4UniformRand(); << 496
457 G4double log_rand_var = std::log(rand_var); << 497 G4double Emin=0.;
458 G4double log_Tcut = std::log(fTcutSecond << 498 G4double Emax=0.;
459 G4double Esec = 0.; << 499 if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
460 G4double log_dE1, log_dE2; << 500 Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
461 G4double log_rand_var1, log_rand_var2; << 501 Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
462 G4double log_E1, log_E2; << 502 G4double dE=0;
463 log_rand_var1 = log_rand_var; << 503 if (Emin < Emax ){
464 log_rand_var2 = log_rand_var; << 504 if (ApplyCutInRange) {
465 << 505 if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
466 G4double Emin = 0.; << 506
467 G4double Emax = 0.; << 507 log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
468 if(theMatrix->IsScatProjToProj()) << 508 log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
469 { // case where Tcut plays a role << 509
470 Emin = GetSecondAdjEnergyMinForScatProjToP << 510 }
471 Emax = GetSecondAdjEnergyMaxForScatProjToP << 511 log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
472 G4double dE = 0.; << 512 log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
473 if(Emin < Emax) << 513 dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
474 { << 514 }
475 if(fApplyCutInRange) << 515
476 { << 516 Esec = aPrimEnergy +dE;
477 if(fSecondPartSameType && fTcutSecond << 517 Esec=std::max(Esec,Emin);
478 return aPrimEnergy; << 518 Esec=std::min(Esec,Emax);
479 << 519
480 log_rand_var1 = log_rand_var + << 520 }
481 theInterpolator->Inter << 521 else { //Tcut condition is already full-filled
482 log_Tcut, *aLogSecon << 522
483 log_rand_var2 = log_rand_var + << 523 log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
484 theInterpolator->Inter << 524 log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
485 log_Tcut, *aLogSecon << 525
486 } << 526 Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
487 log_dE1 = theInterpolator->Interpolate(l << 527 Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
488 * << 528 Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
489 log_dE2 = theInterpolator->Interpolate(l << 529 Esec=std::max(Esec,Emin);
490 * << 530 Esec=std::min(Esec,Emax);
491 dE = std::exp(theInterpolator->Line << 531
492 aLogPrimEnergy, aLogPrimEnergy1, aLogP <<
493 } <<
494 <<
495 Esec = aPrimEnergy + dE; <<
496 Esec = std::max(Esec, Emin); <<
497 Esec = std::min(Esec, Emax); <<
498 } <<
499 else <<
500 { // Tcut condition is already full-filled <<
501 <<
502 log_E1 = theInterpolator->Interpolate(log_ <<
503 *aLo <<
504 log_E2 = theInterpolator->Interpolate(log_ <<
505 *aLo <<
506 <<
507 Esec = std::exp(theInterpolator->LinearInt <<
508 aLogPrimEnergy, aLogPrimEnergy1, aLogPri <<
509 Emin = GetSecondAdjEnergyMinForProdToProj( <<
510 Emax = GetSecondAdjEnergyMaxForProdToProj( <<
511 Esec = std::max(Esec, Emin); <<
512 Esec = std::min(Esec, Emax); <<
513 } 532 }
>> 533
514 return Esec; 534 return Esec;
>> 535
>> 536
>> 537
>> 538
>> 539
515 } 540 }
516 541
517 ////////////////////////////////////////////// 542 //////////////////////////////////////////////////////////////////////////////
518 G4double G4VEmAdjointModel::SampleAdjSecEnergy << 543 //
519 G4double aPrimEnergy, G4bool isScatProjToPro << 544 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
520 { << 545 { SelectCSMatrix(IsScatProjToProjCase);
521 SelectCSMatrix(isScatProjToProj); << 546 return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
522 return SampleAdjSecEnergyFromCSMatrix(fCSMat <<
523 isScat <<
524 } 547 }
525 <<
526 ////////////////////////////////////////////// 548 //////////////////////////////////////////////////////////////////////////////
527 void G4VEmAdjointModel::SelectCSMatrix(G4bool << 549 //
528 { << 550 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
529 fCSMatrixUsed = 0; << 551 {
530 if(!fUseMatrixPerElement) << 552 indexOfUsedCrossSectionMatrix=0;
531 fCSMatrixUsed = fCurrentMaterial->GetIndex << 553 if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
532 else if(!fOneMatrixForAllElements) << 554 else if (!UseOnlyOneMatrixForAllElements) { //Select Material
533 { // Select Material << 555 std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
534 std::vector<G4double>* CS_Vs_Element = &fE << 556 lastCS=lastAdjointCSForScatProjToProjCase;
535 fLastCS = fLa << 557 if ( !IsScatProjToProjCase) {
536 if(!isScatProjToProj) << 558 CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
537 { << 559 lastCS=lastAdjointCSForProdToProjCase;
538 CS_Vs_Element = &fElementCSProdToProj; << 560 }
539 fLastCS = fLastAdjointCSForProdToP << 561 G4double rand_var= G4UniformRand();
540 } << 562 G4double SumCS=0.;
541 G4double SumCS = 0.; << 563 size_t ind=0;
542 std::size_t ind = 0; << 564 for (size_t i=0;i<CS_Vs_Element->size();i++){
543 for(std::size_t i = 0; i < CS_Vs_Element-> << 565 SumCS+=(*CS_Vs_Element)[i];
544 { << 566 if (rand_var<=SumCS/lastCS){
545 SumCS += (*CS_Vs_Element)[i]; << 567 ind=i;
546 if(G4UniformRand() <= SumCS / fLastCS) << 568 break;
547 { << 569 }
548 ind = i; << 570 }
549 break; << 571 indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
550 } <<
551 } <<
552 fCSMatrixUsed = fCurrentMaterial->GetEleme <<
553 } 572 }
554 } 573 }
555 <<
556 ////////////////////////////////////////////// 574 //////////////////////////////////////////////////////////////////////////////
557 G4double G4VEmAdjointModel::SampleAdjSecEnergy << 575 //
558 G4double prim_energy, G4bool isScatProjToPro << 576 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
559 { << 577 {
560 // here we try to use the rejection method << 578 // here we try to use the rejection method
561 constexpr G4int iimax = 1000; << 579 //-----------------------------------------
562 G4double E = 0.; << 580
563 G4double x, xmin, greject; << 581 G4double E=0;
564 if(isScatProjToProj) << 582 G4double x,xmin,greject,q;
565 { << 583 if ( IsScatProjToProjCase){
566 G4double Emax = GetSecondAdjEnergyMaxForSc << 584 G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
567 G4double Emin = prim_energy + fTcutSecond; << 585 G4double Emin= prim_energy+currentTcutForDirectSecond;
568 xmin = Emin / Emax; << 586 xmin=Emin/Emax;
569 G4double grejmax = << 587 G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy;
570 DiffCrossSectionPerAtomPrimToScatPrim(Em << 588
571 << 589 do {
572 G4int ii = 0; << 590 q = G4UniformRand();
573 do << 591 x = 1./(q*(1./xmin -1.) +1.);
574 { << 592 E=x*Emax;
575 // q = G4UniformRand(); << 593 greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy;
576 x = 1. / (G4UniformRand() * (1. / xmin - << 594
577 E = x * Emax; << 595 }
578 greject = << 596
579 DiffCrossSectionPerAtomPrimToScatPrim( << 597 while( greject < G4UniformRand()*grejmax );
580 ++ii; << 598
581 if(ii >= iimax) <<
582 { <<
583 break; <<
584 } <<
585 } <<
586 // Loop checking, 07-Aug-2015, Vladimir Iv <<
587 while(greject < G4UniformRand() * grejmax) <<
588 } <<
589 else <<
590 { <<
591 G4double Emax = GetSecondAdjEnergyMaxForPr <<
592 G4double Emin = GetSecondAdjEnergyMinForPr <<
593 xmin = Emin / Emax; <<
594 G4double grejmax = <<
595 DiffCrossSectionPerAtomPrimToSecond(Emin <<
596 G4int ii = 0; <<
597 do <<
598 { <<
599 x = std::pow(xmin, G4UniformRand() <<
600 E = x * Emax; <<
601 greject = DiffCrossSectionPerAtomPrimToS <<
602 ++ii; <<
603 if(ii >= iimax) <<
604 { <<
605 break; <<
606 } <<
607 } <<
608 // Loop checking, 07-Aug-2015, Vladimir Iv <<
609 while(greject < G4UniformRand() * grejmax) <<
610 } 599 }
611 << 600 else {
>> 601 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
>> 602 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
>> 603 xmin=Emin/Emax;
>> 604 G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1);
>> 605 do {
>> 606 q = G4UniformRand();
>> 607 x = std::pow(xmin, q);
>> 608 E=x*Emax;
>> 609 greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
>> 610
>> 611 }
>> 612
>> 613 while( greject < G4UniformRand()*grejmax );
>> 614
>> 615
>> 616
>> 617 }
>> 618
612 return E; 619 return E;
613 } 620 }
614 621
615 ////////////////////////////////////////////// 622 ////////////////////////////////////////////////////////////////////////////////
616 void G4VEmAdjointModel::CorrectPostStepWeight( << 623 //
617 << 624 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
618 << 625 G4double old_weight,
619 << 626 G4double adjointPrimKinEnergy,
620 << 627 G4double projectileKinEnergy,
621 { << 628 G4bool IsScatProjToProjCase)
622 G4double new_weight = old_weight; << 629 {
623 G4double w_corr = << 630 G4double new_weight=old_weight;
624 fCSManager->GetPostStepWeightCorrection() << 631 G4double w_corr =1./CS_biasing_factor;
625 << 632 w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
626 fLastCS = fLastAdjointCSForScatProjToProj; << 633
627 if(!isScatProjToProj) << 634
628 fLastCS = fLastAdjointCSForProdToProj; << 635 lastCS=lastAdjointCSForScatProjToProjCase;
629 if((adjointPrimKinEnergy - fPreStepEnergy) / << 636 if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
630 { << 637 if ((adjointPrimKinEnergy-preStepEnergy)/preStepEnergy>0.001){ //Is that in all cases needed???
631 G4double post_stepCS = AdjointCrossSection << 638 G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
632 fCurrentCouple, adjointPrimKinEnergy, is << 639 ,IsScatProjToProjCase );
633 if(post_stepCS > 0. && fLastCS > 0.) << 640 if (post_stepCS>0 && lastCS>0) w_corr*=post_stepCS/lastCS;
634 w_corr *= post_stepCS / fLastCS; << 641 }
635 } << 642
636 << 643 new_weight*=w_corr;
637 new_weight *= w_corr; << 644
638 new_weight *= projectileKinEnergy / adjointP << 645 //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
639 // This is needed due to the biasing of diff << 646 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
640 // by the factor adjointPrimKinEnergy/projec << 647 //by the factor adjointPrimKinEnergy/projectileKinEnergy
641 << 648
642 fParticleChange->SetParentWeightByProcess(fa << 649
643 fParticleChange->SetSecondaryWeightByProcess << 650
644 fParticleChange->ProposeParentWeight(new_wei << 651 fParticleChange->SetParentWeightByProcess(false);
>> 652 fParticleChange->SetSecondaryWeightByProcess(false);
>> 653 fParticleChange->ProposeParentWeight(new_weight);
645 } 654 }
646 <<
647 ////////////////////////////////////////////// 655 //////////////////////////////////////////////////////////////////////////////
648 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 656 //
649 G4double kinEnergyScatProj) << 657 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
650 { << 658 { G4double maxEProj= HighEnergyLimit;
651 G4double maxEProj = GetHighEnergyLimit(); << 659 if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
652 if(fSecondPartSameType) <<
653 maxEProj = std::min(kinEnergyScatProj * 2. <<
654 return maxEProj; 660 return maxEProj;
655 } 661 }
656 <<
657 ////////////////////////////////////////////// 662 //////////////////////////////////////////////////////////////////////////////
658 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 663 //
659 G4double primAdjEnergy, G4double tcut) << 664 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
660 { << 665 { G4double Emin=PrimAdjEnergy;
661 G4double Emin = primAdjEnergy; << 666 if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
662 if(fApplyCutInRange) <<
663 Emin += tcut; <<
664 return Emin; 667 return Emin;
665 } 668 }
666 <<
667 ////////////////////////////////////////////// 669 //////////////////////////////////////////////////////////////////////////////
668 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 670 //
669 { << 671 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
670 return fHighEnergyLimit; << 672 { return HighEnergyLimit;
671 } 673 }
672 <<
673 ////////////////////////////////////////////// 674 //////////////////////////////////////////////////////////////////////////////
674 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 675 //
675 G4double primAdjEnergy) << 676 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
676 { << 677 { G4double minEProj=PrimAdjEnergy;
677 G4double minEProj = primAdjEnergy; << 678 if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
678 if(fSecondPartSameType) <<
679 minEProj = primAdjEnergy * 2.; <<
680 return minEProj; 679 return minEProj;
681 } 680 }
682 <<
683 ////////////////////////////////////////////// 681 ////////////////////////////////////////////////////////////////////////////////////////////
684 void G4VEmAdjointModel::DefineCurrentMaterial( << 682 //
685 const G4MaterialCutsCouple* couple) << 683 void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
686 { << 684 { if(couple != currentCouple) {
687 if(couple != fCurrentCouple) << 685 currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
688 { << 686 currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
689 fCurrentCouple = const_cast<G4MaterialCu << 687 currentCoupleIndex = couple->GetIndex();
690 fCurrentMaterial = const_cast<G4Material*> << 688 currentMaterialIndex = currentMaterial->GetIndex();
691 std::size_t idx = 56; << 689 size_t idx=56;
692 fTcutSecond = 1.e-11; << 690 currentTcutForDirectSecond =0.00000000001;
693 if(fAdjEquivDirectSecondPart) << 691 if (theAdjEquivOfDirectSecondPartDef) {
694 { << 692 if (theAdjEquivOfDirectSecondPartDef == G4AdjointGamma::AdjointGamma()) idx = 0;
695 if(fAdjEquivDirectSecondPart == G4Adjoin << 693 else if (theAdjEquivOfDirectSecondPartDef == G4AdjointElectron::AdjointElectron()) idx = 1;
696 idx = 0; << 694 else if (theAdjEquivOfDirectSecondPartDef == G4AdjointPositron::AdjointPositron()) idx = 2;
697 else if(fAdjEquivDirectSecondPart == G4A << 695 if (idx <56){
698 idx = 1; << 696 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
699 else if(fAdjEquivDirectSecondPart == G4A << 697 currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
700 idx = 2; << 698 }
701 if(idx < 56) << 699 }
702 { << 700
703 const std::vector<G4double>* aVec = << 701
704 G4ProductionCutsTable::GetProduction <<
705 idx); <<
706 fTcutSecond = (*aVec)[couple->GetIndex <<
707 } <<
708 } <<
709 } 702 }
710 } 703 }
711 <<
712 ////////////////////////////////////////////// 704 ////////////////////////////////////////////////////////////////////////////////////////////
>> 705 //
713 void G4VEmAdjointModel::SetHighEnergyLimit(G4d 706 void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
714 { << 707 { HighEnergyLimit=aVal;
715 fHighEnergyLimit = aVal; << 708 if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
716 if(fDirectModel) <<
717 fDirectModel->SetHighEnergyLimit(aVal); <<
718 } 709 }
719 <<
720 ////////////////////////////////////////////// 710 ////////////////////////////////////////////////////////////////////////////////////////////
>> 711 //
721 void G4VEmAdjointModel::SetLowEnergyLimit(G4do 712 void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
722 { 713 {
723 fLowEnergyLimit = aVal; << 714 LowEnergyLimit=aVal;
724 if(fDirectModel) << 715 if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
725 fDirectModel->SetLowEnergyLimit(aVal); <<
726 } 716 }
727 <<
728 ////////////////////////////////////////////// 717 ////////////////////////////////////////////////////////////////////////////////////////////
729 void G4VEmAdjointModel::SetAdjointEquivalentOf << 718 //
730 G4ParticleDefinition* aPart) << 719 void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
731 { 720 {
732 fAdjEquivDirectPrimPart = aPart; << 721 theAdjEquivOfDirectPrimPartDef=aPart;
733 if(fAdjEquivDirectPrimPart->GetParticleName( << 722 if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
734 fDirectPrimaryPart = G4Electron::Electron( << 723 theDirectPrimaryPartDef=G4Electron::Electron();
735 else if(fAdjEquivDirectPrimPart->GetParticle << 724 if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
736 fDirectPrimaryPart = G4Gamma::Gamma(); << 725 theDirectPrimaryPartDef=G4Gamma::Gamma();
737 } 726 }
738 727