| 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 // neutron_hp -- source file 26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 09-May-06 fix in Sample by T. Koi 30 // 09-May-06 fix in Sample by T. Koi
31 // 080318 Fix Compilation warnings - gcc-4.3.0 31 // 080318 Fix Compilation warnings - gcc-4.3.0 by T. Koi
32 // (This fix has a real effect to the c << 32 // (This fix has a real effect to the code.)
33 // 080409 Fix div0 error with G4FPE by T. Koi 33 // 080409 Fix div0 error with G4FPE by T. Koi
34 // 080612 Fix contribution from Benoit Pirard 34 // 080612 Fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #1
35 // 080714 Limiting the sum of energy of second 35 // 080714 Limiting the sum of energy of secondary particles by T. Koi
36 // 080801 Fix div0 error wiht G4FPE and memory 36 // 080801 Fix div0 error wiht G4FPE and memory leak by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4Nucle 37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 38 //
39 // P. Arce, June-2014 Conversion neutron_hp to 39 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
40 // 40 //
41 // June-2019 - E. Mendoza --> redefinition of << 41 // June-2019 - E. Mendoza --> redefinition of the residual mass to consider incident particles different than neutrons.
42 // different than neutrons. <<
43 // <<
44 // V. Ivanchenko, July-2023 Basic revision of <<
45 // <<
46 42
47 #include "G4ParticleHPContAngularPar.hh" 43 #include "G4ParticleHPContAngularPar.hh"
48 << 44 #include "G4PhysicalConstants.hh"
49 #include "G4ParticleDefinition.hh" << 45 #include "G4SystemOfUnits.hh"
50 #include "G4Alpha.hh" << 46 #include "G4ParticleHPLegendreStore.hh"
51 #include "G4Deuteron.hh" <<
52 #include "G4Electron.hh" <<
53 #include "G4Gamma.hh" 47 #include "G4Gamma.hh"
54 #include "G4He3.hh" << 48 #include "G4Electron.hh"
55 #include "G4IonTable.hh" << 49 #include "G4Positron.hh"
56 #include "G4Neutron.hh" 50 #include "G4Neutron.hh"
>> 51 #include "G4Proton.hh"
>> 52 #include "G4Deuteron.hh"
>> 53 #include "G4Triton.hh"
>> 54 #include "G4He3.hh"
>> 55 #include "G4Alpha.hh"
>> 56 #include "G4ParticleHPVector.hh"
57 #include "G4NucleiProperties.hh" 57 #include "G4NucleiProperties.hh"
58 #include "G4ParticleHPKallbachMannSyst.hh" 58 #include "G4ParticleHPKallbachMannSyst.hh"
59 #include "G4ParticleHPLegendreStore.hh" << 59 #include "G4IonTable.hh"
60 #include "G4ParticleHPManager.hh" 60 #include "G4ParticleHPManager.hh"
61 #include "G4ParticleHPVector.hh" <<
62 #include "G4PhysicalConstants.hh" <<
63 #include "G4Positron.hh" <<
64 #include "G4Proton.hh" <<
65 #include "G4SystemOfUnits.hh" <<
66 #include "G4Triton.hh" <<
67 <<
68 #include <set> 61 #include <set>
69 #include <vector> 62 #include <vector>
>> 63
>> 64 G4ParticleHPContAngularPar::G4ParticleHPContAngularPar(G4ParticleDefinition* projectile)
>> 65 {
>> 66 theAngular = 0;
>> 67 if (fCache.Get() == 0) cacheInit();
>> 68 fCache.Get()->currentMeanEnergy = -2;
>> 69 fCache.Get()->fresh = true;
>> 70 adjustResult = true;
>> 71 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState() )
>> 72 adjustResult = false;
>> 73
>> 74 theMinEner = DBL_MAX;
>> 75 theMaxEner = -DBL_MAX;
>> 76 theProjectile = projectile;
70 77
71 G4ParticleHPContAngularPar::G4ParticleHPContAn << 78 theEnergy = 0.0;
72 { << 79 nEnergies = 0;
73 theProjectile = (nullptr == p) ? G4Neutron:: << 80 nDiscreteEnergies = 0;
74 toBeCached v; << 81 nAngularParameters = 0;
75 fCache.Put(v); <<
76 if (G4ParticleHPManager::GetInstance()->GetD <<
77 } <<
78 <<
79 G4ParticleHPContAngularPar::G4ParticleHPContAn <<
80 { <<
81 theEnergy = val.theEnergy; <<
82 nEnergies = val.nEnergies; <<
83 nDiscreteEnergies = val.nDiscreteEnergies; <<
84 nAngularParameters = val.nAngularParameters; <<
85 theProjectile = val.theProjectile; <<
86 theManager = val.theManager; <<
87 theInt = val.theInt; <<
88 adjustResult = val.adjustResult; <<
89 theMinEner = val.theMinEner; <<
90 theMaxEner = val.theMaxEner; <<
91 theEnergiesTransformed = val.theEnergiesTran <<
92 theDiscreteEnergies = val.theDiscreteEnergie <<
93 theDiscreteEnergiesOwn = val.theDiscreteEner <<
94 toBeCached v; <<
95 fCache.Put(v); <<
96 const std::size_t esize = nEnergies > 0 ? nE <<
97 theAngular = new G4ParticleHPList[esize]; <<
98 for (G4int ie = 0; ie < nEnergies; ++ie) { <<
99 theAngular[ie].SetLabel(val.theAngular[ie] <<
100 for (G4int ip = 0; ip < nAngularParameters <<
101 theAngular[ie].SetValue(ip, val.theAngul <<
102 } <<
103 } <<
104 } 82 }
105 83
106 G4ParticleHPContAngularPar::~G4ParticleHPContA <<
107 { <<
108 delete[] theAngular; <<
109 } <<
110 84
111 void G4ParticleHPContAngularPar::Init(std::ist << 85 void
112 { << 86 G4ParticleHPContAngularPar::Init(std::istream& aDataFile, G4ParticleDefinition* projectile)
>> 87 {
113 adjustResult = true; 88 adjustResult = true;
114 if (G4ParticleHPManager::GetInstance()->GetD << 89 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState() )
>> 90 adjustResult = false;
115 91
116 theProjectile = (nullptr == p) ? G4Neutron:: << 92 theProjectile = projectile;
117 93
118 aDataFile >> theEnergy >> nEnergies >> nDisc 94 aDataFile >> theEnergy >> nEnergies >> nDiscreteEnergies >> nAngularParameters;
119 theEnergy *= eV; 95 theEnergy *= eV;
120 const std::size_t esize = nEnergies > 0 ? nE << 96 theAngular = new G4ParticleHPList [nEnergies];
121 theAngular = new G4ParticleHPList[esize]; <<
122 G4double sEnergy; 97 G4double sEnergy;
123 for (G4int i = 0; i < nEnergies; ++i) { << 98 for (G4int i = 0; i < nEnergies; i++) {
124 aDataFile >> sEnergy; 99 aDataFile >> sEnergy;
125 sEnergy *= eV; 100 sEnergy *= eV;
126 theAngular[i].SetLabel(sEnergy); 101 theAngular[i].SetLabel(sEnergy);
127 theAngular[i].Init(aDataFile, nAngularPara 102 theAngular[i].Init(aDataFile, nAngularParameters, 1.);
128 theMinEner = std::min(theMinEner, sEnergy) << 103 theMinEner = std::min(theMinEner,sEnergy);
129 theMaxEner = std::max(theMaxEner, sEnergy) << 104 theMaxEner = std::max(theMaxEner,sEnergy);
130 } 105 }
131 } 106 }
132 107
133 G4ReactionProduct* G4ParticleHPContAngularPar: << 108
134 << 109 G4ReactionProduct*
135 << 110 G4ParticleHPContAngularPar::Sample(G4double anEnergy, G4double massCode,
>> 111 G4double /*targetMass*/,
>> 112 G4int angularRep, G4int /*interpolE*/ )
136 { 113 {
137 // The following line is needed because it m 114 // The following line is needed because it may change between runs by UI command
138 adjustResult = true; << 115 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState() )
139 if (G4ParticleHPManager::GetInstance()->GetD << 116 adjustResult = false;
140 << 117
141 auto result = new G4ReactionProduct; << 118 if (fCache.Get() == 0 ) cacheInit();
142 auto Z = static_cast<G4int>(massCode / 1000) << 119 G4ReactionProduct * result = new G4ReactionProduct;
143 auto A = static_cast<G4int>(massCode - 1000 << 120 G4int Z = static_cast<G4int>(massCode/1000);
>> 121 G4int A = static_cast<G4int>(massCode-1000*Z);
144 if (massCode == 0) { 122 if (massCode == 0) {
145 result->SetDefinition(G4Gamma::Gamma()); 123 result->SetDefinition(G4Gamma::Gamma());
146 } << 124 } else if (A == 0) {
147 else if (A == 0) { << 125 result->SetDefinition(G4Electron::Electron());
148 result->SetDefinition(G4Electron::Electron <<
149 if (Z == 1) result->SetDefinition(G4Positr 126 if (Z == 1) result->SetDefinition(G4Positron::Positron());
150 } << 127 } else if (A == 1) {
151 else if (A == 1) { <<
152 result->SetDefinition(G4Neutron::Neutron() 128 result->SetDefinition(G4Neutron::Neutron());
153 if (Z == 1) result->SetDefinition(G4Proton 129 if (Z == 1) result->SetDefinition(G4Proton::Proton());
154 } << 130 } else if (A == 2) {
155 else if (A == 2) { << 131 result->SetDefinition(G4Deuteron::Deuteron());
156 result->SetDefinition(G4Deuteron::Deuteron << 132 } else if (A == 3) {
157 } << 133 result->SetDefinition(G4Triton::Triton());
158 else if (A == 3) { << 134 if(Z == 2) result->SetDefinition(G4He3::He3());
159 result->SetDefinition(G4Triton::Triton()); << 135 } else if (A == 4) {
160 if (Z == 2) result->SetDefinition(G4He3::H <<
161 } <<
162 else if (A == 4) { <<
163 result->SetDefinition(G4Alpha::Alpha()); 136 result->SetDefinition(G4Alpha::Alpha());
164 if (Z != 2) << 137 if (Z != 2) throw G4HadronicException(__FILE__, __LINE__,
165 throw G4HadronicException(__FILE__, __LI << 138 "G4ParticleHPContAngularPar: Unknown ion case 1");
166 "G4ParticleHPC << 139 } else {
167 } << 140 result->SetDefinition(G4IonTable::GetIonTable()->GetIon(Z,A,0) );
168 else { <<
169 result->SetDefinition(G4IonTable::GetIonTa <<
170 } 141 }
171 142
172 G4int i(0); 143 G4int i(0);
173 G4int it(0); 144 G4int it(0);
174 G4double fsEnergy(0); 145 G4double fsEnergy(0);
175 G4double cosTh(0); 146 G4double cosTh(0);
176 /* << 147
177 G4cout << "G4ParticleHPContAngularPar::Sampl <<
178 << " angularRep=" << angularRep << " <<
179 << " Ne=" << nEnergies << G4endl; <<
180 */ <<
181 if (angularRep == 1) { 148 if (angularRep == 1) {
182 if (nDiscreteEnergies != 0) { 149 if (nDiscreteEnergies != 0) {
183 // 1st check remaining_energy << 150 // 1st check remaining_energy
184 // if this is the first set it. (How?) 151 // if this is the first set it. (How?)
185 if (fCache.Get().fresh) { << 152 if (fCache.Get()->fresh == true) {
186 // Discrete Lines, larger energies com << 153 // Discrete Lines, larger energies come first
187 // Continues Emssions, low to high << 154 // Continues Emssions, low to high LAST
188 fCache.Get().remaining_energy = << 155 fCache.Get()->remaining_energy = std::max (theAngular[0].GetLabel(),
189 std::max(theAngular[0].GetLabel(), t << 156 theAngular[nEnergies-1].GetLabel() );
190 fCache.Get().fresh = false; << 157 fCache.Get()->fresh = false;
191 } 158 }
192 159
193 // Cheating for small remaining_energy << 160 // Cheating for small remaining_energy
194 // Temporary solution 161 // Temporary solution
195 if (nDiscreteEnergies == nEnergies) { 162 if (nDiscreteEnergies == nEnergies) {
196 fCache.Get().remaining_energy = << 163 fCache.Get()->remaining_energy = std::max(fCache.Get()->remaining_energy,
197 std::max(fCache.Get().remaining_ener << 164 theAngular[nDiscreteEnergies-1].GetLabel() ); //Minimum Line
198 theAngular[nDiscreteEnergie << 165 } else {
199 } << 166 G4double cont_min = 0.0;
200 else { << 167 for (G4int j = nDiscreteEnergies; j < nEnergies; j++) {
201 G4double cont_min = 0.0; << 168 cont_min = theAngular[j].GetLabel();
202 for (G4int j = nDiscreteEnergies; j < << 169 if (theAngular[j].GetValue(0) != 0.0 ) break;
203 cont_min = theAngular[j].GetLabel(); << 170 }
204 if (theAngular[j].GetValue(0) != 0.0 << 171 fCache.Get()->remaining_energy =
205 } << 172 std::max (fCache.Get()->remaining_energy,
206 fCache.Get().remaining_energy = std::m << 173 std::min(theAngular[nDiscreteEnergies-1].GetLabel(), cont_min) ); //Minimum Line or grid
207 fCache.Get().remaining_energy, std:: <<
208 <<
209 } 174 }
210 175
211 G4double random = G4UniformRand(); 176 G4double random = G4UniformRand();
212 auto running = new G4double[nEnergies + << 177 G4double* running = new G4double[nEnergies+1];
213 running[0] = 0.0; 178 running[0] = 0.0;
214 179
215 G4double delta; 180 G4double delta;
216 for (G4int j = 0; j < nDiscreteEnergies; << 181 for (G4int j = 0; j < nDiscreteEnergies; j++) {
217 delta = 0.0; 182 delta = 0.0;
218 if (theAngular[j].GetLabel() <= fCache << 183 if (theAngular[j].GetLabel() <= fCache.Get()->remaining_energy)
219 delta = theAngular[j].GetValue(0); 184 delta = theAngular[j].GetValue(0);
220 running[j + 1] = running[j] + delta; << 185 running[j+1] = running[j] + delta;
221 } 186 }
222 187
223 G4double tot_prob_DIS = std::max(running << 188 G4double tot_prob_DIS = std::max( running[nDiscreteEnergies], 0.0 );
224 189
225 G4double delta1; 190 G4double delta1;
226 for (G4int j = nDiscreteEnergies; j < nE << 191 for (G4int j = nDiscreteEnergies; j < nEnergies; j++) {
227 delta1 = 0.0; 192 delta1 = 0.0;
228 G4double e_low = 0.0; 193 G4double e_low = 0.0;
229 G4double e_high = 0.0; 194 G4double e_high = 0.0;
230 if (theAngular[j].GetLabel() <= fCache << 195 if (theAngular[j].GetLabel() <= fCache.Get()->remaining_energy)
231 delta1 = theAngular[j].GetValue(0); 196 delta1 = theAngular[j].GetValue(0);
232 197
233 // To calculate Prob. e_low and e_high << 198 // To calculate Prob. e_low and e_high should be in eV
234 // There are two cases: 199 // There are two cases:
235 // 1: theAngular[nDiscreteEnergies].Ge 200 // 1: theAngular[nDiscreteEnergies].GetLabel() != 0.0
236 // delta1 should be used between j- << 201 // delta1 should be used between j-1 and j
237 // At j = nDiscreteEnergies (the fi << 202 // At j = nDiscreteEnergies (the first) e_low should be set explicitly
238 if (theAngular[j].GetLabel() != 0) { 203 if (theAngular[j].GetLabel() != 0) {
239 if (j == nDiscreteEnergies) { 204 if (j == nDiscreteEnergies) {
240 e_low = 0.0 / eV; << 205 e_low = 0.0/eV;
241 } << 206 } else {
242 else { << 207 if ( j < 1 ) j = 1; // Protection against evaluation of arrays at index j-1
243 if (j < 1) j = 1; // Protection a << 208 e_low = theAngular[j-1].GetLabel()/eV;
244 e_low = theAngular[j - 1].GetLabel <<
245 } 209 }
246 e_high = theAngular[j].GetLabel() / << 210 e_high = theAngular[j].GetLabel()/eV;
247 } 211 }
248 212
249 // 2: theAngular[nDiscreteEnergies].Ge 213 // 2: theAngular[nDiscreteEnergies].GetLabel() == 0.0
250 // delta1 should be used between j << 214 // delta1 should be used between j and j+1
251 if (theAngular[j].GetLabel() == 0.0) { 215 if (theAngular[j].GetLabel() == 0.0) {
252 e_low = theAngular[j].GetLabel() / e << 216 e_low = theAngular[j].GetLabel()/eV;
253 if (j != nEnergies - 1) { << 217 if (j != nEnergies-1) {
254 e_high = theAngular[j + 1].GetLabe << 218 e_high = theAngular[j+1].GetLabel()/eV;
255 } << 219 } else {
256 else { << 220 e_high = theAngular[j].GetLabel()/eV;
257 e_high = theAngular[j].GetLabel() << 221 if (theAngular[j].GetValue(0) != 0.0) {
>> 222 throw G4HadronicException(__FILE__, __LINE__,
>> 223 "G4ParticleHPContAngularPar: Unexpected non zero value of theAngular[nEnergies-1].GetValue(0)");
>> 224 }
258 } 225 }
259 } 226 }
260 227
261 running[j + 1] = running[j] + ((e_high << 228 running[j+1] = running[j] + ( ( e_high - e_low ) * delta1);
262 } 229 }
263 G4double tot_prob_CON = std::max(running << 230 G4double tot_prob_CON = std::max( running[ nEnergies ] - running[ nDiscreteEnergies ], 0.0 );
264 231
265 // Give up in the pathological case of n << 232 // Give up in the pathological case of null probabilities
266 if (tot_prob_DIS == 0.0 && tot_prob_CON << 233 if ( tot_prob_DIS == 0.0 && tot_prob_CON == 0.0 ) return result;
267 delete[] running; << 234
268 return result; << 235 // Normalize random
269 } <<
270 // Normalize random <<
271 random *= (tot_prob_DIS + tot_prob_CON); 236 random *= (tot_prob_DIS + tot_prob_CON);
272 // 2nd Judge Discrete or not 237 // 2nd Judge Discrete or not
273 238
274 // This should be relatively close to 1 << 239 // This should be relatively close to 1 For safty
275 if (random <= (tot_prob_DIS / (tot_prob_ << 240 if (random <= (tot_prob_DIS/(tot_prob_DIS + tot_prob_CON) ) || nDiscreteEnergies == nEnergies) {
276 || nDiscreteEnergies == nEnergies) << 241 // Discrete Emission
277 { << 242 for (G4int j = 0; j < nDiscreteEnergies; j++) {
278 // Discrete Emission <<
279 for (G4int j = 0; j < nDiscreteEnergie <<
280 // Here we should use i+1 243 // Here we should use i+1
281 if (random < running[j + 1]) { << 244 if (random < running[ j+1 ] ) {
282 it = j; << 245 it = j;
283 break; 246 break;
284 } 247 }
285 } 248 }
286 fsEnergy = theAngular[it].GetLabel(); << 249 fsEnergy = theAngular[ it ].GetLabel();
287 250
288 G4ParticleHPLegendreStore theStore(1); 251 G4ParticleHPLegendreStore theStore(1);
289 theStore.Init(0, fsEnergy, nAngularPar << 252 theStore.Init(0,fsEnergy,nAngularParameters);
290 for (G4int j = 0; j < nAngularParamete << 253 for (G4int j = 0; j < nAngularParameters; j++) {
291 theStore.SetCoeff(0, j, theAngular[i << 254 theStore.SetCoeff(0,j,theAngular[it].GetValue(j));
292 } 255 }
293 // use it to sample. 256 // use it to sample.
294 cosTh = theStore.SampleMax(fsEnergy); 257 cosTh = theStore.SampleMax(fsEnergy);
295 // Done 258 // Done
296 } << 259
297 else { << 260 } else {
298 // Continuous emission 261 // Continuous emission
299 for (G4int j = nDiscreteEnergies; j < << 262 for (G4int j = nDiscreteEnergies; j < nEnergies; j++) {
300 // Here we should use i 263 // Here we should use i
301 if (random < running[j]) { << 264 if (random < running[ j ] ) {
302 it = j; << 265 it = j;
303 break; 266 break;
304 } 267 }
305 } 268 }
306 269
307 if (it < 1) it = 1; // Protection aga << 270 if ( it < 1 ) it = 1; // Protection against evaluation of arrays at index it-1
308 271
309 G4double x1 = running[it - 1]; << 272 G4double x1 = running[it-1];
310 G4double x2 = running[it]; 273 G4double x2 = running[it];
311 274
312 G4double y1 = 0.0; 275 G4double y1 = 0.0;
313 if (it != nDiscreteEnergies) y1 = theA << 276 if (it != nDiscreteEnergies) y1 = theAngular[it-1].GetLabel();
314 G4double y2 = theAngular[it].GetLabel( 277 G4double y2 = theAngular[it].GetLabel();
315 278
316 fsEnergy = theInt.Interpolate(theManag << 279 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
>> 280 random,x1,x2,y1,y2);
317 281
318 G4ParticleHPLegendreStore theStore(2); 282 G4ParticleHPLegendreStore theStore(2);
319 theStore.Init(0, y1, nAngularParameter 283 theStore.Init(0, y1, nAngularParameters);
320 theStore.Init(1, y2, nAngularParameter 284 theStore.Init(1, y2, nAngularParameters);
321 theStore.SetManager(theManager); 285 theStore.SetManager(theManager);
322 G4int itt; 286 G4int itt;
323 for (G4int j = 0; j < nAngularParamete << 287 for (G4int j = 0; j < nAngularParameters; j++) {
324 itt = it; 288 itt = it;
325 if (it == nDiscreteEnergies) itt = i << 289 if (it == nDiscreteEnergies) itt = it+1;
326 // "This case "it-1" has data for Di << 290 // "This case "it-1" has data for Discrete, so we will use an extrpolated values it and it+1
327 // it+1 << 291 theStore.SetCoeff(0, j, theAngular[itt-1].GetValue(j) );
328 theStore.SetCoeff(0, j, theAngular[i << 292 theStore.SetCoeff(1, j, theAngular[itt].GetValue(j) );
329 theStore.SetCoeff(1, j, theAngular[i <<
330 } 293 }
331 // use it to sample. 294 // use it to sample.
332 cosTh = theStore.SampleMax(fsEnergy); 295 cosTh = theStore.SampleMax(fsEnergy);
333 296
334 // Done << 297 //Done
335 } 298 }
336 299
337 // The remaining energy needs to be lowe 300 // The remaining energy needs to be lowered by the photon energy in *any* case.
338 // Otherwise additional photons with too 301 // Otherwise additional photons with too high energy will be produced - therefore the
339 // adjustResult condition has been remov << 302 // adjustResult condition has been removed
340 fCache.Get().remaining_energy -= fsEnerg << 303 fCache.Get()->remaining_energy -= fsEnergy;
341 delete[] running; 304 delete[] running;
342 305
343 // end (nDiscreteEnergies != 0) branch << 306 // end (nDiscreteEnergies != 0) branch
344 } << 307
345 else { << 308 } else {
346 // Only continue, TK will clean up << 309 // Only continue, TK will clean up
347 if (fCache.Get().fresh) { << 310 if (fCache.Get()->fresh == true) {
348 fCache.Get().remaining_energy = theAng << 311 fCache.Get()->remaining_energy = theAngular[ nEnergies-1 ].GetLabel();
349 fCache.Get().fresh = false; << 312 fCache.Get()->fresh = false;
350 } 313 }
351 314
352 G4double random = G4UniformRand(); 315 G4double random = G4UniformRand();
353 auto running = new G4double[nEnergies]; << 316 G4double* running = new G4double[nEnergies];
354 running[0] = 0; 317 running[0] = 0;
355 G4double weighted = 0; 318 G4double weighted = 0;
356 for (i = 1; i < nEnergies; i++) { 319 for (i = 1; i < nEnergies; i++) {
357 running[i] = running[i - 1]; << 320 running[i]=running[i-1];
358 if (fCache.Get().remaining_energy >= t << 321 if (fCache.Get()->remaining_energy >= theAngular[i].GetLabel() ) {
359 running[i] += theInt.GetBinIntegral( << 322 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
360 theManager.GetScheme(i - 1), theAn << 323 theAngular[i-1].GetLabel(),
361 theAngular[i - 1].GetValue(0), the << 324 theAngular[i].GetLabel(),
362 weighted += theInt.GetWeightedBinInt << 325 theAngular[i-1].GetValue(0),
363 theManager.GetScheme(i - 1), theAn << 326 theAngular[i].GetValue(0) );
364 theAngular[i - 1].GetValue(0), the << 327 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
>> 328 theAngular[i-1].GetLabel(),
>> 329 theAngular[i].GetLabel(),
>> 330 theAngular[i-1].GetValue(0),
>> 331 theAngular[i].GetValue(0));
365 } 332 }
366 } 333 }
367 334
368 // Cache the mean energy in this distrib 335 // Cache the mean energy in this distribution
369 if (nEnergies == 1 || running[nEnergies << 336 if (nEnergies == 1 || running[nEnergies-1] == 0) {
370 fCache.Get().currentMeanEnergy = 0.0; << 337 fCache.Get()->currentMeanEnergy = 0.0;
371 } << 338 } else {
372 else { << 339 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
373 fCache.Get().currentMeanEnergy = weigh << 340 }
374 } << 341
375 << 342 if (nEnergies == 1) it = 0;
376 if (nEnergies == 1) it = 0; << 343 if (running[nEnergies-1] != 0) {
377 if (running[nEnergies - 1] != 0) { <<
378 for (i = 1; i < nEnergies; i++) { 344 for (i = 1; i < nEnergies; i++) {
379 it = i; 345 it = i;
380 if (random < running[i] / running[nE << 346 if (random < running [i]/running[nEnergies-1] ) break;
381 } 347 }
382 } 348 }
383 349
384 if (running[nEnergies - 1] == 0) it = 0; << 350 if (running[nEnergies-1] == 0) it = 0;
385 if (it < nDiscreteEnergies || it == 0) { 351 if (it < nDiscreteEnergies || it == 0) {
386 if (it == 0) { 352 if (it == 0) {
387 fsEnergy = theAngular[it].GetLabel() 353 fsEnergy = theAngular[it].GetLabel();
388 G4ParticleHPLegendreStore theStore(1 354 G4ParticleHPLegendreStore theStore(1);
389 theStore.Init(0, fsEnergy, nAngularP << 355 theStore.Init(0,fsEnergy,nAngularParameters);
390 for (i = 0; i < nAngularParameters; 356 for (i = 0; i < nAngularParameters; i++) {
391 theStore.SetCoeff(0, i, theAngular << 357 theStore.SetCoeff(0, i, theAngular[it].GetValue(i) );
392 } 358 }
393 // use it to sample. 359 // use it to sample.
394 cosTh = theStore.SampleMax(fsEnergy) 360 cosTh = theStore.SampleMax(fsEnergy);
395 } << 361
396 else { << 362 } else {
397 G4double e1, e2; 363 G4double e1, e2;
398 e1 = theAngular[it - 1].GetLabel(); << 364 e1 = theAngular[it-1].GetLabel();
399 e2 = theAngular[it].GetLabel(); 365 e2 = theAngular[it].GetLabel();
400 fsEnergy = theInt.Interpolate(theMan << 366 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
401 runnin << 367 random,
402 runnin << 368 running[it-1]/running[nEnergies-1],
>> 369 running[it]/running[nEnergies-1],
>> 370 e1, e2);
403 // fill a Legendrestore 371 // fill a Legendrestore
404 G4ParticleHPLegendreStore theStore(2 372 G4ParticleHPLegendreStore theStore(2);
405 theStore.Init(0, e1, nAngularParamet << 373 theStore.Init(0,e1,nAngularParameters);
406 theStore.Init(1, e2, nAngularParamet << 374 theStore.Init(1,e2,nAngularParameters);
407 for (i = 0; i < nAngularParameters; 375 for (i = 0; i < nAngularParameters; i++) {
408 theStore.SetCoeff(0, i, theAngular << 376 theStore.SetCoeff(0, i, theAngular[it-1].GetValue(i) );
409 theStore.SetCoeff(1, i, theAngular << 377 theStore.SetCoeff(1, i, theAngular[it].GetValue(i) );
410 } 378 }
411 // use it to sample. 379 // use it to sample.
412 theStore.SetManager(theManager); 380 theStore.SetManager(theManager);
413 cosTh = theStore.SampleMax(fsEnergy) 381 cosTh = theStore.SampleMax(fsEnergy);
414 } 382 }
415 } << 383
416 else { // continuum contribution << 384 } else { // continuum contribution
417 G4double x1 = running[it - 1] / runnin << 385 G4double x1 = running[it-1]/running[nEnergies-1];
418 G4double x2 = running[it] / running[nE << 386 G4double x2 = running[it]/running[nEnergies-1];
419 G4double y1 = theAngular[it - 1].GetLa << 387 G4double y1 = theAngular[it-1].GetLabel();
420 G4double y2 = theAngular[it].GetLabel( 388 G4double y2 = theAngular[it].GetLabel();
421 fsEnergy = theInt.Interpolate(theManag << 389 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
>> 390 random,x1,x2,y1,y2);
422 G4ParticleHPLegendreStore theStore(2); 391 G4ParticleHPLegendreStore theStore(2);
423 theStore.Init(0, y1, nAngularParameter << 392 theStore.Init(0,y1,nAngularParameters);
424 theStore.Init(1, y2, nAngularParameter << 393 theStore.Init(1,y2,nAngularParameters);
425 theStore.SetManager(theManager); 394 theStore.SetManager(theManager);
426 for (i = 0; i < nAngularParameters; i+ 395 for (i = 0; i < nAngularParameters; i++) {
427 theStore.SetCoeff(0, i, theAngular[i << 396 theStore.SetCoeff(0, i, theAngular[it-1].GetValue(i));
428 theStore.SetCoeff(1, i, theAngular[i 397 theStore.SetCoeff(1, i, theAngular[it].GetValue(i));
429 } 398 }
430 // use it to sample. 399 // use it to sample.
431 cosTh = theStore.SampleMax(fsEnergy); 400 cosTh = theStore.SampleMax(fsEnergy);
432 } 401 }
433 delete[] running; << 402 delete [] running;
434 403
435 // The remaining energy needs to be lowe 404 // The remaining energy needs to be lowered by the photon energy in
436 // *any* case. Otherwise additional pho 405 // *any* case. Otherwise additional photons with too much energy will be
437 // produced - therefore the adjustResul 406 // produced - therefore the adjustResult condition has been removed
438 407
439 fCache.Get().remaining_energy -= fsEnerg << 408 fCache.Get()->remaining_energy -= fsEnergy;
440 // end if (nDiscreteEnergies != 0) 409 // end if (nDiscreteEnergies != 0)
441 } 410 }
442 // end of (angularRep == 1) branch 411 // end of (angularRep == 1) branch
443 } << 412
444 else if (angularRep == 2) { << 413 } else if (angularRep == 2) {
445 // first get the energy (already the right 414 // first get the energy (already the right for this incoming energy)
446 G4int j; 415 G4int j;
447 auto running = new G4double[nEnergies]; << 416 G4double* running = new G4double[nEnergies];
448 running[0] = 0; 417 running[0] = 0;
449 G4double weighted = 0; 418 G4double weighted = 0;
450 for (j = 1; j < nEnergies; ++j) { << 419 for (j = 1; j < nEnergies; j++) {
451 if (j != 0) running[j] = running[j - 1]; << 420 if (j != 0) running[j] = running[j-1];
452 running[j] += theInt.GetBinIntegral(theM << 421 running[j] += theInt.GetBinIntegral(theManager.GetScheme(j-1),
453 theA << 422 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
454 theA << 423 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
455 weighted += theInt.GetWeightedBinIntegra << 424 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(j-1),
456 theManager.GetScheme(j - 1), theAngula << 425 theAngular[j-1].GetLabel(), theAngular[j].GetLabel(),
457 theAngular[j - 1].GetValue(0), theAngu << 426 theAngular[j-1].GetValue(0), theAngular[j].GetValue(0));
458 } 427 }
459 428
460 // Cache the mean energy in this distribut 429 // Cache the mean energy in this distribution
461 if (nEnergies == 1) 430 if (nEnergies == 1)
462 fCache.Get().currentMeanEnergy = 0.0; << 431 fCache.Get()->currentMeanEnergy = 0.0;
463 else 432 else
464 fCache.Get().currentMeanEnergy = weighte << 433 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
465 << 434
466 G4int itt(0); 435 G4int itt(0);
467 G4double randkal = G4UniformRand(); 436 G4double randkal = G4UniformRand();
468 for (j = 1; j < nEnergies; ++j) { << 437 for (j = 1; j < nEnergies; j++) {
469 itt = j; 438 itt = j;
470 if (randkal*running[nEnergies - 1] < run << 439 if (randkal < running[j]/running[nEnergies-1] ) break;
471 } 440 }
472 << 441
473 // Interpolate the secondary energy 442 // Interpolate the secondary energy
474 G4double x, x1, x2, y1, y2; 443 G4double x, x1, x2, y1, y2;
475 if (itt == 0) itt = 1; 444 if (itt == 0) itt = 1;
476 x = randkal * running[nEnergies - 1]; << 445 x = randkal*running[nEnergies-1];
477 x1 = running[itt - 1]; << 446 x1 = running[itt-1];
478 x2 = running[itt]; 447 x2 = running[itt];
479 G4double compoundFraction; 448 G4double compoundFraction;
480 // interpolate energy 449 // interpolate energy
481 y1 = theAngular[itt - 1].GetLabel(); << 450 y1 = theAngular[itt-1].GetLabel();
482 y2 = theAngular[itt].GetLabel(); 451 y2 = theAngular[itt].GetLabel();
483 fsEnergy = theInt.Interpolate(theManager.G << 452 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(itt-1),
>> 453 x, x1, x2, y1, y2);
484 454
485 // For theta, interpolate the compoundFrac 455 // For theta, interpolate the compoundFractions
486 G4double cLow = theAngular[itt - 1].GetVal << 456 G4double cLow = theAngular[itt-1].GetValue(1);
487 G4double cHigh = theAngular[itt].GetValue( 457 G4double cHigh = theAngular[itt].GetValue(1);
488 compoundFraction = theInt.Interpolate(theM << 458 compoundFraction = theInt.Interpolate(theManager.GetScheme(itt),
489 << 459 fsEnergy, y1, y2, cLow, cHigh);
490 if (compoundFraction > 1.0) <<
491 compoundFraction = 1.0; // Protection a <<
492 460
493 delete[] running; << 461 if (compoundFraction > 1.0) compoundFraction = 1.0; // Protection against unphysical interpolation
494 462
>> 463 delete [] running;
>> 464
495 // get cosTh 465 // get cosTh
496 G4double incidentEnergy = anEnergy; 466 G4double incidentEnergy = anEnergy;
497 G4double incidentMass = theProjectile->Get 467 G4double incidentMass = theProjectile->GetPDGMass();
498 G4double productEnergy = fsEnergy; 468 G4double productEnergy = fsEnergy;
499 G4double productMass = result->GetMass(); 469 G4double productMass = result->GetMass();
500 auto targetZ = G4int(fCache.Get().theTarge << 470 G4int targetZ = G4int(fCache.Get()->theTargetCode/1000);
501 auto targetA = G4int(fCache.Get().theTarge << 471 G4int targetA = G4int(fCache.Get()->theTargetCode-1000*targetZ);
502 472
503 // To correspond to natural composition (- << 473 // To correspond to natural composition (-nat-) data files.
504 if (targetA == 0) targetA = G4int(fCache.G << 474 if (targetA == 0)
505 G4double targetMass = fCache.Get().theTarg << 475 targetA = G4int ( fCache.Get()->theTarget->GetMass()/amu_c2 + 0.5 );
506 auto incidentA = G4int(incidentMass / amu_ << 476 G4double targetMass = fCache.Get()->theTarget->GetMass();
507 auto incidentZ = G4int(theProjectile->GetP << 477 G4int incidentA=G4int(incidentMass/amu_c2 + 0.5 );
508 G4int residualA = targetA + incidentA - A; << 478 G4int incidentZ=G4int(theProjectile->GetPDGCharge()+ 0.5 );
509 G4int residualZ = targetZ + incidentZ - Z; << 479 G4int residualA = targetA+incidentA-A;
510 G4double residualMass = G4NucleiProperties << 480 G4int residualZ = targetZ+incidentZ-Z;
511 << 481 G4double residualMass =G4NucleiProperties::GetNuclearMass( residualA , residualZ );
512 G4ParticleHPKallbachMannSyst theKallbach( << 482 G4ParticleHPKallbachMannSyst theKallbach(compoundFraction,
513 compoundFraction, incidentEnergy, incide << 483 incidentEnergy, incidentMass,
514 residualA, residualZ, targetMass, target << 484 productEnergy, productMass,
>> 485 residualMass, residualA, residualZ,
>> 486 targetMass, targetA, targetZ,
>> 487 incidentA,incidentZ,A,Z);
515 cosTh = theKallbach.Sample(anEnergy); 488 cosTh = theKallbach.Sample(anEnergy);
516 // end (angularRep == 2) branch 489 // end (angularRep == 2) branch
517 } << 490
518 else if (angularRep > 10 && angularRep < 16) << 491 } else if (angularRep > 10 && angularRep < 16) {
519 G4double random = G4UniformRand(); 492 G4double random = G4UniformRand();
520 auto running = new G4double[nEnergies]; << 493 G4double* running = new G4double[nEnergies];
521 running[0] = 0; << 494 running[0]=0;
522 G4double weighted = 0; 495 G4double weighted = 0;
523 for (i = 1; i < nEnergies; ++i) { << 496 for (i = 1; i < nEnergies; i++) {
524 if (i != 0) running[i] = running[i - 1]; << 497 if (i != 0) running[i] = running[i-1];
525 running[i] += theInt.GetBinIntegral(theM << 498 running[i] += theInt.GetBinIntegral(theManager.GetScheme(i-1),
526 theA << 499 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
527 theA << 500 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
528 weighted += theInt.GetWeightedBinIntegra << 501 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
529 theManager.GetScheme(i - 1), theAngula << 502 theAngular[i-1].GetLabel(), theAngular[i].GetLabel(),
530 theAngular[i - 1].GetValue(0), theAngu << 503 theAngular[i-1].GetValue(0), theAngular[i].GetValue(0));
531 } 504 }
532 505
533 // Cache the mean energy in this distribut 506 // Cache the mean energy in this distribution
534 if (nEnergies == 1) << 507 if (nEnergies == 1)
535 fCache.Get().currentMeanEnergy = 0.0; << 508 fCache.Get()->currentMeanEnergy = 0.0;
536 else 509 else
537 fCache.Get().currentMeanEnergy = weighte << 510 fCache.Get()->currentMeanEnergy = weighted/running[nEnergies-1];
538 << 511
539 if (nEnergies == 1) it = 0; << 512 if (nEnergies == 1) it = 0;
540 for (i = 1; i < nEnergies; i++) { 513 for (i = 1; i < nEnergies; i++) {
541 it = i; 514 it = i;
542 if (random < running[i] / running[nEnerg << 515 if(random<running[i]/running[nEnergies-1]) break;
543 } 516 }
544 517
545 if (it < nDiscreteEnergies || it == 0) { 518 if (it < nDiscreteEnergies || it == 0) {
546 if (it == 0) { 519 if (it == 0) {
547 fsEnergy = theAngular[0].GetLabel(); << 520 fsEnergy = theAngular[0].GetLabel();
548 G4ParticleHPVector theStore; << 521 G4ParticleHPVector theStore;
549 G4int aCounter = 0; 522 G4int aCounter = 0;
550 for (G4int j = 1; j < nAngularParamete << 523 for (G4int j=1; j<nAngularParameters; j+=2) {
551 theStore.SetX(aCounter, theAngular[0 524 theStore.SetX(aCounter, theAngular[0].GetValue(j));
552 theStore.SetY(aCounter, theAngular[0 << 525 theStore.SetY(aCounter, theAngular[0].GetValue(j+1));
553 aCounter++; << 526 aCounter++;
554 } 527 }
555 G4InterpolationManager aMan; 528 G4InterpolationManager aMan;
556 aMan.Init(angularRep - 10, nAngularPar << 529 aMan.Init(angularRep-10, nAngularParameters-1);
557 theStore.SetInterpolationManager(aMan) 530 theStore.SetInterpolationManager(aMan);
558 cosTh = theStore.Sample(); 531 cosTh = theStore.Sample();
559 } << 532 } else {
560 else { <<
561 fsEnergy = theAngular[it].GetLabel(); 533 fsEnergy = theAngular[it].GetLabel();
562 G4ParticleHPVector theStore; << 534 G4ParticleHPVector theStore;
563 G4InterpolationManager aMan; 535 G4InterpolationManager aMan;
564 aMan.Init(angularRep - 10, nAngularPar << 536 aMan.Init(angularRep-10, nAngularParameters-1);
565 theStore.SetInterpolationManager(aMan) << 537 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
566 G4InterpolationScheme currentScheme = 538 G4InterpolationScheme currentScheme = theManager.GetInverseScheme(it);
567 G4int aCounter = 0; 539 G4int aCounter = 0;
568 for (G4int j = 1; j < nAngularParamete << 540 for (G4int j=1; j<nAngularParameters; j+=2) {
569 theStore.SetX(aCounter, theAngular[i 541 theStore.SetX(aCounter, theAngular[it].GetValue(j));
570 theStore.SetY(aCounter, theInt.Inter << 542 theStore.SetY(aCounter,
571 << 543 theInt.Interpolate(currentScheme, random,
572 << 544 running[it-1]/running[nEnergies-1],
573 << 545 running[it]/running[nEnergies-1],
574 << 546 theAngular[it-1].GetValue(j+1),
575 ++aCounter; << 547 theAngular[it].GetValue(j+1) ) );
>> 548 aCounter++;
576 } 549 }
577 cosTh = theStore.Sample(); 550 cosTh = theStore.Sample();
578 } 551 }
579 } << 552
580 else { << 553 } else {
581 G4double x1 = running[it - 1] / running[ << 554 G4double x1 = running[it-1]/running[nEnergies-1];
582 G4double x2 = running[it] / running[nEne << 555 G4double x2 = running[it]/running[nEnergies-1];
583 G4double y1 = theAngular[it - 1].GetLabe << 556 G4double y1 = theAngular[it-1].GetLabel();
584 G4double y2 = theAngular[it].GetLabel(); 557 G4double y2 = theAngular[it].GetLabel();
585 fsEnergy = theInt.Interpolate(theManager << 558 fsEnergy = theInt.Interpolate(theManager.GetInverseScheme(it),
>> 559 random,x1,x2,y1,y2);
586 G4ParticleHPVector theBuff1; 560 G4ParticleHPVector theBuff1;
587 G4ParticleHPVector theBuff2; 561 G4ParticleHPVector theBuff2;
588 G4InterpolationManager aMan; 562 G4InterpolationManager aMan;
589 aMan.Init(angularRep - 10, nAngularParam << 563 aMan.Init(angularRep-10, nAngularParameters-1);
590 564
591 G4int j; 565 G4int j;
592 for (i = 0, j = 1; i < nAngularParameter << 566 for (i = 0, j = 1; i < nAngularParameters; i++, j+=2) {
593 theBuff1.SetX(i, theAngular[it - 1].Ge << 567 theBuff1.SetX(i, theAngular[it-1].GetValue(j));
594 theBuff1.SetY(i, theAngular[it - 1].Ge << 568 theBuff1.SetY(i, theAngular[it-1].GetValue(j+1));
595 theBuff2.SetX(i, theAngular[it].GetVal 569 theBuff2.SetX(i, theAngular[it].GetValue(j));
596 theBuff2.SetY(i, theAngular[it].GetVal << 570 theBuff2.SetY(i, theAngular[it].GetValue(j+1));
597 } 571 }
598 572
599 G4ParticleHPVector theStore; 573 G4ParticleHPVector theStore;
600 theStore.SetInterpolationManager(aMan); << 574 theStore.SetInterpolationManager(aMan); // Store interpolates f(costh)
601 x1 = y1; 575 x1 = y1;
602 x2 = y2; 576 x2 = y2;
603 G4double x, y; 577 G4double x, y;
604 for (i = 0; i < theBuff1.GetVectorLength 578 for (i = 0; i < theBuff1.GetVectorLength(); i++) {
605 x = theBuff1.GetX(i); // costh binnin << 579 x = theBuff1.GetX(i); // costh binning identical
606 y1 = theBuff1.GetY(i); 580 y1 = theBuff1.GetY(i);
607 y2 = theBuff2.GetY(i); 581 y2 = theBuff2.GetY(i);
608 y = theInt.Interpolate(theManager.GetS << 582 y = theInt.Interpolate(theManager.GetScheme(it),
>> 583 fsEnergy, theAngular[it-1].GetLabel(),
609 theAngular[it]. 584 theAngular[it].GetLabel(), y1, y2);
610 theStore.SetX(i, x); 585 theStore.SetX(i, x);
611 theStore.SetY(i, y); 586 theStore.SetY(i, y);
612 } 587 }
613 cosTh = theStore.Sample(); 588 cosTh = theStore.Sample();
614 } 589 }
615 delete[] running; << 590 delete [] running;
616 } << 591
617 else { << 592 } else {
618 throw G4HadronicException(__FILE__, __LINE 593 throw G4HadronicException(__FILE__, __LINE__,
619 "G4ParticleHPCon << 594 "G4ParticleHPContAngularPar::Sample: Unknown angular representation");
620 } 595 }
621 //G4cout << " Efin=" << fsEnergy << G4endl; <<
622 result->SetKineticEnergy(fsEnergy); <<
623 596
624 G4double phi = twopi * G4UniformRand(); << 597 result->SetKineticEnergy(fsEnergy);
625 if(cosTh > 1.0) { cosTh = 1.0; } << 598 G4double phi = twopi*G4UniformRand();
626 else if (cosTh < -1.0) { cosTh = -1.0; } << 599 G4double theta = std::acos(cosTh);
627 G4double sinth = std::sqrt((1.0 - cosTh)*(1. << 600 G4double sinth = std::sin(theta);
628 G4double mtot = result->GetTotalMomentum(); 601 G4double mtot = result->GetTotalMomentum();
629 G4ThreeVector tempVector(mtot * sinth * std: << 602 G4ThreeVector tempVector(mtot*sinth*std::cos(phi), mtot*sinth*std::sin(phi), mtot*std::cos(theta) );
630 result->SetMomentum(tempVector); 603 result->SetMomentum(tempVector);
631 return result; 604 return result;
632 } 605 }
633 606
>> 607
>> 608 #define MERGE_NEW
>> 609
634 void G4ParticleHPContAngularPar::PrepareTableI 610 void G4ParticleHPContAngularPar::PrepareTableInterpolation()
635 { 611 {
636 // Discrete energies: store own energies in 612 // Discrete energies: store own energies in a map for faster searching
637 // 613 //
638 // The data files sometimes have identical d 614 // The data files sometimes have identical discrete energies (likely typos)
639 // which would lead to overwriting the alrea 615 // which would lead to overwriting the already existing index and hence
640 // creating a hole in the lookup table. 616 // creating a hole in the lookup table.
641 // No attempt is made here to correct for th 617 // No attempt is made here to correct for the energies - rather an epsilon
642 // is subtracted from the energy in order to 618 // is subtracted from the energy in order to uniquely identify the line
643 619
644 for (G4int ie = 0; ie < nDiscreteEnergies; i 620 for (G4int ie = 0; ie < nDiscreteEnergies; ie++) {
645 // check if energy is already present and 621 // check if energy is already present and subtract epsilon if that's the case
646 G4double myE = theAngular[ie].GetLabel(); 622 G4double myE = theAngular[ie].GetLabel();
647 while (theDiscreteEnergiesOwn.find(myE) != << 623 while (theDiscreteEnergiesOwn.find(myE) != theDiscreteEnergiesOwn.end() ) {
648 myE -= 1e-6; 624 myE -= 1e-6;
649 } 625 }
650 theDiscreteEnergiesOwn[myE] = ie; 626 theDiscreteEnergiesOwn[myE] = ie;
651 } 627 }
>> 628
>> 629 /*
>> 630 * the approach here makes no sense. It would work only for two sets that
>> 631 * have identical min and max energy. If the 2 sets differ in min, max or
>> 632 * both, the energy inserted would be normalized to its original set but
>> 633 * interpreted with the new - which is not correct.
>> 634 *
>> 635 * Disable the code for now and simply return ...
>> 636 */
>> 637
652 return; 638 return;
>> 639
>> 640 /*
>> 641 *
>> 642
>> 643 if( !angParPrev ) return;
>> 644
>> 645 //----- Discrete energies: use energies that appear in one or another
>> 646 for(ie=0; ie<nDiscreteEnergies; ie++) {
>> 647 theDiscreteEnergies.insert(theAngular[ie].GetLabel());
>> 648 }
>> 649 G4int nDiscreteEnergiesPrev = angParPrev->GetNDiscreteEnergies();
>> 650 for(ie=0; ie<nDiscreteEnergiesPrev; ie++) {
>> 651 theDiscreteEnergies.insert(angParPrev->theAngular[ie].GetLabel());
>> 652 }
>> 653
>> 654 //--- Get the values for which interpolation will be done : all energies of this and previous ContAngularPar
>> 655 for(ie=nDiscreteEnergies; ie<nEnergies; ie++) {
>> 656 G4double ener = theAngular[ie].GetLabel();
>> 657 G4double enerT = (ener-theMinEner)/(theMaxEner-theMinEner);
>> 658 theEnergiesTransformed.insert(enerT);
>> 659 }
>> 660
>> 661 G4int nEnergiesPrev = angParPrev->GetNEnergies();
>> 662 G4double minEnerPrev = angParPrev->GetMinEner();
>> 663 G4double maxEnerPrev = angParPrev->GetMaxEner();
>> 664 for(ie=nDiscreteEnergiesPrev; ie<nEnergiesPrev; ie++) {
>> 665 G4double ener = angParPrev->theAngular[ie].GetLabel();
>> 666 G4double enerT = (ener-minEnerPrev)/(maxEnerPrev-minEnerPrev);
>> 667 theEnergiesTransformed.insert(enerT);
>> 668 }
>> 669 // add the maximum energy
>> 670 //theEnergiesTransformed.insert(1.);
>> 671
>> 672 *
>> 673 */
653 } 674 }
654 675
655 void G4ParticleHPContAngularPar::BuildByInterp << 676
656 << 677 void
657 << 678 G4ParticleHPContAngularPar::BuildByInterpolation(G4double anEnergy,
658 << 679 G4InterpolationScheme aScheme,
>> 680 G4ParticleHPContAngularPar& angpar1,
>> 681 G4ParticleHPContAngularPar& angpar2)
659 { 682 {
660 G4int ie, ie1, ie2, ie1Prev, ie2Prev; 683 G4int ie, ie1, ie2, ie1Prev, ie2Prev;
661 // Only rebuild the interpolation table if t 684 // Only rebuild the interpolation table if there is a new interaction.
662 // For several subsequent samplings of final 685 // For several subsequent samplings of final state particles in the same
663 // interaction the existing table should be 686 // interaction the existing table should be used
664 if (!fCache.Get().fresh) return; << 687 if (fCache.Get()->fresh != true) return;
665 688
666 // Make copies of angpar1 and angpar2. Since 689 // Make copies of angpar1 and angpar2. Since these are given by reference
667 // it can not be excluded that one of them i 690 // it can not be excluded that one of them is "this". Hence this code uses
668 // potentially the old "this" for creating t 691 // potentially the old "this" for creating the new this - which leads to
669 // memory corruption if the old is not store 692 // memory corruption if the old is not stored as separarte object for lookup
670 const G4ParticleHPContAngularPar copyAngpar1 << 693 const G4ParticleHPContAngularPar copyAngpar1(angpar1),copyAngpar2(angpar2);
671 << 694
672 nAngularParameters = copyAngpar1.nAngularPar 695 nAngularParameters = copyAngpar1.nAngularParameters;
673 theManager = copyAngpar1.theManager; 696 theManager = copyAngpar1.theManager;
674 theEnergy = anEnergy; 697 theEnergy = anEnergy;
675 theMinEner = DBL_MAX; // min and max will b << 698 theMinEner = DBL_MAX; // min and max will be re-calculated after interpolation
676 theMaxEner = -DBL_MAX; 699 theMaxEner = -DBL_MAX;
677 700
678 // The two discrete sets must be merged. A v 701 // The two discrete sets must be merged. A vector holds the temporary data to
679 // be copied to the array in the end. Since 702 // be copied to the array in the end. Since the G4ParticleHPList class
680 // contains pointers, can't simply assign el 703 // contains pointers, can't simply assign elements of this type. Each member
681 // needs to call the explicit Set() method i 704 // needs to call the explicit Set() method instead.
682 705
683 // First, average probabilities for those li 706 // First, average probabilities for those lines that are in both sets
684 const std::map<G4double, G4int> discEnerOwn1 << 707 const std::map<G4double,G4int> discEnerOwn1 = copyAngpar1.GetDiscreteEnergiesOwn();
685 const std::map<G4double, G4int> discEnerOwn2 << 708 const std::map<G4double,G4int> discEnerOwn2 = copyAngpar2.GetDiscreteEnergiesOwn();
686 std::map<G4double, G4int>::const_iterator it << 709 std::map<G4double,G4int>::const_iterator itedeo1;
687 std::map<G4double, G4int>::const_iterator it << 710 std::map<G4double,G4int>::const_iterator itedeo2;
688 std::vector<G4ParticleHPList*> vAngular(disc << 711 std::vector<G4ParticleHPList*> vAngular(discEnerOwn1.size() );
689 G4double discEner1; 712 G4double discEner1;
690 for (itedeo1 = discEnerOwn1.cbegin(); itedeo 713 for (itedeo1 = discEnerOwn1.cbegin(); itedeo1 != discEnerOwn1.cend(); ++itedeo1) {
691 discEner1 = itedeo1->first; 714 discEner1 = itedeo1->first;
692 if (discEner1 < theMinEner) { 715 if (discEner1 < theMinEner) {
693 theMinEner = discEner1; 716 theMinEner = discEner1;
694 } 717 }
695 if (discEner1 > theMaxEner) { 718 if (discEner1 > theMaxEner) {
696 theMaxEner = discEner1; 719 theMaxEner = discEner1;
697 } 720 }
698 ie1 = itedeo1->second; 721 ie1 = itedeo1->second;
699 itedeo2 = discEnerOwn2.find(discEner1); 722 itedeo2 = discEnerOwn2.find(discEner1);
700 if (itedeo2 == discEnerOwn2.cend()) { << 723 if( itedeo2 == discEnerOwn2.cend() ) {
701 ie2 = -1; 724 ie2 = -1;
702 } << 725 } else {
703 else { <<
704 ie2 = itedeo2->second; 726 ie2 = itedeo2->second;
705 } 727 }
706 vAngular[ie1] = new G4ParticleHPList(); 728 vAngular[ie1] = new G4ParticleHPList();
707 vAngular[ie1]->SetLabel(copyAngpar1.theAng 729 vAngular[ie1]->SetLabel(copyAngpar1.theAngular[ie1].GetLabel());
708 G4double val1, val2; 730 G4double val1, val2;
709 for (G4int ip = 0; ip < nAngularParameters 731 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
710 val1 = copyAngpar1.theAngular[ie1].GetVa 732 val1 = copyAngpar1.theAngular[ie1].GetValue(ip);
711 if (ie2 != -1) { 733 if (ie2 != -1) {
712 val2 = copyAngpar2.theAngular[ie2].Get 734 val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
713 } << 735 } else {
714 else { <<
715 val2 = 0.; 736 val2 = 0.;
716 } 737 }
717 G4double value = theInt.Interpolate(aSch << 738 G4double value = theInt.Interpolate(aScheme, anEnergy,
718 copy << 739 copyAngpar1.theEnergy, copyAngpar2.theEnergy,
>> 740 val1, val2);
719 vAngular[ie1]->SetValue(ip, value); 741 vAngular[ie1]->SetValue(ip, value);
720 } 742 }
721 } // itedeo1 loop << 743 } // itedeo1 loop
722 744
723 // Add the ones in set2 but not in set1 745 // Add the ones in set2 but not in set1
724 std::vector<G4ParticleHPList*>::const_iterat 746 std::vector<G4ParticleHPList*>::const_iterator itv;
725 G4double discEner2; 747 G4double discEner2;
726 for (itedeo2 = discEnerOwn2.cbegin(); itedeo 748 for (itedeo2 = discEnerOwn2.cbegin(); itedeo2 != discEnerOwn2.cend(); ++itedeo2) {
727 discEner2 = itedeo2->first; 749 discEner2 = itedeo2->first;
728 ie2 = itedeo2->second; 750 ie2 = itedeo2->second;
729 G4bool notFound = true; 751 G4bool notFound = true;
730 itedeo1 = discEnerOwn1.find(discEner2); 752 itedeo1 = discEnerOwn1.find(discEner2);
731 if (itedeo1 != discEnerOwn1.cend()) { << 753 if (itedeo1 != discEnerOwn1.cend() ) {
732 notFound = false; 754 notFound = false;
733 } 755 }
734 if (notFound) { 756 if (notFound) {
735 // not yet in list 757 // not yet in list
736 if (discEner2 < theMinEner) { 758 if (discEner2 < theMinEner) {
737 theMinEner = discEner2; 759 theMinEner = discEner2;
738 } 760 }
739 if (discEner2 > theMaxEner) { 761 if (discEner2 > theMaxEner) {
740 theMaxEner = discEner2; 762 theMaxEner = discEner2;
741 } 763 }
742 // find position to insert 764 // find position to insert
743 G4bool isInserted = false; 765 G4bool isInserted = false;
744 ie = 0; 766 ie = 0;
745 for (itv = vAngular.cbegin(); itv != vAn << 767 for (itv = vAngular.cbegin(); itv != vAngular.cend(); ++itv,++ie) {
746 if (discEner2 > (*itv)->GetLabel()) { << 768 if (discEner2 > (*itv)->GetLabel() ) {
747 itv = vAngular.insert(itv, new G4Par 769 itv = vAngular.insert(itv, new G4ParticleHPList);
748 (*itv)->SetLabel(copyAngpar2.theAngu 770 (*itv)->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
749 isInserted = true; 771 isInserted = true;
750 break; 772 break;
751 } 773 }
752 } 774 }
753 if (!isInserted) { 775 if (!isInserted) {
754 ie = (G4int)vAngular.size(); << 776 ie=(G4int)vAngular.size();
755 vAngular.push_back(new G4ParticleHPLis 777 vAngular.push_back(new G4ParticleHPList);
756 vAngular[ie]->SetLabel(copyAngpar2.the 778 vAngular[ie]->SetLabel(copyAngpar2.theAngular[ie2].GetLabel());
757 isInserted = true; 779 isInserted = true;
758 } 780 }
759 781
760 G4double val1, val2; 782 G4double val1, val2;
761 for (G4int ip = 0; ip < nAngularParamete 783 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
762 val1 = 0; 784 val1 = 0;
763 val2 = copyAngpar2.theAngular[ie2].Get 785 val2 = copyAngpar2.theAngular[ie2].GetValue(ip);
764 G4double value = theInt.Interpolate(aS << 786 G4double value = theInt.Interpolate(aScheme, anEnergy,
765 co << 787 copyAngpar1.theEnergy,
>> 788 copyAngpar2.theEnergy,
>> 789 val1, val2);
766 vAngular[ie]->SetValue(ip, value); 790 vAngular[ie]->SetValue(ip, value);
767 } 791 }
768 } // end if(notFound) << 792 } // end if(notFound)
769 } // end loop on itedeo2 << 793 } // end loop on itedeo2
770 << 794
771 // Store new discrete list << 795 // Store new discrete list
772 nDiscreteEnergies = (G4int)vAngular.size(); 796 nDiscreteEnergies = (G4int)vAngular.size();
773 delete[] theAngular; << 797 if (theAngular != 0) delete [] theAngular;
774 theAngular = nullptr; << 798 theAngular = 0;
775 if (nDiscreteEnergies > 0) { 799 if (nDiscreteEnergies > 0) {
776 theAngular = new G4ParticleHPList[nDiscret << 800 theAngular = new G4ParticleHPList [nDiscreteEnergies];
777 } 801 }
778 theDiscreteEnergiesOwn.clear(); 802 theDiscreteEnergiesOwn.clear();
779 theDiscreteEnergies.clear(); 803 theDiscreteEnergies.clear();
780 for (ie = 0; ie < nDiscreteEnergies; ++ie) { 804 for (ie = 0; ie < nDiscreteEnergies; ++ie) {
781 theAngular[ie].SetLabel(vAngular[ie]->GetL << 805 theAngular[ie].SetLabel(vAngular[ie]->GetLabel() );
782 for (G4int ip = 0; ip < nAngularParameters 806 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
783 theAngular[ie].SetValue(ip, vAngular[ie] 807 theAngular[ie].SetValue(ip, vAngular[ie]->GetValue(ip));
784 } 808 }
785 theDiscreteEnergiesOwn[theAngular[ie].GetL 809 theDiscreteEnergiesOwn[theAngular[ie].GetLabel()] = ie;
786 theDiscreteEnergies.insert(theAngular[ie]. 810 theDiscreteEnergies.insert(theAngular[ie].GetLabel());
787 } 811 }
788 812
789 // The continuous energies need to be made f 813 // The continuous energies need to be made from scratch like the discrete
790 // ones. Therefore the re-assignemnt of theA 814 // ones. Therefore the re-assignemnt of theAngular needs to be done
791 // after the continuous energy set is also f 815 // after the continuous energy set is also finalized. Only then the
792 // total number of nEnergies is known and th 816 // total number of nEnergies is known and the array can be allocated.
793 817
794 // Get minimum and maximum energy interpolat 818 // Get minimum and maximum energy interpolating
795 // Don't use theMinEner or theMaxEner here, 819 // Don't use theMinEner or theMaxEner here, since the transformed energies
796 // need the interpolated range from the orig << 820 // need the interpolated range from the original Angpar
797 G4double interMinEner = copyAngpar1.GetMinEn << 821 G4double interMinEner = copyAngpar1.GetMinEner() + (theEnergy-copyAngpar1.GetEnergy() )
798 + (theEnergy - copyA << 822 * (copyAngpar2.GetMinEner() - copyAngpar1.GetMinEner() )
799 * (copyAngpar2.G << 823 / (copyAngpar2.GetEnergy()-copyAngpar1.GetEnergy() );
800 / (copyAngpar2.G << 824 G4double interMaxEner = copyAngpar1.GetMaxEner() + (theEnergy-copyAngpar1.GetEnergy() )
801 G4double interMaxEner = copyAngpar1.GetMaxEn << 825 * (copyAngpar2.GetMaxEner()-copyAngpar1.GetMaxEner() )
802 + (theEnergy - copyA << 826 / (copyAngpar2.GetEnergy()-copyAngpar1.GetEnergy() );
803 * (copyAngpar2.G <<
804 / (copyAngpar2.G <<
805 827
806 // Loop to energies of new set 828 // Loop to energies of new set
807 theEnergiesTransformed.clear(); 829 theEnergiesTransformed.clear();
808 830
809 G4int nEnergies1 = copyAngpar1.GetNEnergies( 831 G4int nEnergies1 = copyAngpar1.GetNEnergies();
810 G4int nDiscreteEnergies1 = copyAngpar1.GetND 832 G4int nDiscreteEnergies1 = copyAngpar1.GetNDiscreteEnergies();
811 G4double minEner1 = copyAngpar1.GetMinEner() 833 G4double minEner1 = copyAngpar1.GetMinEner();
812 G4double maxEner1 = copyAngpar1.GetMaxEner() 834 G4double maxEner1 = copyAngpar1.GetMaxEner();
813 G4int nEnergies2 = copyAngpar2.GetNEnergies( 835 G4int nEnergies2 = copyAngpar2.GetNEnergies();
814 G4int nDiscreteEnergies2 = copyAngpar2.GetND 836 G4int nDiscreteEnergies2 = copyAngpar2.GetNDiscreteEnergies();
815 G4double minEner2 = copyAngpar2.GetMinEner() 837 G4double minEner2 = copyAngpar2.GetMinEner();
816 G4double maxEner2 = copyAngpar2.GetMaxEner() 838 G4double maxEner2 = copyAngpar2.GetMaxEner();
817 839
818 // First build the list of transformed energ 840 // First build the list of transformed energies normalized
819 // to the new min max by assuming that the m 841 // to the new min max by assuming that the min-max range of
820 // each set would be scalable to the new, in 842 // each set would be scalable to the new, interpolated min
821 // max range 843 // max range
822 844
823 G4double e1(0.); 845 G4double e1(0.);
824 G4double eTNorm1(0.); 846 G4double eTNorm1(0.);
825 for (ie1 = nDiscreteEnergies1; ie1 < nEnergi 847 for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
826 e1 = copyAngpar1.theAngular[ie1].GetLabel( 848 e1 = copyAngpar1.theAngular[ie1].GetLabel();
827 eTNorm1 = (e1 - minEner1); 849 eTNorm1 = (e1 - minEner1);
828 if (maxEner1 != minEner1) eTNorm1 /= (maxE << 850 if (maxEner1 != minEner1) eTNorm1 /= (maxEner1-minEner1);
829 if (eTNorm1 >= 0 && eTNorm1 <= 1) theEnerg 851 if (eTNorm1 >= 0 && eTNorm1 <= 1) theEnergiesTransformed.insert(eTNorm1);
830 } 852 }
831 853
832 G4double e2(0.); 854 G4double e2(0.);
833 G4double eTNorm2(0.); 855 G4double eTNorm2(0.);
834 for (ie2 = nDiscreteEnergies2; ie2 < nEnergi 856 for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
835 e2 = copyAngpar2.theAngular[ie2].GetLabel( 857 e2 = copyAngpar2.theAngular[ie2].GetLabel();
836 eTNorm2 = (e2 - minEner2); 858 eTNorm2 = (e2 - minEner2);
837 if (maxEner2 != minEner2) eTNorm2 /= (maxE << 859 if (maxEner2 != minEner2) eTNorm2 /= (maxEner2-minEner2);
838 if (eTNorm2 >= 0 && eTNorm2 <= 1) theEnerg 860 if (eTNorm2 >= 0 && eTNorm2 <= 1) theEnergiesTransformed.insert(eTNorm2);
839 } 861 }
840 862
841 // Now the list of energies is complete 863 // Now the list of energies is complete
842 nEnergies = nDiscreteEnergies + (G4int)theEn << 864 nEnergies = nDiscreteEnergies+(G4int)theEnergiesTransformed.size();
843 << 865
844 // Create final array of angular parameters 866 // Create final array of angular parameters
845 const std::size_t esize = nEnergies > 0 ? nE << 867 G4ParticleHPList* theNewAngular = new G4ParticleHPList [nEnergies];
846 auto theNewAngular = new G4ParticleHPList[es <<
847 868
848 // Copy discrete energies and interpolated p 869 // Copy discrete energies and interpolated parameters to new array
849 << 870
850 if (theAngular != nullptr) { << 871 if (theAngular != 0) {
851 for (ie = 0; ie < nDiscreteEnergies; ++ie) << 872 for (ie = 0; ie < nDiscreteEnergies; ++ie) {
852 theNewAngular[ie].SetLabel(theAngular[ie 873 theNewAngular[ie].SetLabel(theAngular[ie].GetLabel());
853 for (G4int ip = 0; ip < nAngularParamete 874 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
854 theNewAngular[ie].SetValue(ip, theAngu << 875 theNewAngular[ie].SetValue(ip,theAngular[ie].GetValue(ip));
855 } 876 }
856 } 877 }
857 delete[] theAngular; << 878 delete [] theAngular;
858 } 879 }
859 theAngular = theNewAngular; 880 theAngular = theNewAngular;
860 << 881
861 // Interpolate the continuous energies for n 882 // Interpolate the continuous energies for new array
862 auto iteet = theEnergiesTransformed.begin(); << 883 std::set<G4double>::const_iterator iteet = theEnergiesTransformed.begin();
863 884
864 G4double e1Interp(0.); 885 G4double e1Interp(0.);
865 G4double e2Interp(0.); 886 G4double e2Interp(0.);
866 for (ie = nDiscreteEnergies; ie < nEnergies; << 887 for (ie = nDiscreteEnergies; ie < nEnergies; ++ie, ++iteet) {
867 G4double eT = (*iteet); 888 G4double eT = (*iteet);
868 889
869 //--- Use eT1 = eT: Get energy and paramet 890 //--- Use eT1 = eT: Get energy and parameters of copyAngpar1 for this eT
870 e1Interp = (maxEner1 - minEner1) * eT + mi 891 e1Interp = (maxEner1 - minEner1) * eT + minEner1;
871 //----- Get parameter value corresponding 892 //----- Get parameter value corresponding to this e1Interp
872 for (ie1 = nDiscreteEnergies1; ie1 < nEner 893 for (ie1 = nDiscreteEnergies1; ie1 < nEnergies1; ++ie1) {
873 if ((copyAngpar1.theAngular[ie1].GetLabe << 894 if ((copyAngpar1.theAngular[ie1].GetLabel() - e1Interp) > 1.E-10*e1Interp) break;
874 } 895 }
875 ie1Prev = ie1 - 1; 896 ie1Prev = ie1 - 1;
876 if (ie1 == 0) ++ie1Prev; << 897 if (ie1 == 0) ++ie1Prev;
877 if (ie1 == nEnergies1) { 898 if (ie1 == nEnergies1) {
878 ie1--; 899 ie1--;
879 ie1Prev = ie1; 900 ie1Prev = ie1;
880 } 901 }
881 902
882 //--- Use eT2 = eT: Get energy and paramet 903 //--- Use eT2 = eT: Get energy and parameters of copyAngpar2 for this eT
883 e2Interp = (maxEner2 - minEner2) * eT + mi << 904 e2Interp = (maxEner2-minEner2) * eT + minEner2;
884 //----- Get parameter value corresponding 905 //----- Get parameter value corresponding to this e2Interp
885 for (ie2 = nDiscreteEnergies2; ie2 < nEner 906 for (ie2 = nDiscreteEnergies2; ie2 < nEnergies2; ++ie2) {
886 if ((copyAngpar2.theAngular[ie2].GetLabe << 907 if ((copyAngpar2.theAngular[ie2].GetLabel() - e2Interp) > 1.E-10*e2Interp) break;
887 } 908 }
888 ie2Prev = ie2 - 1; 909 ie2Prev = ie2 - 1;
889 if (ie2 == 0) ++ie2Prev; << 910 if (ie2 == 0) ++ie2Prev;
890 if (ie2 == nEnergies2) { 911 if (ie2 == nEnergies2) {
891 ie2--; 912 ie2--;
892 ie2Prev = ie2; 913 ie2Prev = ie2;
893 } 914 }
894 915
895 //---- Energy corresponding to energy tran << 916 //---- Energy corresponding to energy transformed
896 G4double eN = (interMaxEner - interMinEner << 917 G4double eN = (interMaxEner-interMinEner) * eT + interMinEner;
897 << 918
898 theAngular[ie].SetLabel(eN); 919 theAngular[ie].SetLabel(eN);
899 if (eN < theMinEner) { 920 if (eN < theMinEner) {
900 theMinEner = eN; 921 theMinEner = eN;
901 } 922 }
902 if (eN > theMaxEner) { 923 if (eN > theMaxEner) {
903 theMaxEner = eN; 924 theMaxEner = eN;
904 } 925 }
905 << 926
906 G4double val1(0.); 927 G4double val1(0.);
907 G4double val2(0.); 928 G4double val2(0.);
908 G4double value(0.); 929 G4double value(0.);
909 for (G4int ip = 0; ip < nAngularParameters 930 for (G4int ip = 0; ip < nAngularParameters; ++ip) {
910 val1 = theInt.Interpolate2( << 931 val1 = theInt.Interpolate2(theManager.GetScheme(ie),
911 theManager.GetScheme(ie), e1Int << 932 e1Interp,
912 copyAngpar1.theAngular[ie1].Get << 933 copyAngpar1.theAngular[ie1Prev].GetLabel(),
913 copyAngpar1.theAngular[ie1].Get << 934 copyAngpar1.theAngular[ie1].GetLabel(),
914 * (maxEner1 - minEner1); << 935 copyAngpar1.theAngular[ie1Prev].GetValue(ip),
915 val2 = theInt.Interpolate2( << 936 copyAngpar1.theAngular[ie1].GetValue(ip)) * (maxEner1-minEner1);
916 theManager.GetScheme(ie), e2Int << 937 val2 = theInt.Interpolate2(theManager.GetScheme(ie),
917 copyAngpar2.theAngular[ie2].Get << 938 e2Interp,
918 copyAngpar2.theAngular[ie2].Get << 939 copyAngpar2.theAngular[ie2Prev].GetLabel(),
919 * (maxEner2 - minEner2); << 940 copyAngpar2.theAngular[ie2].GetLabel(),
>> 941 copyAngpar2.theAngular[ie2Prev].GetValue(ip),
>> 942 copyAngpar2.theAngular[ie2].GetValue(ip)) * (maxEner2-minEner2);
920 943
921 value = theInt.Interpolate(aScheme, anEn << 944 value = theInt.Interpolate(aScheme, anEnergy,
>> 945 copyAngpar1.theEnergy, copyAngpar2.theEnergy,
922 val1, val2); 946 val1, val2);
923 if (interMaxEner != interMinEner) { 947 if (interMaxEner != interMinEner) {
924 value /= (interMaxEner - interMinEner) << 948 value /= (interMaxEner-interMinEner);
925 } << 949 } else if (value != 0) {
926 else if (value != 0) { << 950 throw G4HadronicException(__FILE__, __LINE__,
927 throw G4HadronicException(__FILE__, __ << 951 "G4ParticleHPContAngularPar::PrepareTableInterpolation interMaxEner == interMinEner and value != 0.");
928 "G4ParticleH <<
929 "interMaxEne <<
930 } 952 }
931 theAngular[ie].SetValue(ip, value); 953 theAngular[ie].SetValue(ip, value);
932 } 954 }
933 } // end loop on nDiscreteEnergies 955 } // end loop on nDiscreteEnergies
934 956
935 for (itv = vAngular.cbegin(); itv != vAngula << 957 for (itv = vAngular.cbegin(); itv != vAngular.cend(); ++itv) delete (*itv);
936 delete (*itv); << 958
937 } 959 }
938 960
>> 961
939 void G4ParticleHPContAngularPar::Dump() const 962 void G4ParticleHPContAngularPar::Dump() const
940 { 963 {
941 G4cout << theEnergy << " " << nEnergies << " << 964 G4cout << theEnergy << " " << nEnergies << " " << nDiscreteEnergies
942 << G4endl; << 965 << " " << nAngularParameters << G4endl;
943 966
944 for (G4int ii = 0; ii < nEnergies; ++ii) << 967 for (G4int ii = 0; ii < nEnergies; ii++) theAngular[ii].Dump();
945 theAngular[ii].Dump(); <<
946 } 968 }
>> 969
947 970