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