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