| 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 //
27 27
28 #include "G4AdjointCSManager.hh" 28 #include "G4AdjointCSManager.hh"
29 29
30 #include "G4AdjointCSMatrix.hh" 30 #include "G4AdjointCSMatrix.hh"
31 #include "G4AdjointElectron.hh" 31 #include "G4AdjointElectron.hh"
32 #include "G4AdjointGamma.hh" 32 #include "G4AdjointGamma.hh"
33 #include "G4AdjointInterpolator.hh" 33 #include "G4AdjointInterpolator.hh"
34 #include "G4AdjointProton.hh" 34 #include "G4AdjointProton.hh"
35 #include "G4Electron.hh" 35 #include "G4Electron.hh"
36 #include "G4Element.hh" 36 #include "G4Element.hh"
37 #include "G4ElementTable.hh" 37 #include "G4ElementTable.hh"
38 #include "G4Gamma.hh" 38 #include "G4Gamma.hh"
39 #include "G4ParticleDefinition.hh" 39 #include "G4ParticleDefinition.hh"
40 #include "G4PhysicalConstants.hh" 40 #include "G4PhysicalConstants.hh"
41 #include "G4PhysicsLogVector.hh" 41 #include "G4PhysicsLogVector.hh"
42 #include "G4PhysicsTable.hh" 42 #include "G4PhysicsTable.hh"
43 #include "G4ProductionCutsTable.hh" 43 #include "G4ProductionCutsTable.hh"
44 #include "G4Proton.hh" 44 #include "G4Proton.hh"
45 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
46 #include "G4VEmAdjointModel.hh" 46 #include "G4VEmAdjointModel.hh"
47 #include "G4VEmProcess.hh" 47 #include "G4VEmProcess.hh"
48 #include "G4VEnergyLossProcess.hh" 48 #include "G4VEnergyLossProcess.hh"
49 49
50 G4ThreadLocal G4AdjointCSManager* G4AdjointCSM 50 G4ThreadLocal G4AdjointCSManager* G4AdjointCSManager::fInstance = nullptr;
51 51
52 constexpr G4double G4AdjointCSManager::fTmin; 52 constexpr G4double G4AdjointCSManager::fTmin;
53 constexpr G4double G4AdjointCSManager::fTmax; 53 constexpr G4double G4AdjointCSManager::fTmax;
54 constexpr G4int G4AdjointCSManager::fNbins; 54 constexpr G4int G4AdjointCSManager::fNbins;
55 55
56 ////////////////////////////////////////////// 56 ///////////////////////////////////////////////////////
57 G4AdjointCSManager* G4AdjointCSManager::GetAdj 57 G4AdjointCSManager* G4AdjointCSManager::GetAdjointCSManager()
58 { 58 {
59 if(fInstance == nullptr) 59 if(fInstance == nullptr)
60 { 60 {
61 static G4ThreadLocalSingleton<G4AdjointCSM 61 static G4ThreadLocalSingleton<G4AdjointCSManager> inst;
62 fInstance = inst.Instance(); 62 fInstance = inst.Instance();
63 } 63 }
64 return fInstance; 64 return fInstance;
65 } 65 }
66 66
67 ////////////////////////////////////////////// 67 ///////////////////////////////////////////////////////
68 G4AdjointCSManager::G4AdjointCSManager() 68 G4AdjointCSManager::G4AdjointCSManager()
69 { 69 {
70 RegisterAdjointParticle(G4AdjointElectron::A 70 RegisterAdjointParticle(G4AdjointElectron::AdjointElectron());
71 RegisterAdjointParticle(G4AdjointGamma::Adjo 71 RegisterAdjointParticle(G4AdjointGamma::AdjointGamma());
72 RegisterAdjointParticle(G4AdjointProton::Adj 72 RegisterAdjointParticle(G4AdjointProton::AdjointProton());
73 } 73 }
74 74
75 ////////////////////////////////////////////// 75 ///////////////////////////////////////////////////////
76 G4AdjointCSManager::~G4AdjointCSManager() 76 G4AdjointCSManager::~G4AdjointCSManager()
77 { 77 {
78 for (auto& p : fAdjointCSMatricesForProdToPr 78 for (auto& p : fAdjointCSMatricesForProdToProj) {
79 for (auto p1 : p) { 79 for (auto p1 : p) {
80 if (p1) { 80 if (p1) {
81 delete p1; 81 delete p1;
82 p1 = nullptr; 82 p1 = nullptr;
83 } 83 }
84 } 84 }
85 p.clear(); 85 p.clear();
86 } 86 }
87 fAdjointCSMatricesForProdToProj.clear(); 87 fAdjointCSMatricesForProdToProj.clear();
88 88
89 for (auto& p : fAdjointCSMatricesForScatProj 89 for (auto& p : fAdjointCSMatricesForScatProjToProj) {
90 for (auto p1 : p) { 90 for (auto p1 : p) {
91 if (p1) { 91 if (p1) {
92 delete p1; 92 delete p1;
93 p1 = nullptr; 93 p1 = nullptr;
94 } 94 }
95 } 95 }
96 p.clear(); 96 p.clear();
97 } 97 }
98 fAdjointCSMatricesForScatProjToProj.clear(); 98 fAdjointCSMatricesForScatProjToProj.clear();
99 99
100 for (auto p : fAdjointModels) { 100 for (auto p : fAdjointModels) {
101 if (p) { 101 if (p) {
102 delete p; 102 delete p;
103 p = nullptr; 103 p = nullptr;
104 } 104 }
105 } 105 }
106 fAdjointModels.clear(); 106 fAdjointModels.clear();
107 107
108 for (auto p : fTotalAdjSigmaTable) { 108 for (auto p : fTotalAdjSigmaTable) {
109 p->clearAndDestroy(); 109 p->clearAndDestroy();
110 delete p; 110 delete p;
111 p = nullptr; 111 p = nullptr;
112 } 112 }
113 fTotalAdjSigmaTable.clear(); 113 fTotalAdjSigmaTable.clear();
114 114
115 for (auto p : fSigmaTableForAdjointModelScat 115 for (auto p : fSigmaTableForAdjointModelScatProjToProj) {
116 p->clearAndDestroy(); 116 p->clearAndDestroy();
117 delete p; 117 delete p;
118 p = nullptr; 118 p = nullptr;
119 } 119 }
120 fSigmaTableForAdjointModelScatProjToProj.cle 120 fSigmaTableForAdjointModelScatProjToProj.clear();
121 121
122 for (auto p : fSigmaTableForAdjointModelProd 122 for (auto p : fSigmaTableForAdjointModelProdToProj) {
123 p->clearAndDestroy(); 123 p->clearAndDestroy();
124 delete p; 124 delete p;
125 p = nullptr; 125 p = nullptr;
126 } 126 }
127 fSigmaTableForAdjointModelProdToProj.clear() 127 fSigmaTableForAdjointModelProdToProj.clear();
128 128
129 for (auto p : fTotalFwdSigmaTable) { 129 for (auto p : fTotalFwdSigmaTable) {
130 p->clearAndDestroy(); 130 p->clearAndDestroy();
131 delete p; 131 delete p;
132 p = nullptr; 132 p = nullptr;
133 } 133 }
134 fTotalFwdSigmaTable.clear(); 134 fTotalFwdSigmaTable.clear();
135 135
136 for (auto p : fForwardProcesses) { 136 for (auto p : fForwardProcesses) {
137 delete p; 137 delete p;
138 p = nullptr; 138 p = nullptr;
139 } 139 }
140 fForwardProcesses.clear(); 140 fForwardProcesses.clear();
141 141
142 for (auto p : fForwardLossProcesses) { 142 for (auto p : fForwardLossProcesses) {
143 delete p; 143 delete p;
144 p = nullptr; 144 p = nullptr;
145 } 145 }
146 fForwardLossProcesses.clear(); 146 fForwardLossProcesses.clear();
147 } 147 }
148 148
149 ////////////////////////////////////////////// 149 ///////////////////////////////////////////////////////
150 std::size_t G4AdjointCSManager::RegisterEmAdjo 150 std::size_t G4AdjointCSManager::RegisterEmAdjointModel(G4VEmAdjointModel* aModel)
151 { 151 {
152 fAdjointModels.push_back(aModel); 152 fAdjointModels.push_back(aModel);
153 fSigmaTableForAdjointModelScatProjToProj.pus 153 fSigmaTableForAdjointModelScatProjToProj.push_back(new G4PhysicsTable);
154 fSigmaTableForAdjointModelProdToProj.push_ba 154 fSigmaTableForAdjointModelProdToProj.push_back(new G4PhysicsTable);
155 return fAdjointModels.size() - 1; 155 return fAdjointModels.size() - 1;
156 } 156 }
157 157
158 ////////////////////////////////////////////// 158 ///////////////////////////////////////////////////////
159 void G4AdjointCSManager::RegisterEmProcess(G4V 159 void G4AdjointCSManager::RegisterEmProcess(G4VEmProcess* aProcess,
160 G4P 160 G4ParticleDefinition* aFwdPartDef)
161 { 161 {
162 G4ParticleDefinition* anAdjPartDef = 162 G4ParticleDefinition* anAdjPartDef =
163 GetAdjointParticleEquivalent(aFwdPartDef); 163 GetAdjointParticleEquivalent(aFwdPartDef);
164 if(anAdjPartDef && aProcess) 164 if(anAdjPartDef && aProcess)
165 { 165 {
166 RegisterAdjointParticle(anAdjPartDef); 166 RegisterAdjointParticle(anAdjPartDef);
167 167
168 for(std::size_t i = 0; i < fAdjointParticl 168 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
169 { 169 {
170 if(anAdjPartDef->GetParticleName() == 170 if(anAdjPartDef->GetParticleName() ==
171 fAdjointParticlesInAction[i]->GetPart 171 fAdjointParticlesInAction[i]->GetParticleName())
172 fForwardProcesses[i]->push_back(aProce 172 fForwardProcesses[i]->push_back(aProcess);
173 } 173 }
174 } 174 }
175 } 175 }
176 176
177 ////////////////////////////////////////////// 177 ///////////////////////////////////////////////////////
178 void G4AdjointCSManager::RegisterEnergyLossPro 178 void G4AdjointCSManager::RegisterEnergyLossProcess(
179 G4VEnergyLossProcess* aProcess, G4ParticleDe 179 G4VEnergyLossProcess* aProcess, G4ParticleDefinition* aFwdPartDef)
180 { 180 {
181 G4ParticleDefinition* anAdjPartDef = 181 G4ParticleDefinition* anAdjPartDef =
182 GetAdjointParticleEquivalent(aFwdPartDef); 182 GetAdjointParticleEquivalent(aFwdPartDef);
183 if(anAdjPartDef && aProcess) 183 if(anAdjPartDef && aProcess)
184 { 184 {
185 RegisterAdjointParticle(anAdjPartDef); 185 RegisterAdjointParticle(anAdjPartDef);
186 for(std::size_t i = 0; i < fAdjointParticl 186 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
187 { 187 {
188 if(anAdjPartDef->GetParticleName() == 188 if(anAdjPartDef->GetParticleName() ==
189 fAdjointParticlesInAction[i]->GetPart 189 fAdjointParticlesInAction[i]->GetParticleName())
190 fForwardLossProc 190 fForwardLossProcesses[i]->push_back(aProcess);
191 191
192 } 192 }
193 } 193 }
194 } 194 }
195 195
196 ////////////////////////////////////////////// 196 ///////////////////////////////////////////////////////
197 void G4AdjointCSManager::RegisterAdjointPartic 197 void G4AdjointCSManager::RegisterAdjointParticle(G4ParticleDefinition* aPartDef)
198 { 198 {
199 G4bool found = false; 199 G4bool found = false;
200 for(auto p : fAdjointParticlesInAction) 200 for(auto p : fAdjointParticlesInAction)
201 { 201 {
202 if(p->GetParticleName() == aPartDef->GetPa 202 if(p->GetParticleName() == aPartDef->GetParticleName())
203 { 203 {
204 found = true; 204 found = true;
205 } 205 }
206 } 206 }
207 if(!found) 207 if(!found)
208 { 208 {
209 fForwardLossProcesses.push_back(new std::v 209 fForwardLossProcesses.push_back(new std::vector<G4VEnergyLossProcess*>());
210 fTotalFwdSigmaTable.push_back(new G4Physic 210 fTotalFwdSigmaTable.push_back(new G4PhysicsTable);
211 fTotalAdjSigmaTable.push_back(new G4Physic 211 fTotalAdjSigmaTable.push_back(new G4PhysicsTable);
212 fForwardProcesses.push_back(new std::vecto 212 fForwardProcesses.push_back(new std::vector<G4VEmProcess*>());
213 fAdjointParticlesInAction.push_back(aPartD 213 fAdjointParticlesInAction.push_back(aPartDef);
214 fEminForFwdSigmaTables.push_back(std::vect 214 fEminForFwdSigmaTables.push_back(std::vector<G4double>());
215 fEminForAdjSigmaTables.push_back(std::vect 215 fEminForAdjSigmaTables.push_back(std::vector<G4double>());
216 fEkinofFwdSigmaMax.push_back(std::vector<G 216 fEkinofFwdSigmaMax.push_back(std::vector<G4double>());
217 fEkinofAdjSigmaMax.push_back(std::vector<G 217 fEkinofAdjSigmaMax.push_back(std::vector<G4double>());
218 } 218 }
219 } 219 }
220 220
221 ////////////////////////////////////////////// 221 ///////////////////////////////////////////////////////
222 void G4AdjointCSManager::BuildCrossSectionMatr 222 void G4AdjointCSManager::BuildCrossSectionMatrices()
223 { 223 {
224 if(fCSMatricesBuilt) 224 if(fCSMatricesBuilt)
225 return; 225 return;
226 // The Tcut and Tmax matrices will be comput 226 // The Tcut and Tmax matrices will be computed probably just once.
227 // When Tcut changes, some PhysicsTable will 227 // When Tcut changes, some PhysicsTable will be recomputed
228 // for each MaterialCutCouple but not all th 228 // for each MaterialCutCouple but not all the matrices.
229 // The Tcut defines a lower limit in the ene 229 // The Tcut defines a lower limit in the energy of the projectile before
230 // scattering. In the Projectile to Scattere 230 // scattering. In the Projectile to Scattered Projectile case we have
231 // E_ScatProj<E_Proj-Tcut 231 // E_ScatProj<E_Proj-Tcut
232 // Therefore in the adjoint case we have 232 // Therefore in the adjoint case we have
233 // Eproj> E_ScatProj+Tcut 233 // Eproj> E_ScatProj+Tcut
234 // This implies that when computing the adjo 234 // This implies that when computing the adjoint CS we should integrate over
235 // Epro from E_ScatProj+Tcut to Emax 235 // Epro from E_ScatProj+Tcut to Emax
236 // In the Projectile to Secondary case Tcut 236 // In the Projectile to Secondary case Tcut plays a role only in the fact that
237 // Esecond should be greater than Tcut to ha 237 // Esecond should be greater than Tcut to have the possibility to have any
238 // adjoint process. 238 // adjoint process.
239 // To avoid recomputing matrices for all cha 239 // To avoid recomputing matrices for all changes of MaterialCutCouple,
240 // we propose to compute the matrices only o 240 // we propose to compute the matrices only once for the minimum possible Tcut
241 // and then to interpolate the probability f 241 // and then to interpolate the probability for a new Tcut (implemented in
242 // G4VAdjointEmModel) 242 // G4VAdjointEmModel)
243 243
244 fAdjointCSMatricesForScatProjToProj.clear(); 244 fAdjointCSMatricesForScatProjToProj.clear();
245 fAdjointCSMatricesForProdToProj.clear(); 245 fAdjointCSMatricesForProdToProj.clear();
246 const G4ElementTable* theElementTable = G4 246 const G4ElementTable* theElementTable = G4Element::GetElementTable();
247 const G4MaterialTable* theMaterialTable = G4 247 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
248 248
249 G4cout << "========== Computation of cross s 249 G4cout << "========== Computation of cross section matrices for adjoint "
250 "models ==========" 250 "models =========="
251 << G4endl; 251 << G4endl;
252 for(const auto& aModel : fAdjointModels) 252 for(const auto& aModel : fAdjointModels)
253 { 253 {
254 G4cout << "Build adjoint cross section mat 254 G4cout << "Build adjoint cross section matrices for " << aModel->GetName()
255 << G4endl; 255 << G4endl;
256 if(aModel->GetUseMatrix()) 256 if(aModel->GetUseMatrix())
257 { 257 {
258 std::vector<G4AdjointCSMatrix*>* aListOf 258 std::vector<G4AdjointCSMatrix*>* aListOfMat1 =
259 new std::vector<G4AdjointCSMatrix*>(); 259 new std::vector<G4AdjointCSMatrix*>();
260 std::vector<G4AdjointCSMatrix*>* aListOf 260 std::vector<G4AdjointCSMatrix*>* aListOfMat2 =
261 new std::vector<G4AdjointCSMatrix*>(); 261 new std::vector<G4AdjointCSMatrix*>();
262 if(aModel->GetUseMatrixPerElement()) 262 if(aModel->GetUseMatrixPerElement())
263 { 263 {
264 if(aModel->GetUseOnlyOneMatrixForAllEl 264 if(aModel->GetUseOnlyOneMatrixForAllElements())
265 { 265 {
266 std::vector<G4AdjointCSMatrix*> two_ 266 std::vector<G4AdjointCSMatrix*> two_matrices =
267 BuildCrossSectionsModelAndElement( 267 BuildCrossSectionsModelAndElement(aModel, 1, 1, 80);
268 aListOfMat1->push_back(two_matrices[ 268 aListOfMat1->push_back(two_matrices[0]);
269 aListOfMat2->push_back(two_matrices[ 269 aListOfMat2->push_back(two_matrices[1]);
270 } 270 }
271 else 271 else
272 { 272 {
273 for(const auto& anElement : *theElem 273 for(const auto& anElement : *theElementTable)
274 { 274 {
275 G4int Z = G4lrint(anElement->GetZ( 275 G4int Z = G4lrint(anElement->GetZ());
276 G4int A = G4lrint(anElement->GetN( 276 G4int A = G4lrint(anElement->GetN());
277 std::vector<G4AdjointCSMatrix*> tw 277 std::vector<G4AdjointCSMatrix*> two_matrices =
278 BuildCrossSectionsModelAndElemen 278 BuildCrossSectionsModelAndElement(aModel, Z, A, 40);
279 aListOfMat1->push_back(two_matrice 279 aListOfMat1->push_back(two_matrices[0]);
280 aListOfMat2->push_back(two_matrice 280 aListOfMat2->push_back(two_matrices[1]);
281 } 281 }
282 } 282 }
283 } 283 }
284 else 284 else
285 { // Per material case 285 { // Per material case
286 for(const auto& aMaterial : *theMateri 286 for(const auto& aMaterial : *theMaterialTable)
287 { 287 {
288 std::vector<G4AdjointCSMatrix*> two_ 288 std::vector<G4AdjointCSMatrix*> two_matrices =
289 BuildCrossSectionsModelAndMaterial 289 BuildCrossSectionsModelAndMaterial(aModel, aMaterial, 40);
290 aListOfMat1->push_back(two_matrices[ 290 aListOfMat1->push_back(two_matrices[0]);
291 aListOfMat2->push_back(two_matrices[ 291 aListOfMat2->push_back(two_matrices[1]);
292 } 292 }
293 } 293 }
294 fAdjointCSMatricesForProdToProj.push_bac 294 fAdjointCSMatricesForProdToProj.push_back(*aListOfMat1);
295 fAdjointCSMatricesForScatProjToProj.push 295 fAdjointCSMatricesForScatProjToProj.push_back(*aListOfMat2);
296 aModel->SetCSMatrices(aListOfMat1, aList 296 aModel->SetCSMatrices(aListOfMat1, aListOfMat2);
297 } 297 }
298 else 298 else
299 { 299 {
300 G4cout << "The model " << aModel->GetNam 300 G4cout << "The model " << aModel->GetName()
301 << " does not use cross section m 301 << " does not use cross section matrices" << G4endl;
302 std::vector<G4AdjointCSMatrix*> two_empt 302 std::vector<G4AdjointCSMatrix*> two_empty_matrices;
303 fAdjointCSMatricesForProdToProj.push_bac 303 fAdjointCSMatricesForProdToProj.push_back(two_empty_matrices);
304 fAdjointCSMatricesForScatProjToProj.push << 304 fAdjointCSMatricesForScatProjToProj.push_back(two_empty_matrices);
305 } 305 }
306 } 306 }
307 G4cout << " All adjoint cross s 307 G4cout << " All adjoint cross section matrices are computed!"
308 << G4endl; 308 << G4endl;
309 G4cout 309 G4cout
310 << "====================================== 310 << "======================================================================"
311 << G4endl; 311 << G4endl;
312 312
313 fCSMatricesBuilt = true; 313 fCSMatricesBuilt = true;
314 } 314 }
315 315
316 ////////////////////////////////////////////// 316 ///////////////////////////////////////////////////////
317 void G4AdjointCSManager::BuildTotalSigmaTables 317 void G4AdjointCSManager::BuildTotalSigmaTables()
318 { 318 {
319 if(fSigmaTableBuilt) 319 if(fSigmaTableBuilt)
320 return; 320 return;
321 321
322 const G4ProductionCutsTable* theCoupleTable 322 const G4ProductionCutsTable* theCoupleTable =
323 G4ProductionCutsTable::GetProductionCutsTa 323 G4ProductionCutsTable::GetProductionCutsTable();
324 324
325 // Prepare the Sigma table for all AdjointEM 325 // Prepare the Sigma table for all AdjointEMModel, will be filled later on
326 for(std::size_t i = 0; i < fAdjointModels.si 326 for(std::size_t i = 0; i < fAdjointModels.size(); ++i)
327 { 327 {
328 fSigmaTableForAdjointModelScatProjToProj[i 328 fSigmaTableForAdjointModelScatProjToProj[i]->clearAndDestroy();
329 fSigmaTableForAdjointModelProdToProj[i]->c 329 fSigmaTableForAdjointModelProdToProj[i]->clearAndDestroy();
330 for(std::size_t j = 0; j < theCoupleTable- 330 for(std::size_t j = 0; j < theCoupleTable->GetTableSize(); ++j)
331 { 331 {
332 fSigmaTableForAdjointModelScatProjToProj 332 fSigmaTableForAdjointModelScatProjToProj[i]->push_back(
333 new G4PhysicsLogVector(fTmin, fTmax, f 333 new G4PhysicsLogVector(fTmin, fTmax, fNbins));
334 fSigmaTableForAdjointModelProdToProj[i]- 334 fSigmaTableForAdjointModelProdToProj[i]->push_back(
335 new G4PhysicsLogVector(fTmin, fTmax, f 335 new G4PhysicsLogVector(fTmin, fTmax, fNbins));
336 } 336 }
337 } 337 }
338 338
339 for(std::size_t i = 0; i < fAdjointParticles 339 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
340 { 340 {
341 G4ParticleDefinition* thePartDef = fAdjoin 341 G4ParticleDefinition* thePartDef = fAdjointParticlesInAction[i];
342 DefineCurrentParticle(thePartDef); 342 DefineCurrentParticle(thePartDef);
343 fTotalFwdSigmaTable[i]->clearAndDestroy(); 343 fTotalFwdSigmaTable[i]->clearAndDestroy();
344 fTotalAdjSigmaTable[i]->clearAndDestroy(); 344 fTotalAdjSigmaTable[i]->clearAndDestroy();
345 fEminForFwdSigmaTables[i].clear(); 345 fEminForFwdSigmaTables[i].clear();
346 fEminForAdjSigmaTables[i].clear(); 346 fEminForAdjSigmaTables[i].clear();
347 fEkinofFwdSigmaMax[i].clear(); 347 fEkinofFwdSigmaMax[i].clear();
348 fEkinofAdjSigmaMax[i].clear(); 348 fEkinofAdjSigmaMax[i].clear();
349 349
350 for(std::size_t j = 0; j < theCoupleTable- 350 for(std::size_t j = 0; j < theCoupleTable->GetTableSize(); ++j)
351 { 351 {
352 const G4MaterialCutsCouple* couple = 352 const G4MaterialCutsCouple* couple =
353 theCoupleTable->GetMaterialCutsCouple( 353 theCoupleTable->GetMaterialCutsCouple((G4int)j);
354 354
355 // make first the total fwd CS table for 355 // make first the total fwd CS table for FwdProcess
356 G4PhysicsVector* aVector = new G4Physics 356 G4PhysicsVector* aVector = new G4PhysicsLogVector(fTmin, fTmax, fNbins);
357 G4bool Emin_found = false; 357 G4bool Emin_found = false;
358 G4double sigma_max = 0.; 358 G4double sigma_max = 0.;
359 G4double e_sigma_max = 0.; 359 G4double e_sigma_max = 0.;
360 for(std::size_t l = 0; l < fNbins; ++l) 360 for(std::size_t l = 0; l < fNbins; ++l)
361 { 361 {
362 G4double totCS = 0.; 362 G4double totCS = 0.;
363 G4double e = aVector->Energy(l); 363 G4double e = aVector->Energy(l);
364 for(std::size_t k = 0; k < fForwardPro 364 for(std::size_t k = 0; k < fForwardProcesses[i]->size(); ++k)
365 { 365 {
366 totCS += (*fForwardProcesses[i])[k]- 366 totCS += (*fForwardProcesses[i])[k]->GetCrossSection(e, couple);
367 } 367 }
368 for(std::size_t k = 0; k < fForwardLos 368 for(std::size_t k = 0; k < fForwardLossProcesses[i]->size(); ++k)
369 { 369 {
370 if(thePartDef == fAdjIon) 370 if(thePartDef == fAdjIon)
371 { // e is considered already as the 371 { // e is considered already as the scaled energy
372 std::size_t mat_index = couple->Ge 372 std::size_t mat_index = couple->GetIndex();
373 G4VEmModel* currentModel = 373 G4VEmModel* currentModel =
374 (*fForwardLossProcesses[i])[k]-> 374 (*fForwardLossProcesses[i])[k]->SelectModelForMaterial(e,
375 375 mat_index);
376 G4double chargeSqRatio = currentMo 376 G4double chargeSqRatio = currentModel->GetChargeSquareRatio(
377 fFwdIon, couple->GetMaterial(), 377 fFwdIon, couple->GetMaterial(), e / fMassRatio);
378 (*fForwardLossProcesses[i])[k]->Se 378 (*fForwardLossProcesses[i])[k]->SetDynamicMassCharge(fMassRatio,
379 379 chargeSqRatio);
380 } 380 }
381 G4double e1 = e / fMassRatio; 381 G4double e1 = e / fMassRatio;
382 totCS += (*fForwardLossProcesses[i]) 382 totCS += (*fForwardLossProcesses[i])[k]->GetLambda(e1, couple);
383 } 383 }
384 aVector->PutValue(l, totCS); 384 aVector->PutValue(l, totCS);
385 if(totCS > sigma_max) 385 if(totCS > sigma_max)
386 { 386 {
387 sigma_max = totCS; 387 sigma_max = totCS;
388 e_sigma_max = e; 388 e_sigma_max = e;
389 } 389 }
390 if(totCS > 0 && !Emin_found) 390 if(totCS > 0 && !Emin_found)
391 { 391 {
392 fEminForFwdSigmaTables[i].push_back( 392 fEminForFwdSigmaTables[i].push_back(e);
393 Emin_found = true; 393 Emin_found = true;
394 } 394 }
395 } 395 }
396 396
397 fEkinofFwdSigmaMax[i].push_back(e_sigma_ 397 fEkinofFwdSigmaMax[i].push_back(e_sigma_max);
398 398
399 if(!Emin_found) 399 if(!Emin_found)
400 fEminForFwdSigmaTables[i].push_back(fT 400 fEminForFwdSigmaTables[i].push_back(fTmax);
401 401
402 fTotalFwdSigmaTable[i]->push_back(aVecto 402 fTotalFwdSigmaTable[i]->push_back(aVector);
403 403
404 Emin_found = false; 404 Emin_found = false;
405 sigma_max = 0; 405 sigma_max = 0;
406 e_sigma_max = 0.; 406 e_sigma_max = 0.;
407 G4PhysicsVector* aVector1 = new G4Physic 407 G4PhysicsVector* aVector1 = new G4PhysicsLogVector(fTmin, fTmax, fNbins);
408 for(std::size_t eindex = 0; eindex < fNb 408 for(std::size_t eindex = 0; eindex < fNbins; ++eindex)
409 { 409 {
410 G4double e = aVector1->Energy(eind 410 G4double e = aVector1->Energy(eindex);
411 G4double totCS = ComputeTotalAdjointCS 411 G4double totCS = ComputeTotalAdjointCS(
412 couple, thePartDef, 412 couple, thePartDef,
413 e * 0.9999999 / fMassRatio); // fMa 413 e * 0.9999999 / fMassRatio); // fMassRatio needed for ions
414 aVector1->PutValue(eindex, totCS); 414 aVector1->PutValue(eindex, totCS);
415 if(totCS > sigma_max) 415 if(totCS > sigma_max)
416 { 416 {
417 sigma_max = totCS; 417 sigma_max = totCS;
418 e_sigma_max = e; 418 e_sigma_max = e;
419 } 419 }
420 if(totCS > 0 && !Emin_found) 420 if(totCS > 0 && !Emin_found)
421 { 421 {
422 fEminForAdjSigmaTables[i].push_back( 422 fEminForAdjSigmaTables[i].push_back(e);
423 Emin_found = true; 423 Emin_found = true;
424 } 424 }
425 } 425 }
426 fEkinofAdjSigmaMax[i].push_back(e_sigma_ 426 fEkinofAdjSigmaMax[i].push_back(e_sigma_max);
427 if(!Emin_found) 427 if(!Emin_found)
428 fEminForAdjSigmaTables[i].push_back(fT 428 fEminForAdjSigmaTables[i].push_back(fTmax);
429 429
430 fTotalAdjSigmaTable[i]->push_back(aVecto 430 fTotalAdjSigmaTable[i]->push_back(aVector1);
431 } 431 }
432 } 432 }
433 fSigmaTableBuilt = true; 433 fSigmaTableBuilt = true;
434 } 434 }
435 435
436 ////////////////////////////////////////////// 436 ///////////////////////////////////////////////////////
437 G4double G4AdjointCSManager::GetTotalAdjointCS 437 G4double G4AdjointCSManager::GetTotalAdjointCS(
438 G4ParticleDefinition* aPartDef, G4double Eki 438 G4ParticleDefinition* aPartDef, G4double Ekin,
439 const G4MaterialCutsCouple* aCouple) 439 const G4MaterialCutsCouple* aCouple)
440 { 440 {
441 DefineCurrentMaterial(aCouple); 441 DefineCurrentMaterial(aCouple);
442 DefineCurrentParticle(aPartDef); 442 DefineCurrentParticle(aPartDef);
443 return (((*fTotalAdjSigmaTable[fCurrentParti 443 return (((*fTotalAdjSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
444 ->Value(Ekin * fMassRatio)); 444 ->Value(Ekin * fMassRatio));
445 } 445 }
446 446
447 ////////////////////////////////////////////// 447 ///////////////////////////////////////////////////////
448 G4double G4AdjointCSManager::GetTotalForwardCS 448 G4double G4AdjointCSManager::GetTotalForwardCS(
449 G4ParticleDefinition* aPartDef, G4double Eki 449 G4ParticleDefinition* aPartDef, G4double Ekin,
450 const G4MaterialCutsCouple* aCouple) 450 const G4MaterialCutsCouple* aCouple)
451 { 451 {
452 DefineCurrentMaterial(aCouple); 452 DefineCurrentMaterial(aCouple);
453 DefineCurrentParticle(aPartDef); 453 DefineCurrentParticle(aPartDef);
454 return (((*fTotalFwdSigmaTable[fCurrentParti 454 return (((*fTotalFwdSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
455 ->Value(Ekin * fMassRatio)); 455 ->Value(Ekin * fMassRatio));
456 } 456 }
457 457
458 ////////////////////////////////////////////// 458 ///////////////////////////////////////////////////////
459 G4double G4AdjointCSManager::GetAdjointSigma( 459 G4double G4AdjointCSManager::GetAdjointSigma(
460 G4double Ekin_nuc, std::size_t index_model, 460 G4double Ekin_nuc, std::size_t index_model, G4bool is_scat_proj_to_proj,
461 const G4MaterialCutsCouple* aCouple) 461 const G4MaterialCutsCouple* aCouple)
462 { 462 {
463 DefineCurrentMaterial(aCouple); 463 DefineCurrentMaterial(aCouple);
464 if(is_scat_proj_to_proj) 464 if(is_scat_proj_to_proj)
465 return (((*fSigmaTableForAdjointModelScatP 465 return (((*fSigmaTableForAdjointModelScatProjToProj[index_model])
466 [fCurrentMatIndex])->Value(Ekin 466 [fCurrentMatIndex])->Value(Ekin_nuc));
467 else 467 else
468 return ( 468 return (
469 ((*fSigmaTableForAdjointModelProdToProj[ 469 ((*fSigmaTableForAdjointModelProdToProj[index_model])[fCurrentMatIndex])
470 ->Value(Ekin_nuc)); 470 ->Value(Ekin_nuc));
471 } 471 }
472 472
473 ////////////////////////////////////////////// 473 ///////////////////////////////////////////////////////
474 void G4AdjointCSManager::GetEminForTotalCS(G4P 474 void G4AdjointCSManager::GetEminForTotalCS(G4ParticleDefinition* aPartDef,
475 con 475 const G4MaterialCutsCouple* aCouple,
476 G4d 476 G4double& emin_adj,
477 G4d 477 G4double& emin_fwd)
478 { 478 {
479 DefineCurrentMaterial(aCouple); 479 DefineCurrentMaterial(aCouple);
480 DefineCurrentParticle(aPartDef); 480 DefineCurrentParticle(aPartDef);
481 emin_adj = fEminForAdjSigmaTables[fCurrentPa 481 emin_adj = fEminForAdjSigmaTables[fCurrentParticleIndex][fCurrentMatIndex] /
482 fMassRatio; 482 fMassRatio;
483 emin_fwd = fEminForFwdSigmaTables[fCurrentPa 483 emin_fwd = fEminForFwdSigmaTables[fCurrentParticleIndex][fCurrentMatIndex] /
484 fMassRatio; 484 fMassRatio;
485 } 485 }
486 486
487 ////////////////////////////////////////////// 487 ///////////////////////////////////////////////////////
488 void G4AdjointCSManager::GetMaxFwdTotalCS(G4Pa 488 void G4AdjointCSManager::GetMaxFwdTotalCS(G4ParticleDefinition* aPartDef,
489 cons 489 const G4MaterialCutsCouple* aCouple,
490 G4do 490 G4double& e_sigma_max,
491 G4do 491 G4double& sigma_max)
492 { 492 {
493 DefineCurrentMaterial(aCouple); 493 DefineCurrentMaterial(aCouple);
494 DefineCurrentParticle(aPartDef); 494 DefineCurrentParticle(aPartDef);
495 e_sigma_max = fEkinofFwdSigmaMax[fCurrentPar 495 e_sigma_max = fEkinofFwdSigmaMax[fCurrentParticleIndex][fCurrentMatIndex];
496 sigma_max = ((*fTotalFwdSigmaTable[fCurrentP 496 sigma_max = ((*fTotalFwdSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
497 ->Value(e_sigma_max); 497 ->Value(e_sigma_max);
498 e_sigma_max /= fMassRatio; 498 e_sigma_max /= fMassRatio;
499 } 499 }
500 500
501 ////////////////////////////////////////////// 501 ///////////////////////////////////////////////////////
502 void G4AdjointCSManager::GetMaxAdjTotalCS(G4Pa 502 void G4AdjointCSManager::GetMaxAdjTotalCS(G4ParticleDefinition* aPartDef,
503 cons 503 const G4MaterialCutsCouple* aCouple,
504 G4do 504 G4double& e_sigma_max,
505 G4do 505 G4double& sigma_max)
506 { 506 {
507 DefineCurrentMaterial(aCouple); 507 DefineCurrentMaterial(aCouple);
508 DefineCurrentParticle(aPartDef); 508 DefineCurrentParticle(aPartDef);
509 e_sigma_max = fEkinofAdjSigmaMax[fCurrentPar 509 e_sigma_max = fEkinofAdjSigmaMax[fCurrentParticleIndex][fCurrentMatIndex];
510 sigma_max = ((*fTotalAdjSigmaTable[fCurrentP 510 sigma_max = ((*fTotalAdjSigmaTable[fCurrentParticleIndex])[fCurrentMatIndex])
511 ->Value(e_sigma_max); 511 ->Value(e_sigma_max);
512 e_sigma_max /= fMassRatio; 512 e_sigma_max /= fMassRatio;
513 } 513 }
514 514
515 ////////////////////////////////////////////// 515 ///////////////////////////////////////////////////////
516 G4double G4AdjointCSManager::GetCrossSectionCo 516 G4double G4AdjointCSManager::GetCrossSectionCorrection(
517 G4ParticleDefinition* aPartDef, G4double Pre 517 G4ParticleDefinition* aPartDef, G4double PreStepEkin,
518 const G4MaterialCutsCouple* aCouple, G4bool& 518 const G4MaterialCutsCouple* aCouple, G4bool& fwd_is_used)
519 { 519 {
520 static G4double lastEkin = 0.; 520 static G4double lastEkin = 0.;
521 static G4ParticleDefinition* lastPartDef; 521 static G4ParticleDefinition* lastPartDef;
522 522
523 G4double corr_fac = 1.; 523 G4double corr_fac = 1.;
524 if(fForwardCSMode && aPartDef) 524 if(fForwardCSMode && aPartDef)
525 { 525 {
526 if(lastEkin != PreStepEkin || aPartDef != 526 if(lastEkin != PreStepEkin || aPartDef != lastPartDef ||
527 aCouple != fCurrentCouple) 527 aCouple != fCurrentCouple)
528 { 528 {
529 DefineCurrentMaterial(aCouple); 529 DefineCurrentMaterial(aCouple);
530 G4double preadjCS = GetTotalAdjointCS(aP 530 G4double preadjCS = GetTotalAdjointCS(aPartDef, PreStepEkin, aCouple);
531 G4double prefwdCS = GetTotalForwardCS(aP 531 G4double prefwdCS = GetTotalForwardCS(aPartDef, PreStepEkin, aCouple);
532 lastEkin = PreStepEkin; 532 lastEkin = PreStepEkin;
533 lastPartDef = aPartDef; 533 lastPartDef = aPartDef;
534 if(prefwdCS > 0. && preadjCS > 0.) 534 if(prefwdCS > 0. && preadjCS > 0.)
535 { 535 {
536 fForwardCSUsed = true; 536 fForwardCSUsed = true;
537 fLastCSCorrectionFactor = prefwdCS / p 537 fLastCSCorrectionFactor = prefwdCS / preadjCS;
538 } 538 }
539 else 539 else
540 { 540 {
541 fForwardCSUsed = false; 541 fForwardCSUsed = false;
542 fLastCSCorrectionFactor = 1.; 542 fLastCSCorrectionFactor = 1.;
543 } 543 }
544 } 544 }
545 corr_fac = fLastCSCorrectionFactor; 545 corr_fac = fLastCSCorrectionFactor;
546 } 546 }
547 else 547 else
548 { 548 {
549 fForwardCSUsed = false; 549 fForwardCSUsed = false;
550 fLastCSCorrectionFactor = 1.; 550 fLastCSCorrectionFactor = 1.;
551 } 551 }
552 fwd_is_used = fForwardCSUsed; 552 fwd_is_used = fForwardCSUsed;
553 return corr_fac; 553 return corr_fac;
554 } 554 }
555 555
556 ////////////////////////////////////////////// 556 ///////////////////////////////////////////////////////
557 G4double G4AdjointCSManager::GetContinuousWeig 557 G4double G4AdjointCSManager::GetContinuousWeightCorrection(
558 G4ParticleDefinition* aPartDef, G4double Pre 558 G4ParticleDefinition* aPartDef, G4double PreStepEkin, G4double AfterStepEkin,
559 const G4MaterialCutsCouple* aCouple, G4doubl 559 const G4MaterialCutsCouple* aCouple, G4double step_length)
560 { 560 {
561 G4double corr_fac = 1.; 561 G4double corr_fac = 1.;
562 G4double after_fwdCS = GetTotalForwardCS(aPa 562 G4double after_fwdCS = GetTotalForwardCS(aPartDef, AfterStepEkin, aCouple);
563 G4double pre_adjCS = GetTotalAdjointCS(aPa 563 G4double pre_adjCS = GetTotalAdjointCS(aPartDef, PreStepEkin, aCouple);
564 if(!fForwardCSUsed || pre_adjCS == 0. || aft 564 if(!fForwardCSUsed || pre_adjCS == 0. || after_fwdCS == 0.)
565 { 565 {
566 G4double pre_fwdCS = GetTotalForwardCS(aPa 566 G4double pre_fwdCS = GetTotalForwardCS(aPartDef, PreStepEkin, aCouple);
567 corr_fac *= std::exp((pre_adjCS - pre_fwdC 567 corr_fac *= std::exp((pre_adjCS - pre_fwdCS) * step_length);
568 fLastCSCorrectionFactor = 1.; 568 fLastCSCorrectionFactor = 1.;
569 } 569 }
570 else 570 else
571 { 571 {
572 fLastCSCorrectionFactor = after_fwdCS / pr 572 fLastCSCorrectionFactor = after_fwdCS / pre_adjCS;
573 } 573 }
574 return corr_fac; 574 return corr_fac;
575 } 575 }
576 576
577 ////////////////////////////////////////////// 577 ///////////////////////////////////////////////////////
578 G4double G4AdjointCSManager::GetPostStepWeight 578 G4double G4AdjointCSManager::GetPostStepWeightCorrection()
579 { 579 {
580 return 1. / fLastCSCorrectionFactor; 580 return 1. / fLastCSCorrectionFactor;
581 } 581 }
582 582
583 ////////////////////////////////////////////// 583 ///////////////////////////////////////////////////////
584 G4double G4AdjointCSManager::ComputeAdjointCS( 584 G4double G4AdjointCSManager::ComputeAdjointCS(
585 G4Material* aMaterial, G4VEmAdjointModel* aM 585 G4Material* aMaterial, G4VEmAdjointModel* aModel, G4double PrimEnergy,
586 G4double Tcut, G4bool isScatProjToProj, std: 586 G4double Tcut, G4bool isScatProjToProj, std::vector<G4double>& CS_Vs_Element)
587 { 587 {
588 G4double EminSec = 0.; 588 G4double EminSec = 0.;
589 G4double EmaxSec = 0.; 589 G4double EmaxSec = 0.;
590 590
591 static G4double lastPrimaryEnergy = 0.; 591 static G4double lastPrimaryEnergy = 0.;
592 static G4double lastTcut = 0.; 592 static G4double lastTcut = 0.;
593 static G4Material* lastMaterial = nullptr; 593 static G4Material* lastMaterial = nullptr;
594 594
595 if(isScatProjToProj) 595 if(isScatProjToProj)
596 { 596 {
597 EminSec = aModel->GetSecondAdjEnergyMinFor 597 EminSec = aModel->GetSecondAdjEnergyMinForScatProjToProj(PrimEnergy, Tcut);
598 EmaxSec = aModel->GetSecondAdjEnergyMaxFor 598 EmaxSec = aModel->GetSecondAdjEnergyMaxForScatProjToProj(PrimEnergy);
599 } 599 }
600 else if(PrimEnergy > Tcut || !aModel->GetApp 600 else if(PrimEnergy > Tcut || !aModel->GetApplyCutInRange())
601 { 601 {
602 EminSec = aModel->GetSecondAdjEnergyMinFor 602 EminSec = aModel->GetSecondAdjEnergyMinForProdToProj(PrimEnergy);
603 EmaxSec = aModel->GetSecondAdjEnergyMaxFor 603 EmaxSec = aModel->GetSecondAdjEnergyMaxForProdToProj(PrimEnergy);
604 } 604 }
605 if(EminSec >= EmaxSec) 605 if(EminSec >= EmaxSec)
606 return 0.; 606 return 0.;
607 607
608 G4bool need_to_compute = false; 608 G4bool need_to_compute = false;
609 if(aMaterial != lastMaterial || PrimEnergy ! 609 if(aMaterial != lastMaterial || PrimEnergy != lastPrimaryEnergy ||
610 Tcut != lastTcut) 610 Tcut != lastTcut)
611 { 611 {
612 lastMaterial = aMaterial; 612 lastMaterial = aMaterial;
613 lastPrimaryEnergy = PrimEnergy; 613 lastPrimaryEnergy = PrimEnergy;
614 lastTcut = Tcut; 614 lastTcut = Tcut;
615 fIndexOfAdjointEMModelInAction.clear(); 615 fIndexOfAdjointEMModelInAction.clear();
616 fIsScatProjToProj.clear(); 616 fIsScatProjToProj.clear();
617 fLastAdjointCSVsModelsAndElements.clear(); 617 fLastAdjointCSVsModelsAndElements.clear();
618 need_to_compute = true; 618 need_to_compute = true;
619 } 619 }
620 620
621 std::size_t ind = 0; 621 std::size_t ind = 0;
622 if(!need_to_compute) 622 if(!need_to_compute)
623 { 623 {
624 need_to_compute = true; 624 need_to_compute = true;
625 for(std::size_t i = 0; i < fIndexOfAdjoint 625 for(std::size_t i = 0; i < fIndexOfAdjointEMModelInAction.size(); ++i)
626 { 626 {
627 std::size_t ind1 = fIndexOfAdjointEMMode 627 std::size_t ind1 = fIndexOfAdjointEMModelInAction[i];
628 if(aModel == fAdjointModels[ind1] && 628 if(aModel == fAdjointModels[ind1] &&
629 isScatProjToProj == fIsScatProjToProj 629 isScatProjToProj == fIsScatProjToProj[i])
630 { 630 {
631 need_to_compute = false; 631 need_to_compute = false;
632 CS_Vs_Element = fLastAdjointCSVsMode 632 CS_Vs_Element = fLastAdjointCSVsModelsAndElements[ind];
633 } 633 }
634 ++ind; 634 ++ind;
635 } 635 }
636 } 636 }
637 637
638 if(need_to_compute) 638 if(need_to_compute)
639 { 639 {
640 std::size_t ind_model = 0; 640 std::size_t ind_model = 0;
641 for(std::size_t i = 0; i < fAdjointModels. 641 for(std::size_t i = 0; i < fAdjointModels.size(); ++i)
642 { 642 {
643 if(aModel == fAdjointModels[i]) 643 if(aModel == fAdjointModels[i])
644 { 644 {
645 ind_model = i; 645 ind_model = i;
646 break; 646 break;
647 } 647 }
648 } 648 }
649 G4double Tlow = Tcut; 649 G4double Tlow = Tcut;
650 if(!fAdjointModels[ind_model]->GetApplyCut 650 if(!fAdjointModels[ind_model]->GetApplyCutInRange())
651 Tlow = fAdjointModels[ind_model]->GetLow 651 Tlow = fAdjointModels[ind_model]->GetLowEnergyLimit();
652 fIndexOfAdjointEMModelInAction.push_back(i 652 fIndexOfAdjointEMModelInAction.push_back(ind_model);
653 fIsScatProjToProj.push_back(isScatProjToPr 653 fIsScatProjToProj.push_back(isScatProjToProj);
654 CS_Vs_Element.clear(); 654 CS_Vs_Element.clear();
655 if(!aModel->GetUseMatrix()) 655 if(!aModel->GetUseMatrix())
656 { 656 {
657 CS_Vs_Element.push_back(aModel->AdjointC 657 CS_Vs_Element.push_back(aModel->AdjointCrossSection(
658 fCurrentCouple, PrimEnergy, isScatProj 658 fCurrentCouple, PrimEnergy, isScatProjToProj));
659 } 659 }
660 else if(aModel->GetUseMatrixPerElement()) 660 else if(aModel->GetUseMatrixPerElement())
661 { 661 {
662 std::size_t n_el = aMaterial->GetNumberO 662 std::size_t n_el = aMaterial->GetNumberOfElements();
663 if(aModel->GetUseOnlyOneMatrixForAllElem 663 if(aModel->GetUseOnlyOneMatrixForAllElements())
664 { 664 {
665 G4AdjointCSMatrix* theCSMatrix; 665 G4AdjointCSMatrix* theCSMatrix;
666 if(isScatProjToProj) 666 if(isScatProjToProj)
667 { 667 {
668 theCSMatrix = fAdjointCSMatricesForS 668 theCSMatrix = fAdjointCSMatricesForScatProjToProj[ind_model][0];
669 } 669 }
670 else 670 else
671 theCSMatrix = fAdjointCSMatricesForP 671 theCSMatrix = fAdjointCSMatricesForProdToProj[ind_model][0];
672 G4double CS = 0.; 672 G4double CS = 0.;
673 if(PrimEnergy > Tlow) 673 if(PrimEnergy > Tlow)
674 CS = ComputeAdjointCS(PrimEnergy, th 674 CS = ComputeAdjointCS(PrimEnergy, theCSMatrix, Tlow);
675 G4double factor = 0.; 675 G4double factor = 0.;
676 for(G4int i = 0; i < (G4int)n_el; ++i) 676 for(G4int i = 0; i < (G4int)n_el; ++i)
677 { // this could be computed only once 677 { // this could be computed only once
678 factor += aMaterial->GetElement(i)-> 678 factor += aMaterial->GetElement(i)->GetZ() *
679 aMaterial->GetVecNbOfAtoms 679 aMaterial->GetVecNbOfAtomsPerVolume()[i];
680 } 680 }
681 CS *= factor; 681 CS *= factor;
682 CS_Vs_Element.push_back(CS); 682 CS_Vs_Element.push_back(CS);
683 } 683 }
684 else 684 else
685 { 685 {
686 for(G4int i = 0; i < (G4int)n_el; ++i) 686 for(G4int i = 0; i < (G4int)n_el; ++i)
687 { 687 {
688 std::size_t ind_el = aMaterial->GetE 688 std::size_t ind_el = aMaterial->GetElement(i)->GetIndex();
689 G4AdjointCSMatrix* theCSMatrix; 689 G4AdjointCSMatrix* theCSMatrix;
690 if(isScatProjToProj) 690 if(isScatProjToProj)
691 { 691 {
692 theCSMatrix = 692 theCSMatrix =
693 fAdjointCSMatricesForScatProjToP 693 fAdjointCSMatricesForScatProjToProj[ind_model][ind_el];
694 } 694 }
695 else 695 else
696 theCSMatrix = fAdjointCSMatricesFo 696 theCSMatrix = fAdjointCSMatricesForProdToProj[ind_model][ind_el];
697 G4double CS = 0.; 697 G4double CS = 0.;
698 if(PrimEnergy > Tlow) 698 if(PrimEnergy > Tlow)
699 CS = ComputeAdjointCS(PrimEnergy, 699 CS = ComputeAdjointCS(PrimEnergy, theCSMatrix, Tlow);
700 CS_Vs_Element.push_back(CS * 700 CS_Vs_Element.push_back(CS *
701 (aMaterial-> 701 (aMaterial->GetVecNbOfAtomsPerVolume()[i]));
702 } 702 }
703 } 703 }
704 } 704 }
705 else 705 else
706 { 706 {
707 std::size_t ind_mat = aMaterial->GetInde 707 std::size_t ind_mat = aMaterial->GetIndex();
708 G4AdjointCSMatrix* theCSMatrix; 708 G4AdjointCSMatrix* theCSMatrix;
709 if(isScatProjToProj) 709 if(isScatProjToProj)
710 { 710 {
711 theCSMatrix = fAdjointCSMatricesForSca 711 theCSMatrix = fAdjointCSMatricesForScatProjToProj[ind_model][ind_mat];
712 } 712 }
713 else 713 else
714 theCSMatrix = fAdjointCSMatricesForPro 714 theCSMatrix = fAdjointCSMatricesForProdToProj[ind_model][ind_mat];
715 G4double CS = 0.; 715 G4double CS = 0.;
716 if(PrimEnergy > Tlow) 716 if(PrimEnergy > Tlow)
717 CS = ComputeAdjointCS(PrimEnergy, theC 717 CS = ComputeAdjointCS(PrimEnergy, theCSMatrix, Tlow);
718 CS_Vs_Element.push_back(CS); 718 CS_Vs_Element.push_back(CS);
719 } 719 }
720 fLastAdjointCSVsModelsAndElements.push_bac 720 fLastAdjointCSVsModelsAndElements.push_back(CS_Vs_Element);
721 } 721 }
722 722
723 G4double CS = 0.; 723 G4double CS = 0.;
724 for(const auto& cs_vs_el : CS_Vs_Element) 724 for(const auto& cs_vs_el : CS_Vs_Element)
725 { 725 {
726 // We could put the progressive sum of the 726 // We could put the progressive sum of the CS instead of the CS of an
727 // element itself 727 // element itself
728 CS += cs_vs_el; 728 CS += cs_vs_el;
729 } 729 }
730 return CS; 730 return CS;
731 } 731 }
732 732
733 ////////////////////////////////////////////// 733 ///////////////////////////////////////////////////////
734 G4Element* G4AdjointCSManager::SampleElementFr 734 G4Element* G4AdjointCSManager::SampleElementFromCSMatrices(
735 G4Material* aMaterial, G4VEmAdjointModel* aM 735 G4Material* aMaterial, G4VEmAdjointModel* aModel, G4double PrimEnergy,
736 G4double Tcut, G4bool isScatProjToProj) 736 G4double Tcut, G4bool isScatProjToProj)
737 { 737 {
738 std::vector<G4double> CS_Vs_Element; 738 std::vector<G4double> CS_Vs_Element;
739 G4double CS = ComputeAdjointCS(aMaterial, 739 G4double CS = ComputeAdjointCS(aMaterial, aModel, PrimEnergy, Tcut,
740 isScatProjToP 740 isScatProjToProj, CS_Vs_Element);
741 G4double SumCS = 0.; 741 G4double SumCS = 0.;
742 std::size_t ind = 0; 742 std::size_t ind = 0;
743 for(std::size_t i = 0; i < CS_Vs_Element.siz 743 for(std::size_t i = 0; i < CS_Vs_Element.size(); ++i)
744 { 744 {
745 SumCS += CS_Vs_Element[i]; 745 SumCS += CS_Vs_Element[i];
746 if(G4UniformRand() <= SumCS / CS) 746 if(G4UniformRand() <= SumCS / CS)
747 { 747 {
748 ind = i; 748 ind = i;
749 break; 749 break;
750 } 750 }
751 } 751 }
752 752
753 return const_cast<G4Element*>(aMaterial->Get 753 return const_cast<G4Element*>(aMaterial->GetElement((G4int)ind));
754 } 754 }
755 755
756 ////////////////////////////////////////////// 756 ///////////////////////////////////////////////////////
757 G4double G4AdjointCSManager::ComputeTotalAdjoi 757 G4double G4AdjointCSManager::ComputeTotalAdjointCS(
758 const G4MaterialCutsCouple* aCouple, G4Parti 758 const G4MaterialCutsCouple* aCouple, G4ParticleDefinition* aPartDef,
759 G4double Ekin) 759 G4double Ekin)
760 { 760 {
761 G4double TotalCS = 0.; 761 G4double TotalCS = 0.;
762 762
763 DefineCurrentMaterial(aCouple); 763 DefineCurrentMaterial(aCouple);
764 764
765 std::vector<G4double> CS_Vs_Element; 765 std::vector<G4double> CS_Vs_Element;
766 G4double CS; 766 G4double CS;
767 G4VEmAdjointModel* adjModel = nullptr; 767 G4VEmAdjointModel* adjModel = nullptr;
768 for(std::size_t i = 0; i < fAdjointModels.si 768 for(std::size_t i = 0; i < fAdjointModels.size(); ++i)
769 { 769 {
770 G4double Tlow = 0.; 770 G4double Tlow = 0.;
771 adjModel = fAdjointModels[i]; 771 adjModel = fAdjointModels[i];
772 if(!adjModel->GetApplyCutInRange()) 772 if(!adjModel->GetApplyCutInRange())
773 Tlow = adjModel->GetLowEnergyLimit(); 773 Tlow = adjModel->GetLowEnergyLimit();
774 else 774 else
775 { 775 {
776 G4ParticleDefinition* theDirSecondPartDe 776 G4ParticleDefinition* theDirSecondPartDef = GetForwardParticleEquivalent(
777 adjModel->GetAdjointEquivalentOfDirect 777 adjModel->GetAdjointEquivalentOfDirectSecondaryParticleDefinition());
778 std::size_t idx = 56; 778 std::size_t idx = 56;
779 if(theDirSecondPartDef->GetParticleName( 779 if(theDirSecondPartDef->GetParticleName() == "gamma")
780 idx = 0; 780 idx = 0;
781 else if(theDirSecondPartDef->GetParticle 781 else if(theDirSecondPartDef->GetParticleName() == "e-")
782 idx = 1; 782 idx = 1;
783 else if(theDirSecondPartDef->GetParticle 783 else if(theDirSecondPartDef->GetParticleName() == "e+")
784 idx = 2; 784 idx = 2;
785 if(idx < 56) 785 if(idx < 56)
786 { 786 {
787 const std::vector<G4double>* aVec = 787 const std::vector<G4double>* aVec =
788 G4ProductionCutsTable::GetProduction 788 G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(
789 idx); 789 idx);
790 Tlow = (*aVec)[aCouple->GetIndex()]; 790 Tlow = (*aVec)[aCouple->GetIndex()];
791 } 791 }
792 } 792 }
793 if(Ekin <= adjModel->GetHighEnergyLimit() 793 if(Ekin <= adjModel->GetHighEnergyLimit() &&
794 Ekin >= adjModel->GetLowEnergyLimit()) 794 Ekin >= adjModel->GetLowEnergyLimit())
795 { 795 {
796 if(aPartDef == 796 if(aPartDef ==
797 adjModel->GetAdjointEquivalentOfDirec 797 adjModel->GetAdjointEquivalentOfDirectPrimaryParticleDefinition())
798 { 798 {
799 CS = ComputeAdjointCS(fCurrentMaterial 799 CS = ComputeAdjointCS(fCurrentMaterial, adjModel, Ekin, Tlow, true,
800 CS_Vs_Element); 800 CS_Vs_Element);
801 TotalCS += CS; 801 TotalCS += CS;
802 (*fSigmaTableForAdjointModelScatProjTo 802 (*fSigmaTableForAdjointModelScatProjToProj[i])[fCurrentMatIndex]
803 ->PutValue(fNbins, CS); 803 ->PutValue(fNbins, CS);
804 } 804 }
805 if(aPartDef == 805 if(aPartDef ==
806 adjModel->GetAdjointEquivalentOfDirec 806 adjModel->GetAdjointEquivalentOfDirectSecondaryParticleDefinition())
807 { 807 {
808 CS = ComputeAdjointCS(fCurrentMaterial 808 CS = ComputeAdjointCS(fCurrentMaterial, adjModel, Ekin, Tlow, false,
809 CS_Vs_Element); 809 CS_Vs_Element);
810 TotalCS += CS; 810 TotalCS += CS;
811 (*fSigmaTableForAdjointModelProdToProj 811 (*fSigmaTableForAdjointModelProdToProj[i])[fCurrentMatIndex]->PutValue(
812 fNbins, CS); 812 fNbins, CS);
813 } 813 }
814 } 814 }
815 else 815 else
816 { 816 {
817 (*fSigmaTableForAdjointModelScatProjToPr 817 (*fSigmaTableForAdjointModelScatProjToProj[i])[fCurrentMatIndex]
818 ->PutValue(fNbins, 0.); 818 ->PutValue(fNbins, 0.);
819 (*fSigmaTableForAdjointModelProdToProj[i 819 (*fSigmaTableForAdjointModelProdToProj[i])[fCurrentMatIndex]->PutValue(
820 fNbins, 0.); 820 fNbins, 0.);
821 } 821 }
822 } 822 }
823 return TotalCS; 823 return TotalCS;
824 } 824 }
825 825
826 ////////////////////////////////////////////// 826 ///////////////////////////////////////////////////////
827 std::vector<G4AdjointCSMatrix*> 827 std::vector<G4AdjointCSMatrix*>
828 G4AdjointCSManager::BuildCrossSectionsModelAnd 828 G4AdjointCSManager::BuildCrossSectionsModelAndElement(G4VEmAdjointModel* aModel,
829 829 G4int Z, G4int A,
830 830 G4int nbin_pro_decade)
831 { 831 {
832 G4AdjointCSMatrix* theCSMatForProdToProjBack 832 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering =
833 new G4AdjointCSMatrix(false); 833 new G4AdjointCSMatrix(false);
834 G4AdjointCSMatrix* theCSMatForScatProjToProj 834 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering =
835 new G4AdjointCSMatrix(true); 835 new G4AdjointCSMatrix(true);
836 836
837 // make the vector of primary energy of the 837 // make the vector of primary energy of the adjoint particle.
838 G4double EkinMin = aModel->GetLowEner 838 G4double EkinMin = aModel->GetLowEnergyLimit();
839 G4double EkinMaxForScat = aModel->GetHighEne 839 G4double EkinMaxForScat = aModel->GetHighEnergyLimit() * 0.999;
840 G4double EkinMaxForProd = aModel->GetHighEne 840 G4double EkinMaxForProd = aModel->GetHighEnergyLimit() * 0.999;
841 if(aModel->GetSecondPartOfSameType()) 841 if(aModel->GetSecondPartOfSameType())
842 EkinMaxForProd = EkinMaxForProd / 2.; 842 EkinMaxForProd = EkinMaxForProd / 2.;
843 843
844 // Product to projectile backward scattering 844 // Product to projectile backward scattering
845 G4double dE = std::pow(10., 1. / nbin_pro_de 845 G4double dE = std::pow(10., 1. / nbin_pro_decade);
846 G4double E2 = 846 G4double E2 =
847 std::pow(10., double(int(std::log10(EkinMi 847 std::pow(10., double(int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
848 nbin_pro_decade) / dE; 848 nbin_pro_decade) / dE;
849 G4double E1 = EkinMin; 849 G4double E1 = EkinMin;
850 // Loop checking, 07-Aug-2015, Vladimir Ivan 850 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
851 while(E1 < EkinMaxForProd) 851 while(E1 < EkinMaxForProd)
852 { 852 {
853 E1 = std::max(EkinMin, E2); 853 E1 = std::max(EkinMin, E2);
854 E1 = std::min(EkinMaxForProd, E1); 854 E1 = std::min(EkinMaxForProd, E1);
855 std::vector<std::vector<double>*> aMat = 855 std::vector<std::vector<double>*> aMat =
856 aModel->ComputeAdjointCrossSectionVector 856 aModel->ComputeAdjointCrossSectionVectorPerAtomForSecond(E1, Z, A,
857 857 nbin_pro_decade);
858 if(aMat.size() >= 2) 858 if(aMat.size() >= 2)
859 { 859 {
860 std::vector<double>* log_ESecVec = aMat[ 860 std::vector<double>* log_ESecVec = aMat[0];
861 std::vector<double>* log_CSVec = aMat[ 861 std::vector<double>* log_CSVec = aMat[1];
862 G4double log_adjointCS = log_C 862 G4double log_adjointCS = log_CSVec->back();
863 // normalise CSVec such that it becomes 863 // normalise CSVec such that it becomes a probability vector
864 for(std::size_t j = 0; j < log_CSVec->si 864 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
865 { 865 {
866 if(j == 0) 866 if(j == 0)
867 (*log_CSVec)[j] = 0.; 867 (*log_CSVec)[j] = 0.;
868 else 868 else
869 (*log_CSVec)[j] = 869 (*log_CSVec)[j] =
870 std::log(1. - std::exp((*log_CSVec 870 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS) + 1e-50);
871 } 871 }
872 (*log_CSVec)[log_CSVec->size() - 1] = 872 (*log_CSVec)[log_CSVec->size() - 1] =
873 (*log_CSVec)[log_CSVec->size() - 2] - 873 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
874 theCSMatForProdToProjBackwardScattering- 874 theCSMatForProdToProjBackwardScattering->AddData(
875 std::log(E1), log_adjointCS, log_ESecV 875 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
876 } 876 }
877 E1 = E2; 877 E1 = E2;
878 E2 *= dE; 878 E2 *= dE;
879 } 879 }
880 880
881 // Scattered projectile to projectile backwa 881 // Scattered projectile to projectile backward scattering
882 E2 = std::pow(10., G4double(int(std::log10(E 882 E2 = std::pow(10., G4double(int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
883 nbin_pro_decade) / dE; 883 nbin_pro_decade) / dE;
884 E1 = EkinMin; 884 E1 = EkinMin;
885 // Loop checking, 07-Aug-2015, Vladimir Ivan 885 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
886 while(E1 < EkinMaxForScat) 886 while(E1 < EkinMaxForScat)
887 { 887 {
888 E1 = std::max(EkinMin, E2); 888 E1 = std::max(EkinMin, E2);
889 E1 = std::min(EkinMaxForScat, E1); 889 E1 = std::min(EkinMaxForScat, E1);
890 std::vector<std::vector<G4double>*> aMat = 890 std::vector<std::vector<G4double>*> aMat =
891 aModel->ComputeAdjointCrossSectionVector 891 aModel->ComputeAdjointCrossSectionVectorPerAtomForScatProj(
892 E1, Z, A, nbin_pro_decade); 892 E1, Z, A, nbin_pro_decade);
893 if(aMat.size() >= 2) 893 if(aMat.size() >= 2)
894 { 894 {
895 std::vector<G4double>* log_ESecVec = aMa 895 std::vector<G4double>* log_ESecVec = aMat[0];
896 std::vector<G4double>* log_CSVec = aMa 896 std::vector<G4double>* log_CSVec = aMat[1];
897 G4double log_adjointCS = log 897 G4double log_adjointCS = log_CSVec->back();
898 // normalise CSVec such that it becomes 898 // normalise CSVec such that it becomes a probability vector
899 for(std::size_t j = 0; j < log_CSVec->si 899 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
900 { 900 {
901 if(j == 0) 901 if(j == 0)
902 (*log_CSVec)[j] = 0.; 902 (*log_CSVec)[j] = 0.;
903 else 903 else
904 (*log_CSVec)[j] = 904 (*log_CSVec)[j] =
905 std::log(1. - std::exp((*log_CSVec 905 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS) + 1e-50);
906 } 906 }
907 (*log_CSVec)[log_CSVec->size() - 1] = 907 (*log_CSVec)[log_CSVec->size() - 1] =
908 (*log_CSVec)[log_CSVec->size() - 2] - 908 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
909 theCSMatForScatProjToProjBackwardScatter 909 theCSMatForScatProjToProjBackwardScattering->AddData(
910 std::log(E1), log_adjointCS, log_ESecV 910 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
911 } 911 }
912 E1 = E2; 912 E1 = E2;
913 E2 *= dE; 913 E2 *= dE;
914 } 914 }
915 915
916 std::vector<G4AdjointCSMatrix*> res; 916 std::vector<G4AdjointCSMatrix*> res;
917 res.push_back(theCSMatForProdToProjBackwardS 917 res.push_back(theCSMatForProdToProjBackwardScattering);
918 res.push_back(theCSMatForScatProjToProjBackw 918 res.push_back(theCSMatForScatProjToProjBackwardScattering);
919 919
920 return res; 920 return res;
921 } 921 }
922 922
923 ////////////////////////////////////////////// 923 ///////////////////////////////////////////////////////
924 std::vector<G4AdjointCSMatrix*> 924 std::vector<G4AdjointCSMatrix*>
925 G4AdjointCSManager::BuildCrossSectionsModelAnd 925 G4AdjointCSManager::BuildCrossSectionsModelAndMaterial(
926 G4VEmAdjointModel* aModel, G4Material* aMate 926 G4VEmAdjointModel* aModel, G4Material* aMaterial, G4int nbin_pro_decade)
927 { 927 {
928 G4AdjointCSMatrix* theCSMatForProdToProjBack 928 G4AdjointCSMatrix* theCSMatForProdToProjBackwardScattering =
929 new G4AdjointCSMatrix(false); 929 new G4AdjointCSMatrix(false);
930 G4AdjointCSMatrix* theCSMatForScatProjToProj 930 G4AdjointCSMatrix* theCSMatForScatProjToProjBackwardScattering =
931 new G4AdjointCSMatrix(true); 931 new G4AdjointCSMatrix(true);
932 932
933 G4double EkinMin = aModel->GetLowEner 933 G4double EkinMin = aModel->GetLowEnergyLimit();
934 G4double EkinMaxForScat = aModel->GetHighEne 934 G4double EkinMaxForScat = aModel->GetHighEnergyLimit() * 0.999;
935 G4double EkinMaxForProd = aModel->GetHighEne 935 G4double EkinMaxForProd = aModel->GetHighEnergyLimit() * 0.999;
936 if(aModel->GetSecondPartOfSameType()) 936 if(aModel->GetSecondPartOfSameType())
937 EkinMaxForProd /= 2.; 937 EkinMaxForProd /= 2.;
938 938
939 // Product to projectile backward scattering 939 // Product to projectile backward scattering
940 G4double dE = std::pow(10., 1. / nbin_pro_de 940 G4double dE = std::pow(10., 1. / nbin_pro_decade);
941 G4double E2 = 941 G4double E2 =
942 std::pow(10., double(int(std::log10(EkinMi 942 std::pow(10., double(int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
943 nbin_pro_decade) / dE; 943 nbin_pro_decade) / dE;
944 G4double E1 = EkinMin; 944 G4double E1 = EkinMin;
945 // Loop checking, 07-Aug-2015, Vladimir Ivan 945 // Loop checking, 07-Aug-2015, Vladimir Ivanchenko
946 while(E1 < EkinMaxForProd) 946 while(E1 < EkinMaxForProd)
947 { 947 {
948 E1 = std::max(EkinMin, E2); 948 E1 = std::max(EkinMin, E2);
949 E1 = std::min(EkinMaxForProd, E1); 949 E1 = std::min(EkinMaxForProd, E1);
950 std::vector<std::vector<G4double>*> aMat = 950 std::vector<std::vector<G4double>*> aMat =
951 aModel->ComputeAdjointCrossSectionVector 951 aModel->ComputeAdjointCrossSectionVectorPerVolumeForSecond(
952 aMaterial, E1, nbin_pro_decade); 952 aMaterial, E1, nbin_pro_decade);
953 if(aMat.size() >= 2) 953 if(aMat.size() >= 2)
954 { 954 {
955 std::vector<G4double>* log_ESecVec = aMa 955 std::vector<G4double>* log_ESecVec = aMat[0];
956 std::vector<G4double>* log_CSVec = aMa 956 std::vector<G4double>* log_CSVec = aMat[1];
957 G4double log_adjointCS = log 957 G4double log_adjointCS = log_CSVec->back();
958 958
959 // normalise CSVec such that it becomes 959 // normalise CSVec such that it becomes a probability vector
960 for(std::size_t j = 0; j < log_CSVec->si 960 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
961 { 961 {
962 if(j == 0) 962 if(j == 0)
963 (*log_CSVec)[j] = 0.; 963 (*log_CSVec)[j] = 0.;
964 else 964 else
965 (*log_CSVec)[j] = 965 (*log_CSVec)[j] =
966 std::log(1. - std::exp((*log_CSVec 966 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS));
967 } 967 }
968 (*log_CSVec)[log_CSVec->size() - 1] = 968 (*log_CSVec)[log_CSVec->size() - 1] =
969 (*log_CSVec)[log_CSVec->size() - 2] - 969 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
970 theCSMatForProdToProjBackwardScattering- 970 theCSMatForProdToProjBackwardScattering->AddData(
971 std::log(E1), log_adjointCS, log_ESecV 971 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
972 } 972 }
973 973
974 E1 = E2; 974 E1 = E2;
975 E2 *= dE; 975 E2 *= dE;
976 } 976 }
977 977
978 // Scattered projectile to projectile backwa 978 // Scattered projectile to projectile backward scattering
979 E2 = 979 E2 =
980 std::pow(10., G4double(G4int(std::log10(Ek 980 std::pow(10., G4double(G4int(std::log10(EkinMin) * nbin_pro_decade) + 1) /
981 nbin_pro_decade) / 981 nbin_pro_decade) /
982 dE; 982 dE;
983 E1 = EkinMin; 983 E1 = EkinMin;
984 while(E1 < EkinMaxForScat) 984 while(E1 < EkinMaxForScat)
985 { 985 {
986 E1 = std::max(EkinMin, E2); 986 E1 = std::max(EkinMin, E2);
987 E1 = std::min(EkinMaxForScat, E1); 987 E1 = std::min(EkinMaxForScat, E1);
988 std::vector<std::vector<G4double>*> aMat = 988 std::vector<std::vector<G4double>*> aMat =
989 aModel->ComputeAdjointCrossSectionVector 989 aModel->ComputeAdjointCrossSectionVectorPerVolumeForScatProj(
990 aMaterial, E1, nbin_pro_decade); 990 aMaterial, E1, nbin_pro_decade);
991 if(aMat.size() >= 2) 991 if(aMat.size() >= 2)
992 { 992 {
993 std::vector<G4double>* log_ESecVec = aMa 993 std::vector<G4double>* log_ESecVec = aMat[0];
994 std::vector<G4double>* log_CSVec = aMa 994 std::vector<G4double>* log_CSVec = aMat[1];
995 G4double log_adjointCS = log 995 G4double log_adjointCS = log_CSVec->back();
996 996
997 for(std::size_t j = 0; j < log_CSVec->si 997 for(std::size_t j = 0; j < log_CSVec->size(); ++j)
998 { 998 {
999 if(j == 0) 999 if(j == 0)
1000 (*log_CSVec)[j] = 0.; 1000 (*log_CSVec)[j] = 0.;
1001 else 1001 else
1002 (*log_CSVec)[j] = 1002 (*log_CSVec)[j] =
1003 std::log(1. - std::exp((*log_CSVe 1003 std::log(1. - std::exp((*log_CSVec)[j] - log_adjointCS));
1004 } 1004 }
1005 (*log_CSVec)[log_CSVec->size() - 1] = 1005 (*log_CSVec)[log_CSVec->size() - 1] =
1006 (*log_CSVec)[log_CSVec->size() - 2] - 1006 (*log_CSVec)[log_CSVec->size() - 2] - std::log(1000.);
1007 1007
1008 theCSMatForScatProjToProjBackwardScatte 1008 theCSMatForScatProjToProjBackwardScattering->AddData(
1009 std::log(E1), log_adjointCS, log_ESec 1009 std::log(E1), log_adjointCS, log_ESecVec, log_CSVec, 0);
1010 } 1010 }
1011 E1 = E2; 1011 E1 = E2;
1012 E2 *= dE; 1012 E2 *= dE;
1013 } 1013 }
1014 1014
1015 std::vector<G4AdjointCSMatrix*> res; 1015 std::vector<G4AdjointCSMatrix*> res;
1016 res.push_back(theCSMatForProdToProjBackward 1016 res.push_back(theCSMatForProdToProjBackwardScattering);
1017 res.push_back(theCSMatForScatProjToProjBack 1017 res.push_back(theCSMatForScatProjToProjBackwardScattering);
1018 1018
1019 return res; 1019 return res;
1020 } 1020 }
1021 1021
1022 ///////////////////////////////////////////// 1022 ///////////////////////////////////////////////////////
1023 G4ParticleDefinition* G4AdjointCSManager::Get 1023 G4ParticleDefinition* G4AdjointCSManager::GetAdjointParticleEquivalent(
1024 G4ParticleDefinition* theFwdPartDef) 1024 G4ParticleDefinition* theFwdPartDef)
1025 { 1025 {
1026 if(theFwdPartDef->GetParticleName() == "e-" 1026 if(theFwdPartDef->GetParticleName() == "e-")
1027 return G4AdjointElectron::AdjointElectron 1027 return G4AdjointElectron::AdjointElectron();
1028 else if(theFwdPartDef->GetParticleName() == 1028 else if(theFwdPartDef->GetParticleName() == "gamma")
1029 return G4AdjointGamma::AdjointGamma(); 1029 return G4AdjointGamma::AdjointGamma();
1030 else if(theFwdPartDef->GetParticleName() == 1030 else if(theFwdPartDef->GetParticleName() == "proton")
1031 return G4AdjointProton::AdjointProton(); 1031 return G4AdjointProton::AdjointProton();
1032 else if(theFwdPartDef == fFwdIon) 1032 else if(theFwdPartDef == fFwdIon)
1033 return fAdjIon; 1033 return fAdjIon;
1034 return nullptr; 1034 return nullptr;
1035 } 1035 }
1036 1036
1037 ///////////////////////////////////////////// 1037 ///////////////////////////////////////////////////////
1038 G4ParticleDefinition* G4AdjointCSManager::Get 1038 G4ParticleDefinition* G4AdjointCSManager::GetForwardParticleEquivalent(
1039 G4ParticleDefinition* theAdjPartDef) 1039 G4ParticleDefinition* theAdjPartDef)
1040 { 1040 {
1041 if(theAdjPartDef->GetParticleName() == "adj 1041 if(theAdjPartDef->GetParticleName() == "adj_e-")
1042 return G4Electron::Electron(); 1042 return G4Electron::Electron();
1043 else if(theAdjPartDef->GetParticleName() == 1043 else if(theAdjPartDef->GetParticleName() == "adj_gamma")
1044 return G4Gamma::Gamma(); 1044 return G4Gamma::Gamma();
1045 else if(theAdjPartDef->GetParticleName() == 1045 else if(theAdjPartDef->GetParticleName() == "adj_proton")
1046 return G4Proton::Proton(); 1046 return G4Proton::Proton();
1047 else if(theAdjPartDef == fAdjIon) 1047 else if(theAdjPartDef == fAdjIon)
1048 return fFwdIon; 1048 return fFwdIon;
1049 return nullptr; 1049 return nullptr;
1050 } 1050 }
1051 1051
1052 ///////////////////////////////////////////// 1052 ///////////////////////////////////////////////////////
1053 void G4AdjointCSManager::DefineCurrentMateria 1053 void G4AdjointCSManager::DefineCurrentMaterial(
1054 const G4MaterialCutsCouple* couple) 1054 const G4MaterialCutsCouple* couple)
1055 { 1055 {
1056 if(couple != fCurrentCouple) 1056 if(couple != fCurrentCouple)
1057 { 1057 {
1058 fCurrentCouple = const_cast<G4Ma 1058 fCurrentCouple = const_cast<G4MaterialCutsCouple*>(couple);
1059 fCurrentMaterial = const_cast<G4Ma 1059 fCurrentMaterial = const_cast<G4Material*>(couple->GetMaterial());
1060 fCurrentMatIndex = couple->GetInde 1060 fCurrentMatIndex = couple->GetIndex();
1061 fLastCSCorrectionFactor = 1.; 1061 fLastCSCorrectionFactor = 1.;
1062 } 1062 }
1063 } 1063 }
1064 1064
1065 ///////////////////////////////////////////// 1065 ///////////////////////////////////////////////////////
1066 void G4AdjointCSManager::DefineCurrentParticl 1066 void G4AdjointCSManager::DefineCurrentParticle(
1067 const G4ParticleDefinition* aPartDef) 1067 const G4ParticleDefinition* aPartDef)
1068 { 1068 {
1069 << 1069 static G4ParticleDefinition* currentParticleDef = nullptr;
1070 if(aPartDef != fCurrentParticleDef) << 1070
>> 1071 if(aPartDef != currentParticleDef)
1071 { 1072 {
1072 fCurrentParticleDef = const_cast<G4Partic << 1073 currentParticleDef = const_cast<G4ParticleDefinition*>(aPartDef);
1073 fMassRatio = 1.; 1074 fMassRatio = 1.;
1074 if(aPartDef == fAdjIon) 1075 if(aPartDef == fAdjIon)
1075 fMassRatio = proton_mass_c2 / aPartDef- 1076 fMassRatio = proton_mass_c2 / aPartDef->GetPDGMass();
1076 fCurrentParticleIndex = 1000000; 1077 fCurrentParticleIndex = 1000000;
1077 for(std::size_t i = 0; i < fAdjointPartic 1078 for(std::size_t i = 0; i < fAdjointParticlesInAction.size(); ++i)
1078 { 1079 {
1079 if(aPartDef == fAdjointParticlesInActio 1080 if(aPartDef == fAdjointParticlesInAction[i])
1080 fCurrentParticleIndex = i; 1081 fCurrentParticleIndex = i;
1081 } 1082 }
1082 } 1083 }
1083 } 1084 }
1084 1085
1085 ///////////////////////////////////////////// 1086 ////////////////////////////////////////////////////////////////////////////////
1086 G4double G4AdjointCSManager::ComputeAdjointCS 1087 G4double G4AdjointCSManager::ComputeAdjointCS(
1087 G4double aPrimEnergy, G4AdjointCSMatrix* an 1088 G4double aPrimEnergy, G4AdjointCSMatrix* anAdjointCSMatrix, G4double Tcut)
1088 { 1089 {
1089 std::vector<double>* theLogPrimEnergyVector 1090 std::vector<double>* theLogPrimEnergyVector =
1090 anAdjointCSMatrix->GetLogPrimEnergyVector 1091 anAdjointCSMatrix->GetLogPrimEnergyVector();
1091 if(theLogPrimEnergyVector->empty()) 1092 if(theLogPrimEnergyVector->empty())
1092 { 1093 {
1093 G4cout << "No data are contained in the g 1094 G4cout << "No data are contained in the given AdjointCSMatrix!" << G4endl;
1094 return 0.; 1095 return 0.;
1095 } 1096 }
1096 G4double log_Tcut = std::log(Tcut); 1097 G4double log_Tcut = std::log(Tcut);
1097 G4double log_E = std::log(aPrimEnergy); 1098 G4double log_E = std::log(aPrimEnergy);
1098 1099
1099 if(aPrimEnergy <= Tcut || log_E > theLogPri 1100 if(aPrimEnergy <= Tcut || log_E > theLogPrimEnergyVector->back())
1100 return 0.; 1101 return 0.;
1101 1102
1102 G4AdjointInterpolator* theInterpolator = G4 1103 G4AdjointInterpolator* theInterpolator = G4AdjointInterpolator::GetInstance();
1103 1104
1104 std::size_t ind = 1105 std::size_t ind =
1105 theInterpolator->FindPositionForLogVector 1106 theInterpolator->FindPositionForLogVector(log_E, *theLogPrimEnergyVector);
1106 G4double aLogPrimEnergy1, aLogPrimEnergy2; 1107 G4double aLogPrimEnergy1, aLogPrimEnergy2;
1107 G4double aLogCS1, aLogCS2; 1108 G4double aLogCS1, aLogCS2;
1108 G4double log01, log02; 1109 G4double log01, log02;
1109 std::vector<G4double>* aLogSecondEnergyVect 1110 std::vector<G4double>* aLogSecondEnergyVector1 = nullptr;
1110 std::vector<G4double>* aLogSecondEnergyVect 1111 std::vector<G4double>* aLogSecondEnergyVector2 = nullptr;
1111 std::vector<G4double>* aLogProbVector1 1112 std::vector<G4double>* aLogProbVector1 = nullptr;
1112 std::vector<G4double>* aLogProbVector2 1113 std::vector<G4double>* aLogProbVector2 = nullptr;
1113 std::vector<std::size_t>* aLogProbVectorInd 1114 std::vector<std::size_t>* aLogProbVectorIndex1 = nullptr;
1114 std::vector<std::size_t>* aLogProbVectorInd 1115 std::vector<std::size_t>* aLogProbVectorIndex2 = nullptr;
1115 1116
1116 anAdjointCSMatrix->GetData((G4int)ind, aLog 1117 anAdjointCSMatrix->GetData((G4int)ind, aLogPrimEnergy1, aLogCS1, log01,
1117 aLogSecondEnergy 1118 aLogSecondEnergyVector1, aLogProbVector1,
1118 aLogProbVectorIn 1119 aLogProbVectorIndex1);
1119 anAdjointCSMatrix->GetData(G4int(ind + 1), 1120 anAdjointCSMatrix->GetData(G4int(ind + 1), aLogPrimEnergy2, aLogCS2, log02,
1120 aLogSecondEnergy 1121 aLogSecondEnergyVector2, aLogProbVector2,
1121 aLogProbVectorIn 1122 aLogProbVectorIndex2);
1122 if (! (aLogProbVector1 && aLogProbVector2 & 1123 if (! (aLogProbVector1 && aLogProbVector2 &&
1123 aLogSecondEnergyVector1 && aLogS 1124 aLogSecondEnergyVector1 && aLogSecondEnergyVector2)){
1124 return 0.; 1125 return 0.;
1125 } 1126 }
1126 1127
1127 if(anAdjointCSMatrix->IsScatProjToProj()) 1128 if(anAdjointCSMatrix->IsScatProjToProj())
1128 { // case where the Tcut plays a role 1129 { // case where the Tcut plays a role
1129 G4double log_minimum_prob1, log_minimum_p 1130 G4double log_minimum_prob1, log_minimum_prob2;
1130 log_minimum_prob1 = theInterpolator->Inte 1131 log_minimum_prob1 = theInterpolator->InterpolateForLogVector(
1131 log_Tcut, *aLogSecondEnergyVector1, *aL 1132 log_Tcut, *aLogSecondEnergyVector1, *aLogProbVector1);
1132 log_minimum_prob2 = theInterpolator->Inte 1133 log_minimum_prob2 = theInterpolator->InterpolateForLogVector(
1133 log_Tcut, *aLogSecondEnergyVector2, *aL 1134 log_Tcut, *aLogSecondEnergyVector2, *aLogProbVector2);
1134 aLogCS1 += log_minimum_prob1; 1135 aLogCS1 += log_minimum_prob1;
1135 aLogCS2 += log_minimum_prob2; 1136 aLogCS2 += log_minimum_prob2;
1136 } 1137 }
1137 1138
1138 G4double log_adjointCS = theInterpolator->L 1139 G4double log_adjointCS = theInterpolator->LinearInterpolation(
1139 log_E, aLogPrimEnergy1, aLogPrimEnergy2, 1140 log_E, aLogPrimEnergy1, aLogPrimEnergy2, aLogCS1, aLogCS2);
1140 return std::exp(log_adjointCS); 1141 return std::exp(log_adjointCS);
1141 } 1142 }
1142 1143