| 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 93358 2015-10-19 13:41:21Z gcosmo $
>> 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); <<
224 <<
225 G4double int_cross_section = 0.; <<
226 // Loop checking, 07-Aug-2015, Vladimir Ivan 268 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
227 while(E1 < maxEProj * 0.9999999) << 269 while (E1 <maxEProj*0.9999999){
228 { << 270 //G4cout<<E1<<'\t'<<E2<<G4endl;
229 int_cross_section += << 271
230 integral.Simpson(this, &G4VEmAdjointMode << 272 int_cross_section +=integral.Simpson(this,
231 std::min(E2, maxEProj * << 273 &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
232 log_ESec_vector->push_back(std::log(std::m << 274 log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
233 log_Prob_vector->push_back(std::log(int_cr << 275 log_Prob_vector->push_back(std::log(int_cross_section));
234 E1 = E2; << 276 E1=E2;
235 E2 *= fE; << 277 E2*=fE;
236 } << 278
237 std::vector<std::vector<G4double>*> res_mat; << 279 }
238 if(int_cross_section > 0.) << 280 std::vector< std::vector<G4double>* > res_mat;
239 { << 281 res_mat.clear();
240 res_mat.push_back(log_ESec_vector); << 282 if (int_cross_section >0.) {
241 res_mat.push_back(log_Prob_vector); << 283 res_mat.push_back(log_ESec_vector);
242 } << 284 res_mat.push_back(log_Prob_vector);
243 else { << 285 }
244 delete log_ESec_vector; << 286
245 delete log_Prob_vector; <<
246 } <<
247 return res_mat; 287 return res_mat;
248 } 288 }
249 289
250 ////////////////////////////////////////////// 290 /////////////////////////////////////////////////////////////////////////////////////
251 std::vector<std::vector<G4double>*> << 291 //
252 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 292 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerAtomForScatProj(
253 G4double kinEnergyScatProj, G4double Z, G4do << 293 G4double kinEnergyScatProj,
254 G4int nbin_pro_decade) // nb bins pro order << 294 G4double Z,
255 { << 295 G4double A ,
256 G4Integrator<G4VEmAdjointModel, G4double (G4 << 296 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
257 integral; << 297 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
258 fASelectedNucleus = G4lrint(A << 298 ASelectedNucleus=int(A);
259 fZSelectedNucleus = G4lrint(Z << 299 ZSelectedNucleus=int(Z);
260 fKinEnergyScatProjForIntegration = kinEnergy << 300 kinEnergyScatProjForIntegration = kinEnergyScatProj;
261 << 301
262 // compute the vector of integrated cross se << 302 //compute the vector of integrated cross sections
263 G4double minEProj = GetSecondAdjEnergyMinFor << 303 //-------------------
264 G4double maxEProj = GetSecondAdjEnergyMaxFor << 304
265 G4double dEmax = maxEProj - kinEnergyScat << 305 G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
266 G4double dEmin = GetLowEnergyLimit(); << 306 G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
267 G4double dE1 = dEmin; << 307 G4double dEmax=maxEProj-kinEnergyScatProj;
268 G4double dE2 = dEmin; << 308 G4double dEmin=GetLowEnergyLimit();
269 << 309 G4double dE1=dEmin;
270 std::vector<G4double>* log_ESec_vector = new << 310 G4double dE2=dEmin;
271 std::vector<G4double>* log_Prob_vector = new << 311
>> 312
>> 313 std::vector< double>* log_ESec_vector = new std::vector< double>();
>> 314 std::vector< double>* log_Prob_vector = new std::vector< double>();
272 log_ESec_vector->push_back(std::log(dEmin)); 315 log_ESec_vector->push_back(std::log(dEmin));
273 log_Prob_vector->push_back(-50.); 316 log_Prob_vector->push_back(-50.);
274 G4int nbins = std::max(G4int(std::log10(dEma << 317 G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
275 G4double fE = std::pow(dEmax / dEmin, 1. / n << 318 G4double fE=std::pow(dEmax/dEmin,1./nbins);
276 << 319
277 G4double int_cross_section = 0.; << 320
>> 321
>> 322
>> 323
>> 324 G4double int_cross_section=0.;
>> 325
278 // Loop checking, 07-Aug-2015, Vladimir Ivan 326 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
279 while(dE1 < dEmax * 0.9999999999999) << 327 while (dE1 <dEmax*0.9999999999999){
280 { << 328 dE2=dE1*fE;
281 dE2 = dE1 * fE; << 329 int_cross_section +=integral.Simpson(this,
282 int_cross_section += << 330 &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
283 integral.Simpson(this, &G4VEmAdjointMode << 331 //G4cout<<"int_cross_section "<<minEProj+dE1<<'\t'<<int_cross_section<<G4endl;
284 minEProj + dE1, std::mi << 332 log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
285 log_ESec_vector->push_back(std::log(std::m << 333 log_Prob_vector->push_back(std::log(int_cross_section));
286 log_Prob_vector->push_back(std::log(int_cr << 334 dE1=dE2;
287 dE1 = dE2; << 335
288 } << 336 }
289 << 337
290 std::vector<std::vector<G4double>*> res_mat; << 338
291 if(int_cross_section > 0.) << 339 std::vector< std::vector<G4double> *> res_mat;
292 { << 340 res_mat.clear();
293 res_mat.push_back(log_ESec_vector); << 341 if (int_cross_section >0.) {
294 res_mat.push_back(log_Prob_vector); << 342 res_mat.push_back(log_ESec_vector);
295 } << 343 res_mat.push_back(log_Prob_vector);
296 else { << 344 }
297 delete log_ESec_vector; << 345
298 delete log_Prob_vector; <<
299 } <<
300 <<
301 return res_mat; 346 return res_mat;
302 } 347 }
303 <<
304 ////////////////////////////////////////////// 348 ////////////////////////////////////////////////////////////////////////////////
305 std::vector<std::vector<G4double>*> << 349 //
306 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 350 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForSecond(
307 G4Material* aMaterial, G4double kinEnergyPro << 351 G4Material* aMaterial,
308 G4int nbin_pro_decade) // nb bins pro order << 352 G4double kinEnergyProd,
309 { << 353 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
310 G4Integrator<G4VEmAdjointModel, G4double (G4 << 354 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
311 integral; << 355 SelectedMaterial= aMaterial;
312 fSelectedMaterial = aMaterial; << 356 kinEnergyProdForIntegration = kinEnergyProd;
313 fKinEnergyProdForIntegration = kinEnergyProd << 357 //compute the vector of integrated cross sections
314 << 358 //-------------------
315 // compute the vector of integrated cross se << 359
316 G4double minEProj = GetSecondAdjEnergyMinFor << 360 G4double minEProj= GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd);
317 G4double maxEProj = GetSecondAdjEnergyMaxFor << 361 G4double maxEProj= GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd);
318 G4double E1 = minEProj; << 362 G4double E1=minEProj;
319 std::vector<G4double>* log_ESec_vector = new << 363 std::vector< double>* log_ESec_vector = new std::vector< double>();
320 std::vector<G4double>* log_Prob_vector = new << 364 std::vector< double>* log_Prob_vector = new std::vector< double>();
>> 365 log_ESec_vector->clear();
>> 366 log_Prob_vector->clear();
321 log_ESec_vector->push_back(std::log(E1)); 367 log_ESec_vector->push_back(std::log(E1));
322 log_Prob_vector->push_back(-50.); 368 log_Prob_vector->push_back(-50.);
323 << 369
324 G4double E2 = << 370 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 << 371 G4double fE=std::pow(10.,1./nbin_pro_decade);
326 nbin_pro_decade); << 372 G4double int_cross_section=0.;
327 G4double fE = std::pow(10., 1. / nbin_pro_de << 373
328 << 374 if (std::pow(fE,5.)>(maxEProj/minEProj)) fE = std::pow(maxEProj/minEProj,0.2);
329 if(std::pow(fE, 5.) > (maxEProj / minEProj)) << 375
330 fE = std::pow(maxEProj / minEProj, 0.2); <<
331 <<
332 G4double int_cross_section = 0.; <<
333 // Loop checking, 07-Aug-2015, Vladimir Ivan 376 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
334 while(E1 < maxEProj * 0.9999999) << 377 while (E1 <maxEProj*0.9999999){
335 { << 378
336 int_cross_section += << 379 int_cross_section +=integral.Simpson(this,
337 integral.Simpson(this, &G4VEmAdjointMode << 380 &G4VEmAdjointModel::DiffCrossSectionFunction1,E1,std::min(E2,maxEProj*0.99999999), 5);
338 std::min(E2, maxEProj * << 381 log_ESec_vector->push_back(std::log(std::min(E2,maxEProj)));
339 log_ESec_vector->push_back(std::log(std::m << 382 log_Prob_vector->push_back(std::log(int_cross_section));
340 log_Prob_vector->push_back(std::log(int_cr << 383 E1=E2;
341 E1 = E2; << 384 E2*=fE;
342 E2 *= fE; << 385
343 } << 386 }
344 std::vector<std::vector<G4double>*> res_mat; << 387 std::vector< std::vector<G4double>* > res_mat;
345 << 388 res_mat.clear();
346 if(int_cross_section > 0.) << 389
347 { << 390 if (int_cross_section >0.) {
348 res_mat.push_back(log_ESec_vector); << 391 res_mat.push_back(log_ESec_vector);
349 res_mat.push_back(log_Prob_vector); << 392 res_mat.push_back(log_Prob_vector);
350 } << 393 }
351 else { << 394
352 delete log_ESec_vector; << 395
353 delete log_Prob_vector; << 396
354 } <<
355 return res_mat; 397 return res_mat;
356 } 398 }
357 399
358 ////////////////////////////////////////////// << 400 /////////////////////////////////////////////////////////////////////////////////////
359 std::vector<std::vector<G4double>*> << 401 //
360 G4VEmAdjointModel::ComputeAdjointCrossSectionV << 402 std::vector< std::vector<G4double>* > G4VEmAdjointModel::ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
361 G4Material* aMaterial, G4double kinEnergySca << 403 G4Material* aMaterial,
362 G4int nbin_pro_decade) // nb bins pro order << 404 G4double kinEnergyScatProj,
363 { << 405 G4int nbin_pro_decade) //nb bins pro order of magnitude of energy
364 G4Integrator<G4VEmAdjointModel, G4double (G4 << 406 { G4Integrator<G4VEmAdjointModel, double(G4VEmAdjointModel::*)(double)> integral;
365 integral; << 407 SelectedMaterial= aMaterial;
366 fSelectedMaterial = aMaterial << 408 kinEnergyScatProjForIntegration = kinEnergyScatProj;
367 fKinEnergyScatProjForIntegration = kinEnergy << 409
368 << 410 //compute the vector of integrated cross sections
369 // compute the vector of integrated cross se << 411 //-------------------
370 G4double minEProj = GetSecondAdjEnergyMinFor << 412
371 G4double maxEProj = GetSecondAdjEnergyMaxFor << 413 G4double minEProj= GetSecondAdjEnergyMinForScatProjToProjCase(kinEnergyScatProj);
372 << 414 G4double maxEProj= GetSecondAdjEnergyMaxForScatProjToProjCase(kinEnergyScatProj);
373 G4double dEmax = maxEProj - kinEnergyScatPro << 415
374 G4double dEmin = GetLowEnergyLimit(); << 416
375 G4double dE1 = dEmin; << 417 G4double dEmax=maxEProj-kinEnergyScatProj;
376 G4double dE2 = dEmin; << 418 G4double dEmin=GetLowEnergyLimit();
377 << 419 G4double dE1=dEmin;
378 std::vector<G4double>* log_ESec_vector = new << 420 G4double dE2=dEmin;
379 std::vector<G4double>* log_Prob_vector = new << 421
>> 422
>> 423 std::vector< double>* log_ESec_vector = new std::vector< double>();
>> 424 std::vector< double>* log_Prob_vector = new std::vector< double>();
380 log_ESec_vector->push_back(std::log(dEmin)); 425 log_ESec_vector->push_back(std::log(dEmin));
381 log_Prob_vector->push_back(-50.); 426 log_Prob_vector->push_back(-50.);
382 G4int nbins = std::max(int(std::log10(dEmax << 427 G4int nbins=std::max( int(std::log10(dEmax/dEmin))*nbin_pro_decade,5);
383 G4double fE = std::pow(dEmax / dEmin, 1. / n << 428 G4double fE=std::pow(dEmax/dEmin,1./nbins);
384 << 429
385 G4double int_cross_section = 0.; << 430 G4double int_cross_section=0.;
>> 431
386 // Loop checking, 07-Aug-2015, Vladimir Ivan 432 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
387 while(dE1 < dEmax * 0.9999999999999) << 433 while (dE1 <dEmax*0.9999999999999){
388 { << 434 dE2=dE1*fE;
389 dE2 = dE1 * fE; << 435 int_cross_section +=integral.Simpson(this,
390 int_cross_section += << 436 &G4VEmAdjointModel::DiffCrossSectionFunction2,minEProj+dE1,std::min(minEProj+dE2,maxEProj), 5);
391 integral.Simpson(this, &G4VEmAdjointMode << 437 log_ESec_vector->push_back(std::log(std::min(dE2,maxEProj-minEProj)));
392 minEProj + dE1, std::mi << 438 log_Prob_vector->push_back(std::log(int_cross_section));
393 log_ESec_vector->push_back(std::log(std::m << 439 dE1=dE2;
394 log_Prob_vector->push_back(std::log(int_cr << 440
395 dE1 = dE2; << 441 }
396 } << 442
397 << 443
398 std::vector<std::vector<G4double>*> res_mat; << 444
399 if(int_cross_section > 0.) << 445
400 { << 446
401 res_mat.push_back(log_ESec_vector); << 447 std::vector< std::vector<G4double> *> res_mat;
402 res_mat.push_back(log_Prob_vector); << 448 res_mat.clear();
403 } << 449 if (int_cross_section >0.) {
404 else { << 450 res_mat.push_back(log_ESec_vector);
405 delete log_ESec_vector; << 451 res_mat.push_back(log_Prob_vector);
406 delete log_Prob_vector; << 452 }
407 } << 453
408 <<
409 return res_mat; 454 return res_mat;
410 } 455 }
411 <<
412 ////////////////////////////////////////////// 456 //////////////////////////////////////////////////////////////////////////////
413 G4double G4VEmAdjointModel::SampleAdjSecEnergy << 457 //
414 std::size_t MatrixIndex, G4double aPrimEnerg << 458 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(size_t MatrixIndex,G4double aPrimEnergy,G4bool IsScatProjToProjCase)
415 { << 459 {
416 G4AdjointCSMatrix* theMatrix = (*fCSMatrixPr << 460
417 if(isScatProjToProj) << 461
418 theMatrix = (*fCSMatrixProjToProjBackScat) << 462 G4AdjointCSMatrix* theMatrix= (*pOnCSMatrixForProdToProjBackwardScattering)[MatrixIndex];
419 std::vector<G4double>* theLogPrimEnergyVecto << 463 if (IsScatProjToProjCase) theMatrix= (*pOnCSMatrixForScatProjToProjBackwardScattering)[MatrixIndex];
420 theMatrix->GetLogPrimEnergyVector(); << 464 std::vector< double>* theLogPrimEnergyVector = theMatrix->GetLogPrimEnergyVector();
421 << 465
422 if(theLogPrimEnergyVector->empty()) << 466 if (theLogPrimEnergyVector->size() ==0){
423 { << 467 G4cout<<"No data are contained in the given AdjointCSMatrix!"<<G4endl;
424 G4cout << "No data are contained in the gi << 468 G4cout<<"The sampling procedure will be stopped."<<G4endl;
425 G4cout << "The sampling procedure will be << 469 return 0.;
426 return 0.; << 470
427 } << 471 }
428 << 472
429 G4AdjointInterpolator* theInterpolator = G4A << 473 G4AdjointInterpolator* theInterpolator=G4AdjointInterpolator::GetInstance();
430 G4double aLogPrimEnergy = std << 474 G4double aLogPrimEnergy = std::log(aPrimEnergy);
431 G4int ind = (G4int)theInterpolator->FindPosi << 475 size_t ind =theInterpolator->FindPositionForLogVector(aLogPrimEnergy,*theLogPrimEnergyVector);
432 aLogPrimEnergy, *theLogPrimEnergyVector); << 476
433 << 477
434 G4double aLogPrimEnergy1, aLogPrimEnergy2; << 478 G4double aLogPrimEnergy1,aLogPrimEnergy2;
435 G4double aLogCS1, aLogCS2; << 479 G4double aLogCS1,aLogCS2;
436 G4double log01, log02; << 480 G4double log01,log02;
437 std::vector<G4double>* aLogSecondEnergyVecto << 481 std::vector< double>* aLogSecondEnergyVector1 =0;
438 std::vector<G4double>* aLogSecondEnergyVecto << 482 std::vector< double>* aLogSecondEnergyVector2 =0;
439 std::vector<G4double>* aLogProbVector1 << 483 std::vector< double>* aLogProbVector1=0;
440 std::vector<G4double>* aLogProbVector2 << 484 std::vector< double>* aLogProbVector2=0;
441 std::vector<std::size_t>* aLogProbVectorInde << 485 std::vector< size_t>* aLogProbVectorIndex1=0;
442 std::vector<std::size_t>* aLogProbVectorInde << 486 std::vector< size_t>* aLogProbVectorIndex2=0;
443 << 487
444 theMatrix->GetData(ind, aLogPrimEnergy1, aLo << 488 theMatrix->GetData(ind, aLogPrimEnergy1,aLogCS1,log01, aLogSecondEnergyVector1,aLogProbVector1,aLogProbVectorIndex1);
445 aLogSecondEnergyVector1, << 489 theMatrix->GetData(ind+1, aLogPrimEnergy2,aLogCS2,log02, aLogSecondEnergyVector2,aLogProbVector2,aLogProbVectorIndex2);
446 aLogProbVectorIndex1 ); << 490
447 theMatrix->GetData(ind + 1, aLogPrimEnergy2, << 491 G4double rand_var = G4UniformRand();
448 aLogSecondEnergyVector2, << 492 G4double log_rand_var= std::log(rand_var);
449 aLogProbVectorIndex2); << 493 G4double log_Tcut =std::log(currentTcutForDirectSecond);
450 << 494 G4double Esec=0;
451 if (! (aLogProbVector1 && aLogProbVector2 && << 495 G4double log_dE1,log_dE2;
452 aLogSecondEnergyVector1 && aLogSeco << 496 G4double log_rand_var1,log_rand_var2;
453 return 0.; << 497 G4double log_E1,log_E2;
454 } << 498 log_rand_var1=log_rand_var;
455 << 499 log_rand_var2=log_rand_var;
456 G4double rand_var = G4UniformRand(); << 500
457 G4double log_rand_var = std::log(rand_var); << 501 G4double Emin=0.;
458 G4double log_Tcut = std::log(fTcutSecond << 502 G4double Emax=0.;
459 G4double Esec = 0.; << 503 if (theMatrix->IsScatProjToProjCase()){ //case where Tcut plays a role
460 G4double log_dE1, log_dE2; << 504 Emin=GetSecondAdjEnergyMinForScatProjToProjCase(aPrimEnergy,currentTcutForDirectSecond);
461 G4double log_rand_var1, log_rand_var2; << 505 Emax=GetSecondAdjEnergyMaxForScatProjToProjCase(aPrimEnergy);
462 G4double log_E1, log_E2; << 506 G4double dE=0;
463 log_rand_var1 = log_rand_var; << 507 if (Emin < Emax ){
464 log_rand_var2 = log_rand_var; << 508 if (ApplyCutInRange) {
465 << 509 if (second_part_of_same_type && currentTcutForDirectSecond>aPrimEnergy) return aPrimEnergy;
466 G4double Emin = 0.; << 510
467 G4double Emax = 0.; << 511 log_rand_var1=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector1,*aLogProbVector1);
468 if(theMatrix->IsScatProjToProj()) << 512 log_rand_var2=log_rand_var+theInterpolator->InterpolateForLogVector(log_Tcut,*aLogSecondEnergyVector2,*aLogProbVector2);
469 { // case where Tcut plays a role << 513
470 Emin = GetSecondAdjEnergyMinForScatProjToP << 514 }
471 Emax = GetSecondAdjEnergyMaxForScatProjToP << 515 log_dE1 = theInterpolator->Interpolate(log_rand_var1,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
472 G4double dE = 0.; << 516 log_dE2 = theInterpolator->Interpolate(log_rand_var2,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
473 if(Emin < Emax) << 517 dE=std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_dE1,log_dE2));
474 { << 518 }
475 if(fApplyCutInRange) << 519
476 { << 520 Esec = aPrimEnergy +dE;
477 if(fSecondPartSameType && fTcutSecond << 521 Esec=std::max(Esec,Emin);
478 return aPrimEnergy; << 522 Esec=std::min(Esec,Emax);
479 << 523
480 log_rand_var1 = log_rand_var + << 524 }
481 theInterpolator->Inter << 525 else { //Tcut condition is already full-filled
482 log_Tcut, *aLogSecon << 526
483 log_rand_var2 = log_rand_var + << 527 log_E1 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector1,*aLogSecondEnergyVector1,"Lin");
484 theInterpolator->Inter << 528 log_E2 = theInterpolator->Interpolate(log_rand_var,*aLogProbVector2,*aLogSecondEnergyVector2,"Lin");
485 log_Tcut, *aLogSecon << 529
486 } << 530 Esec = std::exp(theInterpolator->LinearInterpolation(aLogPrimEnergy,aLogPrimEnergy1,aLogPrimEnergy2,log_E1,log_E2));
487 log_dE1 = theInterpolator->Interpolate(l << 531 Emin=GetSecondAdjEnergyMinForProdToProjCase(aPrimEnergy);
488 * << 532 Emax=GetSecondAdjEnergyMaxForProdToProjCase(aPrimEnergy);
489 log_dE2 = theInterpolator->Interpolate(l << 533 Esec=std::max(Esec,Emin);
490 * << 534 Esec=std::min(Esec,Emax);
491 dE = std::exp(theInterpolator->Line << 535
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 } 536 }
>> 537
514 return Esec; 538 return Esec;
>> 539
>> 540
>> 541
>> 542
>> 543
515 } 544 }
516 545
517 ////////////////////////////////////////////// 546 //////////////////////////////////////////////////////////////////////////////
518 G4double G4VEmAdjointModel::SampleAdjSecEnergy << 547 //
519 G4double aPrimEnergy, G4bool isScatProjToPro << 548 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromCSMatrix(G4double aPrimEnergy,G4bool IsScatProjToProjCase)
520 { << 549 { SelectCSMatrix(IsScatProjToProjCase);
521 SelectCSMatrix(isScatProjToProj); << 550 return SampleAdjSecEnergyFromCSMatrix(indexOfUsedCrossSectionMatrix, aPrimEnergy, IsScatProjToProjCase);
522 return SampleAdjSecEnergyFromCSMatrix(fCSMat <<
523 isScat <<
524 } 551 }
525 <<
526 ////////////////////////////////////////////// 552 //////////////////////////////////////////////////////////////////////////////
527 void G4VEmAdjointModel::SelectCSMatrix(G4bool << 553 //
528 { << 554 void G4VEmAdjointModel::SelectCSMatrix(G4bool IsScatProjToProjCase)
529 fCSMatrixUsed = 0; << 555 {
530 if(!fUseMatrixPerElement) << 556 indexOfUsedCrossSectionMatrix=0;
531 fCSMatrixUsed = fCurrentMaterial->GetIndex << 557 if (!UseMatrixPerElement) indexOfUsedCrossSectionMatrix = currentMaterialIndex;
532 else if(!fOneMatrixForAllElements) << 558 else if (!UseOnlyOneMatrixForAllElements) { //Select Material
533 { // Select Material << 559 std::vector<G4double>* CS_Vs_Element = &CS_Vs_ElementForScatProjToProjCase;
534 std::vector<G4double>* CS_Vs_Element = &fE << 560 lastCS=lastAdjointCSForScatProjToProjCase;
535 fLastCS = fLa << 561 if ( !IsScatProjToProjCase) {
536 if(!isScatProjToProj) << 562 CS_Vs_Element = &CS_Vs_ElementForProdToProjCase;
537 { << 563 lastCS=lastAdjointCSForProdToProjCase;
538 CS_Vs_Element = &fElementCSProdToProj; << 564 }
539 fLastCS = fLastAdjointCSForProdToP << 565 G4double rand_var= G4UniformRand();
540 } << 566 G4double SumCS=0.;
541 G4double SumCS = 0.; << 567 size_t ind=0;
542 std::size_t ind = 0; << 568 for (size_t i=0;i<CS_Vs_Element->size();i++){
543 for(std::size_t i = 0; i < CS_Vs_Element-> << 569 SumCS+=(*CS_Vs_Element)[i];
544 { << 570 if (rand_var<=SumCS/lastCS){
545 SumCS += (*CS_Vs_Element)[i]; << 571 ind=i;
546 if(G4UniformRand() <= SumCS / fLastCS) << 572 break;
547 { << 573 }
548 ind = i; << 574 }
549 break; << 575 indexOfUsedCrossSectionMatrix = currentMaterial->GetElement(ind)->GetIndex();
550 } <<
551 } <<
552 fCSMatrixUsed = fCurrentMaterial->GetEleme <<
553 } 576 }
554 } 577 }
555 <<
556 ////////////////////////////////////////////// 578 //////////////////////////////////////////////////////////////////////////////
557 G4double G4VEmAdjointModel::SampleAdjSecEnergy << 579 //
558 G4double prim_energy, G4bool isScatProjToPro << 580 G4double G4VEmAdjointModel::SampleAdjSecEnergyFromDiffCrossSectionPerAtom(G4double prim_energy,G4bool IsScatProjToProjCase)
559 { << 581 {
560 // here we try to use the rejection method << 582 // here we try to use the rejection method
561 constexpr G4int iimax = 1000; << 583 //-----------------------------------------
562 G4double E = 0.; << 584
563 G4double x, xmin, greject; << 585 const G4int iimax = 1000;
564 if(isScatProjToProj) << 586 G4double E=0;
565 { << 587 G4double x,xmin,greject,q;
566 G4double Emax = GetSecondAdjEnergyMaxForSc << 588 if ( IsScatProjToProjCase){
567 G4double Emin = prim_energy + fTcutSecond; << 589 G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(prim_energy);
568 xmin = Emin / Emax; << 590 G4double Emin= prim_energy+currentTcutForDirectSecond;
569 G4double grejmax = << 591 xmin=Emin/Emax;
570 DiffCrossSectionPerAtomPrimToScatPrim(Em << 592 G4double grejmax = DiffCrossSectionPerAtomPrimToScatPrim(Emin,prim_energy,1)*prim_energy;
571 << 593
572 G4int ii = 0; << 594 G4int ii =0;
573 do << 595 do {
574 { << 596 q = G4UniformRand();
575 // q = G4UniformRand(); << 597 x = 1./(q*(1./xmin -1.) +1.);
576 x = 1. / (G4UniformRand() * (1. / xmin - << 598 E=x*Emax;
577 E = x * Emax; << 599 greject = DiffCrossSectionPerAtomPrimToScatPrim( E,prim_energy ,1)*prim_energy;
578 greject = << 600 ++ii;
579 DiffCrossSectionPerAtomPrimToScatPrim( << 601 if(ii >= iimax) { break; }
580 ++ii; << 602 }
581 if(ii >= iimax) << 603 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
582 { << 604 while( greject < G4UniformRand()*grejmax );
583 break; << 605
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 } 606 }
611 << 607 else {
>> 608 G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(prim_energy);
>> 609 G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(prim_energy);;
>> 610 xmin=Emin/Emax;
>> 611 G4double grejmax = DiffCrossSectionPerAtomPrimToSecond(Emin,prim_energy,1);
>> 612 G4int ii =0;
>> 613 do {
>> 614 q = G4UniformRand();
>> 615 x = std::pow(xmin, q);
>> 616 E=x*Emax;
>> 617 greject = DiffCrossSectionPerAtomPrimToSecond( E,prim_energy ,1);
>> 618 ++ii;
>> 619 if(ii >= iimax) { break; }
>> 620 }
>> 621 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
>> 622 while( greject < G4UniformRand()*grejmax );
>> 623
>> 624
>> 625
>> 626 }
>> 627
612 return E; 628 return E;
613 } 629 }
614 630
615 ////////////////////////////////////////////// 631 ////////////////////////////////////////////////////////////////////////////////
616 void G4VEmAdjointModel::CorrectPostStepWeight( << 632 //
617 << 633 void G4VEmAdjointModel::CorrectPostStepWeight(G4ParticleChange* fParticleChange,
618 << 634 G4double old_weight,
619 << 635 G4double adjointPrimKinEnergy,
620 << 636 G4double projectileKinEnergy,
621 { << 637 G4bool IsScatProjToProjCase)
622 G4double new_weight = old_weight; << 638 {
623 G4double w_corr = << 639 G4double new_weight=old_weight;
624 fCSManager->GetPostStepWeightCorrection() << 640 G4double w_corr =1./CS_biasing_factor;
625 << 641 w_corr*=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection();
626 fLastCS = fLastAdjointCSForScatProjToProj; << 642
627 if(!isScatProjToProj) << 643
628 fLastCS = fLastAdjointCSForProdToProj; << 644 lastCS=lastAdjointCSForScatProjToProjCase;
629 if((adjointPrimKinEnergy - fPreStepEnergy) / << 645 if ( !IsScatProjToProjCase) lastCS=lastAdjointCSForProdToProjCase;
630 { << 646 if ((adjointPrimKinEnergy-preStepEnergy)/preStepEnergy>0.001){ //Is that in all cases needed???
631 G4double post_stepCS = AdjointCrossSection << 647 G4double post_stepCS=AdjointCrossSection(currentCouple, adjointPrimKinEnergy
632 fCurrentCouple, adjointPrimKinEnergy, is << 648 ,IsScatProjToProjCase );
633 if(post_stepCS > 0. && fLastCS > 0.) << 649 if (post_stepCS>0 && lastCS>0) w_corr*=post_stepCS/lastCS;
634 w_corr *= post_stepCS / fLastCS; << 650 }
635 } << 651
636 << 652 new_weight*=w_corr;
637 new_weight *= w_corr; << 653
638 new_weight *= projectileKinEnergy / adjointP << 654 //G4cout<<"Post step "<<new_weight<<'\t'<<w_corr<<'\t'<<old_weight<<G4endl;
639 // This is needed due to the biasing of diff << 655 new_weight*=projectileKinEnergy/adjointPrimKinEnergy;//This is needed due to the biasing of diff CS
640 // by the factor adjointPrimKinEnergy/projec << 656 //by the factor adjointPrimKinEnergy/projectileKinEnergy
641 << 657
642 fParticleChange->SetParentWeightByProcess(fa << 658
643 fParticleChange->SetSecondaryWeightByProcess << 659
644 fParticleChange->ProposeParentWeight(new_wei << 660 fParticleChange->SetParentWeightByProcess(false);
>> 661 fParticleChange->SetSecondaryWeightByProcess(false);
>> 662 fParticleChange->ProposeParentWeight(new_weight);
645 } 663 }
646 <<
647 ////////////////////////////////////////////// 664 //////////////////////////////////////////////////////////////////////////////
648 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 665 //
649 G4double kinEnergyScatProj) << 666 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForScatProjToProjCase(G4double kinEnergyScatProj)
650 { << 667 { G4double maxEProj= HighEnergyLimit;
651 G4double maxEProj = GetHighEnergyLimit(); << 668 if (second_part_of_same_type) maxEProj=std::min(kinEnergyScatProj*2.,HighEnergyLimit);
652 if(fSecondPartSameType) <<
653 maxEProj = std::min(kinEnergyScatProj * 2. <<
654 return maxEProj; 669 return maxEProj;
655 } 670 }
656 <<
657 ////////////////////////////////////////////// 671 //////////////////////////////////////////////////////////////////////////////
658 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 672 //
659 G4double primAdjEnergy, G4double tcut) << 673 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForScatProjToProjCase(G4double PrimAdjEnergy,G4double Tcut)
660 { << 674 { G4double Emin=PrimAdjEnergy;
661 G4double Emin = primAdjEnergy; << 675 if (ApplyCutInRange) Emin=PrimAdjEnergy+Tcut;
662 if(fApplyCutInRange) <<
663 Emin += tcut; <<
664 return Emin; 676 return Emin;
665 } 677 }
666 <<
667 ////////////////////////////////////////////// 678 //////////////////////////////////////////////////////////////////////////////
668 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 679 //
669 { << 680 G4double G4VEmAdjointModel::GetSecondAdjEnergyMaxForProdToProjCase(G4double )
670 return fHighEnergyLimit; << 681 { return HighEnergyLimit;
671 } 682 }
672 <<
673 ////////////////////////////////////////////// 683 //////////////////////////////////////////////////////////////////////////////
674 G4double G4VEmAdjointModel::GetSecondAdjEnergy << 684 //
675 G4double primAdjEnergy) << 685 G4double G4VEmAdjointModel::GetSecondAdjEnergyMinForProdToProjCase(G4double PrimAdjEnergy)
676 { << 686 { G4double minEProj=PrimAdjEnergy;
677 G4double minEProj = primAdjEnergy; << 687 if (second_part_of_same_type) minEProj=PrimAdjEnergy*2.;
678 if(fSecondPartSameType) <<
679 minEProj = primAdjEnergy * 2.; <<
680 return minEProj; 688 return minEProj;
681 } 689 }
682 <<
683 ////////////////////////////////////////////// 690 ////////////////////////////////////////////////////////////////////////////////////////////
684 void G4VEmAdjointModel::DefineCurrentMaterial( << 691 //
685 const G4MaterialCutsCouple* couple) << 692 void G4VEmAdjointModel::DefineCurrentMaterial(const G4MaterialCutsCouple* couple)
686 { << 693 { if(couple != currentCouple) {
687 if(couple != fCurrentCouple) << 694 currentCouple = const_cast<G4MaterialCutsCouple*> (couple);
688 { << 695 currentMaterial = const_cast<G4Material*> (couple->GetMaterial());
689 fCurrentCouple = const_cast<G4MaterialCu << 696 currentCoupleIndex = couple->GetIndex();
690 fCurrentMaterial = const_cast<G4Material*> << 697 currentMaterialIndex = currentMaterial->GetIndex();
691 std::size_t idx = 56; << 698 size_t idx=56;
692 fTcutSecond = 1.e-11; << 699 currentTcutForDirectSecond =0.00000000001;
693 if(fAdjEquivDirectSecondPart) << 700 if (theAdjEquivOfDirectSecondPartDef) {
694 { << 701 if (theAdjEquivOfDirectSecondPartDef == G4AdjointGamma::AdjointGamma()) idx = 0;
695 if(fAdjEquivDirectSecondPart == G4Adjoin << 702 else if (theAdjEquivOfDirectSecondPartDef == G4AdjointElectron::AdjointElectron()) idx = 1;
696 idx = 0; << 703 else if (theAdjEquivOfDirectSecondPartDef == G4AdjointPositron::AdjointPositron()) idx = 2;
697 else if(fAdjEquivDirectSecondPart == G4A << 704 if (idx <56){
698 idx = 1; << 705 const std::vector<G4double>* aVec = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(idx);
699 else if(fAdjEquivDirectSecondPart == G4A << 706 currentTcutForDirectSecond=(*aVec)[currentCoupleIndex];
700 idx = 2; << 707 }
701 if(idx < 56) << 708 }
702 { << 709
703 const std::vector<G4double>* aVec = << 710
704 G4ProductionCutsTable::GetProduction <<
705 idx); <<
706 fTcutSecond = (*aVec)[couple->GetIndex <<
707 } <<
708 } <<
709 } 711 }
710 } 712 }
711 <<
712 ////////////////////////////////////////////// 713 ////////////////////////////////////////////////////////////////////////////////////////////
>> 714 //
713 void G4VEmAdjointModel::SetHighEnergyLimit(G4d 715 void G4VEmAdjointModel::SetHighEnergyLimit(G4double aVal)
714 { << 716 { HighEnergyLimit=aVal;
715 fHighEnergyLimit = aVal; << 717 if (theDirectEMModel) theDirectEMModel->SetHighEnergyLimit( aVal);
716 if(fDirectModel) <<
717 fDirectModel->SetHighEnergyLimit(aVal); <<
718 } 718 }
719 <<
720 ////////////////////////////////////////////// 719 ////////////////////////////////////////////////////////////////////////////////////////////
>> 720 //
721 void G4VEmAdjointModel::SetLowEnergyLimit(G4do 721 void G4VEmAdjointModel::SetLowEnergyLimit(G4double aVal)
722 { 722 {
723 fLowEnergyLimit = aVal; << 723 LowEnergyLimit=aVal;
724 if(fDirectModel) << 724 if (theDirectEMModel) theDirectEMModel->SetLowEnergyLimit( aVal);
725 fDirectModel->SetLowEnergyLimit(aVal); <<
726 } 725 }
727 <<
728 ////////////////////////////////////////////// 726 ////////////////////////////////////////////////////////////////////////////////////////////
729 void G4VEmAdjointModel::SetAdjointEquivalentOf << 727 //
730 G4ParticleDefinition* aPart) << 728 void G4VEmAdjointModel::SetAdjointEquivalentOfDirectPrimaryParticleDefinition(G4ParticleDefinition* aPart)
731 { 729 {
732 fAdjEquivDirectPrimPart = aPart; << 730 theAdjEquivOfDirectPrimPartDef=aPart;
733 if(fAdjEquivDirectPrimPart->GetParticleName( << 731 if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_e-")
734 fDirectPrimaryPart = G4Electron::Electron( << 732 theDirectPrimaryPartDef=G4Electron::Electron();
735 else if(fAdjEquivDirectPrimPart->GetParticle << 733 if (theAdjEquivOfDirectPrimPartDef->GetParticleName() =="adj_gamma")
736 fDirectPrimaryPart = G4Gamma::Gamma(); << 734 theDirectPrimaryPartDef=G4Gamma::Gamma();
737 } 735 }
738 736