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