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