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