| 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 // 070523 bug fix for G4FPE_DEBUG on by A. How 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 080808 bug fix in Sample() and GetXsec() by 31 // 080808 bug fix in Sample() and GetXsec() by T. Koi
32 // 32 //
33 // P. Arce, June-2014 Conversion neutron_hp to 33 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
34 // 34 //
35 #include "G4ParticleHPVector.hh" 35 #include "G4ParticleHPVector.hh"
36 <<
37 #include "G4SystemOfUnits.hh" 36 #include "G4SystemOfUnits.hh"
38 #include "G4Threading.hh" 37 #include "G4Threading.hh"
39 38
40 // if the ranges do not match, constant extrap << 39 // if the ranges do not match, constant extrapolation is used.
41 G4ParticleHPVector& operator+(G4ParticleHPVect << 40 G4ParticleHPVector & operator + (G4ParticleHPVector & left, G4ParticleHPVector & right)
42 { << 41 {
43 auto result = new G4ParticleHPVector; << 42 G4ParticleHPVector * result = new G4ParticleHPVector;
44 G4int j = 0; << 43 G4int j=0;
45 G4double x; << 44 G4double x;
46 G4double y; << 45 G4double y;
47 G4int running = 0; << 46 G4int running = 0;
48 for (G4int i = 0; i < left.GetVectorLength() << 47 for(G4int i=0; i<left.GetVectorLength(); i++)
49 while (j < right.GetVectorLength()) // Lo << 48 {
50 { << 49 while(j<right.GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
51 if (right.GetX(j) < left.GetX(i) * 1.001 <<
52 x = right.GetX(j); <<
53 y = right.GetY(j) + left.GetY(x); <<
54 result->SetData(running++, x, y); <<
55 j++; <<
56 } <<
57 // else if(std::abs((right.GetX(j)-left. <<
58 else if (left.GetX(i) + right.GetX(j) == <<
59 || std::abs((right.GetX(j) - le <<
60 { 50 {
61 x = left.GetX(i); << 51 if(right.GetX(j)<left.GetX(i)*1.001)
62 y = left.GetY(i) + right.GetY(x); << 52 {
63 result->SetData(running++, x, y); << 53 x = right.GetX(j);
64 break; << 54 y = right.GetY(j)+left.GetY(x);
65 } << 55 result->SetData(running++, x, y);
66 else { << 56 j++;
67 break; << 57 }
68 } << 58 //else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
69 } << 59 else if( left.GetX(i)+right.GetX(j) == 0
70 if (j == right.GetVectorLength()) { << 60 || std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j))) > 0.001 )
71 x = left.GetX(i); << 61 {
72 y = left.GetY(i) + right.GetY(x); << 62 x = left.GetX(i);
73 result->SetData(running++, x, y); << 63 y = left.GetY(i)+right.GetY(x);
74 } << 64 result->SetData(running++, x, y);
75 } <<
76 result->ThinOut(0.02); <<
77 return *result; <<
78 } <<
79 <<
80 G4ParticleHPVector::G4ParticleHPVector() <<
81 { <<
82 theData = new G4ParticleHPDataPoint[20]; <<
83 nPoints = 20; <<
84 nEntries = 0; <<
85 Verbose = 0; <<
86 theIntegral = nullptr; <<
87 totalIntegral = -1; <<
88 isFreed = 0; <<
89 maxValue = -DBL_MAX; <<
90 the15percentBorderCash = -DBL_MAX; <<
91 the50percentBorderCash = -DBL_MAX; <<
92 label = -DBL_MAX; <<
93 } <<
94 <<
95 G4ParticleHPVector::G4ParticleHPVector(G4int n <<
96 { <<
97 nPoints = std::max(n, 20); <<
98 theData = new G4ParticleHPDataPoint[nPoints] <<
99 nEntries = 0; <<
100 Verbose = 0; <<
101 theIntegral = nullptr; <<
102 totalIntegral = -1; <<
103 isFreed = 0; <<
104 maxValue = -DBL_MAX; <<
105 the15percentBorderCash = -DBL_MAX; <<
106 the50percentBorderCash = -DBL_MAX; <<
107 label = -DBL_MAX; <<
108 } <<
109 <<
110 G4ParticleHPVector::~G4ParticleHPVector() <<
111 { <<
112 // if(Verbose==1)G4cout <<"G4ParticleHPVe <<
113 delete[] theData; <<
114 // if(Verbose==1)G4cout <<"Vector: delete <<
115 delete[] theIntegral; <<
116 // if(Verbose==1)G4cout <<"Vector: delete <<
117 theHash.Clear(); <<
118 isFreed = 1; <<
119 } <<
120 <<
121 G4ParticleHPVector& G4ParticleHPVector::operat <<
122 { <<
123 if (&right == this) return *this; <<
124 <<
125 G4int i; <<
126 <<
127 totalIntegral = right.totalIntegral; <<
128 if (right.theIntegral != nullptr) theIntegra <<
129 for (i = 0; i < right.nEntries; i++) { <<
130 SetPoint(i, right.GetPoint(i)); // copy t <<
131 if (right.theIntegral != nullptr) theInteg <<
132 } <<
133 theManager = right.theManager; <<
134 label = right.label; <<
135 <<
136 Verbose = right.Verbose; <<
137 the15percentBorderCash = right.the15percentB <<
138 the50percentBorderCash = right.the50percentB <<
139 theHash = right.theHash; <<
140 return *this; <<
141 } <<
142 <<
143 G4double G4ParticleHPVector::GetXsec(G4double <<
144 { <<
145 if (nEntries == 0) return 0; <<
146 // if(!theHash.Prepared()) Hash(); <<
147 if (!theHash.Prepared()) { <<
148 if (G4Threading::IsWorkerThread()) { <<
149 ; <<
150 } <<
151 else { <<
152 Hash(); <<
153 } <<
154 } <<
155 G4int min = theHash.GetMinIndex(e); <<
156 G4int i; <<
157 for (i = min; i < nEntries; i++) { <<
158 // if(theData[i].GetX()>e) break; <<
159 if (theData[i].GetX() >= e) break; <<
160 } <<
161 G4int low = i - 1; <<
162 G4int high = i; <<
163 if (i == 0) { <<
164 low = 0; <<
165 high = 1; <<
166 } <<
167 else if (i == nEntries) { <<
168 low = nEntries - 2; <<
169 high = nEntries - 1; <<
170 } <<
171 G4double y; <<
172 if (e < theData[nEntries - 1].GetX()) { <<
173 // Protect against doubled-up x values <<
174 // if( (theData[high].GetX()-theData[low]. <<
175 if (theData[high].GetX() != 0 <<
176 // 080808 TKDB <<
177 //&&( theData[high].GetX()-theData[low <<
178 && (std::abs((theData[high].GetX() - t <<
179 < 0.000001)) <<
180 { <<
181 y = theData[low].GetY(); <<
182 } <<
183 else { <<
184 y = theInt.Interpolate(theManager.GetSch <<
185 theData[high].Get <<
186 } <<
187 } <<
188 else { <<
189 y = theData[nEntries - 1].GetY(); <<
190 } <<
191 return y; <<
192 } <<
193 <<
194 G4int G4ParticleHPVector::GetEnergyIndex(G4dou <<
195 { <<
196 // returns energy index of the bin in which <<
197 // returns energy index of emax for NJOY <<
198 if (nEntries == 0) return 0; <<
199 if ( ( ! theHash.Prepared() ) && ( ! G4Threa <<
200 G4int min = theHash.GetMinIndex(e); <<
201 G4int i = 0; <<
202 for (i=min ; i < nEntries; i++) if (theData[ <<
203 return i; <<
204 } <<
205 <<
206 void G4ParticleHPVector::Dump() <<
207 { <<
208 G4cout << nEntries << G4endl; <<
209 for (G4int i = 0; i < nEntries; i++) { <<
210 G4cout << theData[i].GetX() << " "; <<
211 G4cout << theData[i].GetY() << " "; <<
212 // if (i!=1&&i==5*(i/5)) G4cout << G4 <<
213 G4cout << G4endl; <<
214 } <<
215 G4cout << G4endl; <<
216 } <<
217 <<
218 void G4ParticleHPVector::Check(G4int i) <<
219 { <<
220 if (i > nEntries) <<
221 throw G4HadronicException(__FILE__, __LINE <<
222 "Skipped some in <<
223 if (i == nPoints) { <<
224 nPoints = static_cast<G4int>(1.2 * nPoints <<
225 auto buff = new G4ParticleHPDataPoint[nPoi <<
226 for (G4int j = 0; j < nEntries; j++) <<
227 buff[j] = theData[j]; <<
228 delete[] theData; <<
229 theData = buff; <<
230 } <<
231 if (i == nEntries) nEntries = i + 1; <<
232 } <<
233 <<
234 void G4ParticleHPVector::Merge(G4Interpolation <<
235 G4ParticleHPVec <<
236 { <<
237 // interpolate between labels according to a <<
238 // continue in unknown areas by substraction <<
239 <<
240 CleanUp(); <<
241 G4int s_tmp = 0, n = 0, m_tmp = 0; <<
242 G4ParticleHPVector* tmp; <<
243 G4int a = s_tmp, p = n, t; <<
244 while (a < active->GetVectorLength()) // Lo <<
245 { <<
246 if (active->GetEnergy(a) <= passive->GetEn <<
247 G4double xa = active->GetEnergy(a); <<
248 G4double yy = theInt.Interpolate(aScheme <<
249 active- <<
250 SetData(m_tmp, xa, yy); <<
251 theManager.AppendScheme(m_tmp, active->G <<
252 m_tmp++; <<
253 a++; <<
254 G4double xp = passive->GetEnergy(p); <<
255 // if( std::abs(std::abs(xp-xa)/xa)<0.00 <<
256 if (xa != 0 && std::abs(std::abs(xp - xa <<
257 { <<
258 p++; <<
259 tmp = active; <<
260 t = a; <<
261 active = passive; <<
262 a = p; <<
263 passive = tmp; <<
264 p = t; <<
265 } <<
266 } <<
267 else { <<
268 tmp = active; <<
269 t = a; <<
270 active = passive; <<
271 a = p; <<
272 passive = tmp; <<
273 p = t; <<
274 } <<
275 } <<
276 <<
277 G4double deltaX = passive->GetXsec(GetEnergy <<
278 while (p != passive->GetVectorLength() <<
279 && passive->GetEnergy(p) <= aValue) <<
280 { <<
281 G4double anX; <<
282 anX = passive->GetXsec(p) - deltaX; <<
283 if (anX > 0) { <<
284 // if(std::abs(GetEnergy(m-1)-passive->G <<
285 if (passive->GetEnergy(p) == 0 <<
286 || std::abs(GetEnergy(m_tmp - 1) - p <<
287 > 0.0000001) <<
288 { <<
289 SetData(m_tmp, passive->GetEnergy(p), <<
290 theManager.AppendScheme(m_tmp++, passi <<
291 } <<
292 } <<
293 p++; <<
294 } <<
295 // Rebuild the Hash; <<
296 if (theHash.Prepared()) { <<
297 ReHash(); <<
298 } <<
299 } <<
300 <<
301 void G4ParticleHPVector::ThinOut(G4double prec <<
302 { <<
303 // anything in there? <<
304 if (GetVectorLength() == 0) return; <<
305 // make the new vector <<
306 auto aBuff = new G4ParticleHPDataPoint[nPoin <<
307 G4double x, x1, x2, y, y1, y2; <<
308 G4int count = 0, current = 2, start = 1; <<
309 <<
310 // First element always goes and is never te <<
311 aBuff[0] = theData[0]; <<
312 <<
313 // Find the rest <<
314 while (current < GetVectorLength()) // Loop <<
315 { <<
316 x1 = aBuff[count].GetX(); <<
317 y1 = aBuff[count].GetY(); <<
318 x2 = theData[current].GetX(); <<
319 y2 = theData[current].GetY(); <<
320 <<
321 if (x1 - x2 == 0) { <<
322 // Following block added for avoiding di <<
323 for (G4int j = start; j < current; j++) <<
324 y = (y2 + y1) / 2.; <<
325 if (std::abs(y - theData[j].GetY()) > <<
326 aBuff[++count] = theData[current - 1 <<
327 start = current; // the next candid <<
328 break; 65 break;
329 } 66 }
330 } <<
331 } <<
332 else { <<
333 for (G4int j = start; j < current; j++) <<
334 x = theData[j].GetX(); <<
335 if (x1 - x2 == 0) <<
336 y = (y2 + y1) / 2.; <<
337 else 67 else
338 y = theInt.Lin(x, x1, x2, y1, y2); << 68 {
339 if (std::abs(y - theData[j].GetY()) > <<
340 aBuff[++count] = theData[current - 1 <<
341 start = current; // the next candid <<
342 break; 69 break;
343 } 70 }
344 } 71 }
>> 72 if(j==right.GetVectorLength())
>> 73 {
>> 74 x = left.GetX(i);
>> 75 y = left.GetY(i)+right.GetY(x);
>> 76 result->SetData(running++, x, y);
>> 77 }
345 } 78 }
346 current++; << 79 result->ThinOut(0.02);
>> 80 return *result;
347 } 81 }
348 // The last one also always goes, and is nev << 82
349 aBuff[++count] = theData[GetVectorLength() - << 83 G4ParticleHPVector::G4ParticleHPVector()
350 delete[] theData; << 84 {
351 theData = aBuff; << 85 theData = new G4ParticleHPDataPoint[20];
352 nEntries = count + 1; << 86 nPoints=20;
353 << 87 nEntries=0;
354 // Rebuild the Hash; << 88 Verbose=0;
355 if (theHash.Prepared()) { << 89 theIntegral=0;
356 ReHash(); << 90 totalIntegral=-1;
357 } << 91 isFreed = 0;
358 } << 92 maxValue = -DBL_MAX;
359 << 93 the15percentBorderCash = -DBL_MAX;
360 G4bool G4ParticleHPVector::IsBlocked(G4double << 94 the50percentBorderCash = -DBL_MAX;
361 { << 95 label = -DBL_MAX;
362 G4bool result = false; << 96 }
363 std::vector<G4double>::iterator i; << 97
364 for (i = theBlocked.begin(); i != theBlocked << 98 G4ParticleHPVector::G4ParticleHPVector(G4int n)
365 G4double aBlock = *i; << 99 {
366 if (std::abs(aX - aBlock) < 0.1 * MeV) { << 100 nPoints=std::max(n, 20);
367 result = true; << 101 theData = new G4ParticleHPDataPoint[nPoints];
368 theBlocked.erase(i); << 102 nEntries=0;
369 break; << 103 Verbose=0;
370 } << 104 theIntegral=0;
371 } << 105 totalIntegral=-1;
372 return result; << 106 isFreed = 0;
373 } << 107 maxValue = -DBL_MAX;
374 << 108 the15percentBorderCash = -DBL_MAX;
375 G4double G4ParticleHPVector::Sample() // Samp << 109 the50percentBorderCash = -DBL_MAX;
376 { << 110 label = -DBL_MAX;
377 G4double result = 0.; << 111 }
378 G4int j; << 112
379 for (j = 0; j < GetVectorLength(); j++) { << 113 G4ParticleHPVector::~G4ParticleHPVector()
380 if (GetY(j) < 0) SetY(j, 0); << 114 {
381 } << 115 // if(Verbose==1)G4cout <<"G4ParticleHPVector::~G4ParticleHPVector"<<G4endl;
382 << 116 delete [] theData;
383 if (!theBuffered.empty() && G4UniformRand() << 117 // if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
384 result = theBuffered[0]; << 118 delete [] theIntegral;
385 theBuffered.erase(theBuffered.begin()); << 119 // if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
386 if (result < GetX(GetVectorLength() - 1)) << 120 theHash.Clear();
387 } << 121 isFreed = 1;
388 if (GetVectorLength() == 1) { << 122 }
389 result = theData[0].GetX(); << 123
390 } << 124 G4ParticleHPVector & G4ParticleHPVector::
391 else { << 125 operator = (const G4ParticleHPVector & right)
392 if (theIntegral == nullptr) { << 126 {
393 IntegrateAndNormalise(); << 127 if(&right == this) return *this;
394 } << 128
395 G4int icounter = 0; << 129 G4int i;
396 G4int icounter_max = 1024; << 130
397 do { << 131 totalIntegral = right.totalIntegral;
398 icounter++; << 132 if(right.theIntegral!=0) theIntegral = new G4double[right.nEntries];
399 if (icounter > icounter_max) { << 133 for(i=0; i<right.nEntries; i++)
400 G4cout << "Loop-counter exceeded the t << 134 {
401 << __FILE__ << "." << G4endl; << 135 SetPoint(i, right.GetPoint(i)); // copy theData
402 break; << 136 if(right.theIntegral!=0) theIntegral[i] = right.theIntegral[i];
403 } << 137 }
404 // 080808 << 138 theManager = right.theManager;
405 /* << 139 label = right.label;
406 G4double rand; << 140
407 G4double value, test, baseline; << 141 Verbose = right.Verbose;
408 baseline = theData[GetVectorLeng << 142 the15percentBorderCash = right.the15percentBorderCash;
409 do << 143 the50percentBorderCash = right.the50percentBorderCash;
410 { << 144 theHash = right.theHash;
411 value = baseline*G4UniformRand << 145 return *this;
412 value += theData[0].GetX(); << 146 }
413 test = GetY(value)/maxValue; << 147
414 rand = G4UniformRand(); << 148
415 } << 149 G4double G4ParticleHPVector::GetXsec(G4double e)
416 //while(test<rand); << 150 {
417 while( test < rand && test > 0 ) << 151 if(nEntries == 0) return 0;
418 result = value; << 152 //if(!theHash.Prepared()) Hash();
419 */ << 153 if ( !theHash.Prepared() ) {
420 G4double rand; << 154 if ( G4Threading::IsWorkerThread() ) {
421 G4double value = 0., test; << 155 ;
422 G4int jcounter = 0; << 156 } else {
423 G4int jcounter_max = 1024; << 157 Hash();
424 do { << 158 }
425 jcounter++; << 159 }
426 if (jcounter > jcounter_max) { << 160 G4int min = theHash.GetMinIndex(e);
427 G4cout << "Loop-counter exceeded the << 161 G4int i;
428 << __FILE__ << "." << G4endl; << 162 for(i=min ; i<nEntries; i++)
429 break; << 163 {
430 } << 164 //if(theData[i].GetX()>e) break;
431 rand = G4UniformRand(); << 165 if(theData[i].GetX() >= e) break;
432 G4int ibin = -1; << 166 }
433 for (G4int i = 0; i < GetVectorLength( << 167 G4int low = i-1;
434 if (rand < theIntegral[i]) { << 168 G4int high = i;
435 ibin = i; << 169 if(i==0)
436 break; << 170 {
437 } << 171 low = 0;
>> 172 high = 1;
>> 173 }
>> 174 else if(i==nEntries)
>> 175 {
>> 176 low = nEntries-2;
>> 177 high = nEntries-1;
>> 178 }
>> 179 G4double y;
>> 180 if(e<theData[nEntries-1].GetX())
>> 181 {
>> 182 // Protect against doubled-up x values
>> 183 //if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
>> 184 if ( theData[high].GetX() !=0
>> 185 //080808 TKDB
>> 186 //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
>> 187 &&( std::abs( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() ) < 0.000001 ) )
>> 188 {
>> 189 y = theData[low].GetY();
>> 190 }
>> 191 else
>> 192 {
>> 193 y = theInt.Interpolate(theManager.GetScheme(high), e,
>> 194 theData[low].GetX(), theData[high].GetX(),
>> 195 theData[low].GetY(), theData[high].GetY());
>> 196 }
>> 197 }
>> 198 else
>> 199 {
>> 200 y=theData[nEntries-1].GetY();
>> 201 }
>> 202 return y;
>> 203 }
>> 204
>> 205 void G4ParticleHPVector::Dump()
>> 206 {
>> 207 G4cout << nEntries<<G4endl;
>> 208 for(G4int i=0; i<nEntries; i++)
>> 209 {
>> 210 G4cout << theData[i].GetX()<<" ";
>> 211 G4cout << theData[i].GetY()<<" ";
>> 212 // if (i!=1&&i==5*(i/5)) G4cout << G4endl;
>> 213 G4cout << G4endl;
>> 214 }
>> 215 G4cout << G4endl;
>> 216 }
>> 217
>> 218 void G4ParticleHPVector::Check(G4int i)
>> 219 {
>> 220 if(i>nEntries) throw G4HadronicException(__FILE__, __LINE__, "Skipped some index numbers in G4ParticleHPVector");
>> 221 if(i==nPoints)
>> 222 {
>> 223 nPoints = static_cast<G4int>(1.2*nPoints);
>> 224 G4ParticleHPDataPoint * buff = new G4ParticleHPDataPoint[nPoints];
>> 225 for (G4int j=0; j<nEntries; j++) buff[j] = theData[j];
>> 226 delete [] theData;
>> 227 theData = buff;
>> 228 }
>> 229 if(i==nEntries) nEntries=i+1;
>> 230 }
>> 231
>> 232 void G4ParticleHPVector::
>> 233 Merge(G4InterpolationScheme aScheme, G4double aValue,
>> 234 G4ParticleHPVector * active, G4ParticleHPVector * passive)
>> 235 {
>> 236 // interpolate between labels according to aScheme, cut at aValue,
>> 237 // continue in unknown areas by substraction of the last difference.
>> 238
>> 239 CleanUp();
>> 240 G4int s_tmp = 0, n=0, m_tmp=0;
>> 241 G4ParticleHPVector * tmp;
>> 242 G4int a = s_tmp, p = n, t;
>> 243 while ( a<active->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi
>> 244 {
>> 245 if(active->GetEnergy(a) <= passive->GetEnergy(p))
>> 246 {
>> 247 G4double xa = active->GetEnergy(a);
>> 248 G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
>> 249 active->GetXsec(a), passive->GetXsec(xa));
>> 250 SetData(m_tmp, xa, yy);
>> 251 theManager.AppendScheme(m_tmp, active->GetScheme(a));
>> 252 m_tmp++;
>> 253 a++;
>> 254 G4double xp = passive->GetEnergy(p);
>> 255 //if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() )
>> 256 if ( xa != 0
>> 257 && std::abs(std::abs(xp-xa)/xa) < 0.0000001
>> 258 && a < active->GetVectorLength() )
>> 259 {
>> 260 p++;
>> 261 tmp = active; t=a;
>> 262 active = passive; a=p;
>> 263 passive = tmp; p=t;
438 } 264 }
439 if (ibin < 0) G4cout << "TKDB 080807 " << 265 } else {
440 // result << 266 tmp = active; t=a;
441 rand = G4UniformRand(); << 267 active = passive; a=p;
442 G4double x1, x2; << 268 passive = tmp; p=t;
443 if (ibin == 0) { << 269 }
444 x1 = theData[ibin].GetX(); << 270 }
445 value = x1; << 271
446 break; << 272 G4double deltaX = passive->GetXsec(GetEnergy(m_tmp-1)) - GetXsec(m_tmp-1);
>> 273 while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue) // Loop checking, 11.05.2015, T. Koi
>> 274 {
>> 275 G4double anX;
>> 276 anX = passive->GetXsec(p)-deltaX;
>> 277 if(anX>0)
>> 278 {
>> 279 //if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
>> 280 if ( passive->GetEnergy(p) == 0
>> 281 || std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p) > 0.0000001 )
>> 282 {
>> 283 SetData(m_tmp, passive->GetEnergy(p), anX);
>> 284 theManager.AppendScheme(m_tmp++, passive->GetScheme(p));
447 } 285 }
448 x1 = theData[ibin - 1].GetX(); << 286 }
>> 287 p++;
>> 288 }
>> 289 // Rebuild the Hash;
>> 290 if(theHash.Prepared())
>> 291 {
>> 292 ReHash();
>> 293 }
>> 294 }
>> 295
>> 296 void G4ParticleHPVector::ThinOut(G4double precision)
>> 297 {
>> 298 // anything in there?
>> 299 if(GetVectorLength()==0) return;
>> 300 // make the new vector
>> 301 G4ParticleHPDataPoint * aBuff = new G4ParticleHPDataPoint[nPoints];
>> 302 G4double x, x1, x2, y, y1, y2;
>> 303 G4int count = 0, current = 2, start = 1;
>> 304
>> 305 // First element always goes and is never tested.
>> 306 aBuff[0] = theData[0];
>> 307
>> 308 // Find the rest
>> 309 while(current < GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
>> 310 {
>> 311 x1=aBuff[count].GetX();
>> 312 y1=aBuff[count].GetY();
>> 313 x2=theData[current].GetX();
>> 314 y2=theData[current].GetY();
>> 315
>> 316 if ( x1-x2 == 0 ) {
>> 317 //Following block added for avoiding div 0 error on Release + G4FPE_DEBUG
>> 318 for ( G4int j=start; j<current; j++ ) {
>> 319 y = (y2+y1)/2.;
>> 320 if ( std::abs( y-theData[j].GetY() ) > precision*y ) {
>> 321 aBuff[++count] = theData[current-1]; // for this one, everything was fine
>> 322 start = current; // the next candidate
>> 323 break;
>> 324 }
>> 325 }
>> 326 } else {
>> 327 for(G4int j=start; j<current; j++)
>> 328 {
>> 329 x = theData[j].GetX();
>> 330 if(x1-x2 == 0) y = (y2+y1)/2.;
>> 331 else y = theInt.Lin(x, x1, x2, y1, y2);
>> 332 if (std::abs(y-theData[j].GetY())>precision*y)
>> 333 {
>> 334 aBuff[++count] = theData[current-1]; // for this one, everything was fine
>> 335 start = current; // the next candidate
>> 336 break;
>> 337 }
>> 338 }
>> 339 }
>> 340 current++ ;
>> 341 }
>> 342 // The last one also always goes, and is never tested.
>> 343 aBuff[++count] = theData[GetVectorLength()-1];
>> 344 delete [] theData;
>> 345 theData = aBuff;
>> 346 nEntries = count+1;
>> 347
>> 348 // Rebuild the Hash;
>> 349 if(theHash.Prepared())
>> 350 {
>> 351 ReHash();
>> 352 }
>> 353 }
449 354
450 x2 = theData[ibin].GetX(); << 355 G4bool G4ParticleHPVector::IsBlocked(G4double aX)
451 value = rand * (x2 - x1) + x1; << 356 {
452 //************************************ << 357 G4bool result = false;
453 /* << 358 std::vector<G4double>::iterator i;
454 test = GetY ( value ) / std::max << 359 for(i=theBlocked.begin(); i!=theBlocked.end(); i++)
455 */ << 360 {
456 //************************************ << 361 G4double aBlock = *i;
457 // EMendoza - Always linear interpolat << 362 if(std::abs(aX-aBlock) < 0.1*MeV)
458 G4double y1 = theData[ibin - 1].GetY() << 363 {
459 G4double y2 = theData[ibin].GetY(); << 364 result = true;
460 G4double mval = (y2 - y1) / (x2 - x1); << 365 theBlocked.erase(i);
461 G4double bval = y1 - mval * x1; << 366 break;
462 test = (mval * value + bval) / std::ma << 367 }
463 //************************************ <<
464 } while (G4UniformRand() > test); // Lo <<
465 result = value; <<
466 // 080807 <<
467 } while (IsBlocked(result)); // Loop chec <<
468 } <<
469 return result; <<
470 } <<
471 <<
472 G4double G4ParticleHPVector::Get15percentBorde <<
473 { <<
474 if (the15percentBorderCash > -DBL_MAX / 2.) <<
475 G4double result; <<
476 if (GetVectorLength() == 1) { <<
477 result = theData[0].GetX(); <<
478 the15percentBorderCash = result; <<
479 } <<
480 else { <<
481 if (theIntegral == nullptr) { <<
482 IntegrateAndNormalise(); <<
483 } 368 }
484 G4int i; << 369 return result;
485 result = theData[GetVectorLength() - 1].Ge << 370 }
486 for (i = 0; i < GetVectorLength(); i++) { << 371
487 if (theIntegral[i] / theIntegral[GetVect << 372 G4double G4ParticleHPVector::Sample() // Samples X according to distribution Y
488 result = theData[std::min(i + 1, GetVe << 373 {
489 the15percentBorderCash = result; << 374 G4double result;
490 break; << 375 G4int j;
491 } << 376 for(j=0; j<GetVectorLength(); j++)
492 } << 377 {
493 the15percentBorderCash = result; << 378 if(GetY(j)<0) SetY(j, 0);
494 } <<
495 return result; <<
496 } <<
497 <<
498 G4double G4ParticleHPVector::Get50percentBorde <<
499 { <<
500 if (the50percentBorderCash > -DBL_MAX / 2.) <<
501 G4double result; <<
502 if (GetVectorLength() == 1) { <<
503 result = theData[0].GetX(); <<
504 the50percentBorderCash = result; <<
505 } <<
506 else { <<
507 if (theIntegral == nullptr) { <<
508 IntegrateAndNormalise(); <<
509 } 379 }
510 G4int i; << 380
511 G4double x = 0.5; << 381 if(theBuffered.size() !=0 && G4UniformRand()<0.5)
512 result = theData[GetVectorLength() - 1].Ge << 382 {
513 for (i = 0; i < GetVectorLength(); i++) { << 383 result = theBuffered[0];
514 if (theIntegral[i] / theIntegral[GetVect << 384 theBuffered.erase(theBuffered.begin());
515 G4int it; << 385 if(result < GetX(GetVectorLength()-1) ) return result;
516 it = i; << 386 }
517 if (it == GetVectorLength() - 1) { << 387 if(GetVectorLength()==1)
518 result = theData[GetVectorLength() - << 388 {
>> 389 result = theData[0].GetX();
>> 390 }
>> 391 else
>> 392 {
>> 393 if(theIntegral==0) { IntegrateAndNormalise(); }
>> 394 G4int icounter=0;
>> 395 G4int icounter_max=1024;
>> 396 do
>> 397 {
>> 398 icounter++;
>> 399 if ( icounter > icounter_max ) {
>> 400 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
>> 401 break;
519 } 402 }
520 else { << 403 //080808
521 G4double x1, x2, y1, y2; << 404 /*
522 x1 = theIntegral[i - 1] / theIntegra << 405 G4double rand;
523 x2 = theIntegral[i] / theIntegral[Ge << 406 G4double value, test, baseline;
524 y1 = theData[i - 1].GetX(); << 407 baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
525 y2 = theData[i].GetX(); << 408 do
526 result = theLin.Lin(x, x1, x2, y1, y << 409 {
>> 410 value = baseline*G4UniformRand();
>> 411 value += theData[0].GetX();
>> 412 test = GetY(value)/maxValue;
>> 413 rand = G4UniformRand();
>> 414 }
>> 415 //while(test<rand);
>> 416 while( test < rand && test > 0 );
>> 417 result = value;
>> 418 */
>> 419 G4double rand;
>> 420 G4double value, test;
>> 421 G4int jcounter=0;
>> 422 G4int jcounter_max=1024;
>> 423 do
>> 424 {
>> 425 jcounter++;
>> 426 if ( jcounter > jcounter_max ) {
>> 427 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
>> 428 break;
>> 429 }
>> 430 rand = G4UniformRand();
>> 431 G4int ibin = -1;
>> 432 for ( G4int i = 0 ; i < GetVectorLength() ; i++ )
>> 433 {
>> 434 if ( rand < theIntegral[i] )
>> 435 {
>> 436 ibin = i;
>> 437 break;
>> 438 }
>> 439 }
>> 440 if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl;
>> 441 // result
>> 442 rand = G4UniformRand();
>> 443 G4double x1, x2;
>> 444 if ( ibin == 0 )
>> 445 {
>> 446 x1 = theData[ ibin ].GetX();
>> 447 value = x1;
>> 448 break;
>> 449 }
>> 450 else
>> 451 {
>> 452 x1 = theData[ ibin-1 ].GetX();
>> 453 }
>> 454
>> 455 x2 = theData[ ibin ].GetX();
>> 456 value = rand * ( x2 - x1 ) + x1;
>> 457 //***********************************************************************
>> 458 /*
>> 459 test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
>> 460 */
>> 461 //***********************************************************************
>> 462 //EMendoza - Always linear interpolation:
>> 463 G4double y1=theData[ ibin-1 ].GetY();
>> 464 G4double y2=theData[ ibin ].GetY();
>> 465 G4double mval=(y2-y1)/(x2-x1);
>> 466 G4double bval=y1-mval*x1;
>> 467 test =(mval*value+bval)/std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
>> 468 //***********************************************************************
527 } 469 }
528 the50percentBorderCash = result; << 470 while ( G4UniformRand() > test ); // Loop checking, 11.05.2015, T. Koi
529 break; << 471 result = value;
>> 472 //080807
530 } 473 }
>> 474 while(IsBlocked(result)); // Loop checking, 11.05.2015, T. Koi
531 } 475 }
532 the50percentBorderCash = result; << 476 return result;
533 } 477 }
534 return result; <<
535 } <<
536 478
537 G4double G4ParticleHPVector::GetMaxY(G4double << 479 G4double G4ParticleHPVector::Get15percentBorder()
538 { << 480 {
539 G4double xsmax = 0.; << 481 if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
540 if (emin > emax || nEntries == 0) return xsm << 482 G4double result;
541 if (emin >= theData[nEntries - 1].GetX()) { << 483 if(GetVectorLength()==1)
542 xsmax = theData[nEntries - 1].GetY(); << 484 {
543 return xsmax; << 485 result = theData[0].GetX();
544 } << 486 the15percentBorderCash = result;
545 if (emax <= theData[0].GetX()) { << 487 }
546 xsmax = theData[0].GetY(); << 488 else
547 return xsmax; << 489 {
548 } << 490 if(theIntegral==0) { IntegrateAndNormalise(); }
549 if (!theHash.Prepared() && !G4Threading::IsW << 491 G4int i;
550 // Find the lowest index, low, where x(energ << 492 result = theData[GetVectorLength()-1].GetX();
551 G4int i = theHash.GetMinIndex(emin); << 493 for(i=0;i<GetVectorLength();i++)
552 for (; i < nEntries; ++i) { << 494 {
553 if (theData[i].GetX() >= emin) break; << 495 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15)
554 } << 496 {
555 G4int low = i; << 497 result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
556 // Find the lowest index, high, where x(ener << 498 the15percentBorderCash = result;
557 i = theHash.GetMinIndex(emax); << 499 break;
558 for (; i < nEntries; ++i) { << 500 }
559 if (theData[i].GetX() >= emax) break; << 501 }
560 } << 502 the15percentBorderCash = result;
561 G4int high = i; << 503 }
562 xsmax = GetXsec(emin); // Set xsmax as low << 504 return result;
563 // Find the highest cross-section <<
564 for (i = low; i < high; ++i) { <<
565 if (xsmax < theData[i].GetY()) xsmax = the <<
566 } <<
567 // Check if it is smaller than the high bord <<
568 G4double highborder = GetXsec(emax); <<
569 if (xsmax < highborder) xsmax = highborder; <<
570 <<
571 if (xsmax == 0.) { <<
572 throw G4HadronicException(__FILE__, __LINE <<
573 "G4ParticleHPVec <<
574 "G4Nucleus::GetB <<
575 } 505 }
576 506
577 return xsmax; << 507 G4double G4ParticleHPVector::Get50percentBorder()
578 } << 508 {
>> 509 if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
>> 510 G4double result;
>> 511 if(GetVectorLength()==1)
>> 512 {
>> 513 result = theData[0].GetX();
>> 514 the50percentBorderCash = result;
>> 515 }
>> 516 else
>> 517 {
>> 518 if(theIntegral==0) { IntegrateAndNormalise(); }
>> 519 G4int i;
>> 520 G4double x = 0.5;
>> 521 result = theData[GetVectorLength()-1].GetX();
>> 522 for(i=0;i<GetVectorLength();i++)
>> 523 {
>> 524 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x)
>> 525 {
>> 526 G4int it;
>> 527 it = i;
>> 528 if(it == GetVectorLength()-1)
>> 529 {
>> 530 result = theData[GetVectorLength()-1].GetX();
>> 531 }
>> 532 else
>> 533 {
>> 534 G4double x1, x2, y1, y2;
>> 535 x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
>> 536 x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
>> 537 y1 = theData[i-1].GetX();
>> 538 y2 = theData[i].GetX();
>> 539 result = theLin.Lin(x, x1, x2, y1, y2);
>> 540 }
>> 541 the50percentBorderCash = result;
>> 542 break;
>> 543 }
>> 544 }
>> 545 the50percentBorderCash = result;
>> 546 }
>> 547 return result;
>> 548 }
579 549