source: trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPVector.cc@ 1350

Last change on this file since 1350 was 962, checked in by garnier, 17 years ago

update processes
File size: 14.1 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// neutron_hp -- source file
27// J.P. Wellisch, Nov-1996
28// A prototype of the low energy neutron transport model.
29//
30// 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31// 080808 bug fix in Sample() and GetXsec() by T. Koi
32//
33#include "G4NeutronHPVector.hh"
34
35 // if the ranges do not match, constant extrapolation is used.
36 G4NeutronHPVector & operator + (G4NeutronHPVector & left, G4NeutronHPVector & right)
37 {
38 G4NeutronHPVector * result = new G4NeutronHPVector;
39 G4int j=0;
40 G4double x;
41 G4double y;
42 G4int running = 0;
43 for(G4int i=0; i<left.GetVectorLength(); i++)
44 {
45 while(j<right.GetVectorLength())
46 {
47 if(right.GetX(j)<left.GetX(i)*1.001)
48 {
49 x = right.GetX(j);
50 y = right.GetY(j)+left.GetY(x);
51 result->SetData(running++, x, y);
52 j++;
53 }
54 //else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
55 else if( left.GetX(i)+right.GetX(j) == 0
56 || std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j))) > 0.001 )
57 {
58 x = left.GetX(i);
59 y = left.GetY(i)+right.GetY(x);
60 result->SetData(running++, x, y);
61 break;
62 }
63 else
64 {
65 break;
66 }
67 }
68 if(j==right.GetVectorLength())
69 {
70 x = left.GetX(i);
71 y = left.GetY(i)+right.GetY(x);
72 result->SetData(running++, x, y);
73 }
74 }
75 result->ThinOut(0.02);
76 return *result;
77 }
78
79 G4NeutronHPVector::G4NeutronHPVector()
80 {
81 theData = new G4NeutronHPDataPoint[20];
82 nPoints=20;
83 nEntries=0;
84 Verbose=0;
85 theIntegral=0;
86 totalIntegral=-1;
87 isFreed = 0;
88 maxValue = -DBL_MAX;
89 the15percentBorderCash = -DBL_MAX;
90 the50percentBorderCash = -DBL_MAX;
91 label = -DBL_MAX;
92
93 }
94
95 G4NeutronHPVector::G4NeutronHPVector(G4int n)
96 {
97 nPoints=std::max(n, 20);
98 theData = new G4NeutronHPDataPoint[nPoints];
99 nEntries=0;
100 Verbose=0;
101 theIntegral=0;
102 totalIntegral=-1;
103 isFreed = 0;
104 maxValue = -DBL_MAX;
105 the15percentBorderCash = -DBL_MAX;
106 the50percentBorderCash = -DBL_MAX;
107 }
108
109 G4NeutronHPVector::~G4NeutronHPVector()
110 {
111// if(Verbose==1)G4cout <<"G4NeutronHPVector::~G4NeutronHPVector"<<G4endl;
112 delete [] theData;
113// if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
114 delete [] theIntegral;
115// if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
116 isFreed = 1;
117 }
118
119 G4NeutronHPVector & G4NeutronHPVector::
120 operator = (const G4NeutronHPVector & right)
121 {
122 if(&right == this) return *this;
123
124 G4int i;
125
126 totalIntegral = right.totalIntegral;
127 if(right.theIntegral!=0) theIntegral = new G4double[right.nEntries];
128 for(i=0; i<right.nEntries; i++)
129 {
130 SetPoint(i, right.GetPoint(i)); // copy theData
131 if(right.theIntegral!=0) theIntegral[i] = right.theIntegral[i];
132 }
133 theManager = right.theManager;
134 label = right.label;
135
136 Verbose = right.Verbose;
137 the15percentBorderCash = right.the15percentBorderCash;
138 the50percentBorderCash = right.the50percentBorderCash;
139 theHash = right.theHash;
140 return *this;
141 }
142
143
144 G4double G4NeutronHPVector::GetXsec(G4double e)
145 {
146 if(nEntries == 0) return 0;
147 if(!theHash.Prepared()) Hash();
148 G4int min = theHash.GetMinIndex(e);
149 G4int i;
150 for(i=min ; i<nEntries; i++)
151 {
152 //if(theData[i].GetX()>e) break;
153 if(theData[i].GetX() >= e) break;
154 }
155 G4int low = i-1;
156 G4int high = i;
157 if(i==0)
158 {
159 low = 0;
160 high = 1;
161 }
162 else if(i==nEntries)
163 {
164 low = nEntries-2;
165 high = nEntries-1;
166 }
167 G4double y;
168 if(e<theData[nEntries-1].GetX())
169 {
170 // Protect against doubled-up x values
171 //if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
172 if ( theData[high].GetX() !=0
173 //080808 TKDB
174 //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
175 &&( std::abs( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() ) < 0.000001 ) )
176 {
177 y = theData[low].GetY();
178 }
179 else
180 {
181 y = theInt.Interpolate(theManager.GetScheme(high), e,
182 theData[low].GetX(), theData[high].GetX(),
183 theData[low].GetY(), theData[high].GetY());
184 }
185 }
186 else
187 {
188 y=theData[nEntries-1].GetY();
189 }
190 return y;
191 }
192
193 void G4NeutronHPVector::Dump()
194 {
195 G4cout << nEntries<<G4endl;
196 for(G4int i=0; i<nEntries; i++)
197 {
198 G4cout << theData[i].GetX()<<" ";
199 G4cout << theData[i].GetY()<<" ";
200// if (i!=1&&i==5*(i/5)) G4cout << G4endl;
201 G4cout << G4endl;
202 }
203 G4cout << G4endl;
204 }
205
206 void G4NeutronHPVector::Check(G4int i)
207 {
208 if(i>nEntries) throw G4HadronicException(__FILE__, __LINE__, "Skipped some index numbers in G4NeutronHPVector");
209 if(i==nPoints)
210 {
211 nPoints = static_cast<G4int>(1.2*nPoints);
212 G4NeutronHPDataPoint * buff = new G4NeutronHPDataPoint[nPoints];
213 for (G4int j=0; j<nEntries; j++) buff[j] = theData[j];
214 delete [] theData;
215 theData = buff;
216 }
217 if(i==nEntries) nEntries=i+1;
218 }
219
220 void G4NeutronHPVector::
221 Merge(G4InterpolationScheme aScheme, G4double aValue,
222 G4NeutronHPVector * active, G4NeutronHPVector * passive)
223 {
224 // interpolate between labels according to aScheme, cut at aValue,
225 // continue in unknown areas by substraction of the last difference.
226
227 CleanUp();
228 G4int s = 0, n=0, m=0;
229 G4NeutronHPVector * tmp;
230 G4int a = s, p = n, t;
231 while ( a<active->GetVectorLength() )
232 {
233 if(active->GetEnergy(a) <= passive->GetEnergy(p))
234 {
235 G4double xa = active->GetEnergy(a);
236 G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
237 active->GetXsec(a), passive->GetXsec(xa));
238 SetData(m, xa, yy);
239 theManager.AppendScheme(m, active->GetScheme(a));
240 m++;
241 a++;
242 G4double xp = passive->GetEnergy(p);
243 //if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() )
244 if ( xa != 0
245 && std::abs(std::abs(xp-xa)/xa) < 0.0000001
246 && a < active->GetVectorLength() )
247 {
248 p++;
249 tmp = active; t=a;
250 active = passive; a=p;
251 passive = tmp; p=t;
252 }
253 } else {
254 tmp = active; t=a;
255 active = passive; a=p;
256 passive = tmp; p=t;
257 }
258 }
259
260 G4double deltaX = passive->GetXsec(GetEnergy(m-1)) - GetXsec(m-1);
261 while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue)
262 {
263 G4double anX;
264 anX = passive->GetXsec(p)-deltaX;
265 if(anX>0)
266 {
267 //if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
268 if ( passive->GetEnergy(p) == 0
269 || std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p) > 0.0000001 )
270 {
271 SetData(m, passive->GetEnergy(p), anX);
272 theManager.AppendScheme(m++, passive->GetScheme(p));
273 }
274 }
275 p++;
276 }
277 // Rebuild the Hash;
278 if(theHash.Prepared())
279 {
280 ReHash();
281 }
282 }
283
284 void G4NeutronHPVector::ThinOut(G4double precision)
285 {
286 // anything in there?
287 if(GetVectorLength()==0) return;
288 // make the new vector
289 G4NeutronHPDataPoint * aBuff = new G4NeutronHPDataPoint[nPoints];
290 G4double x, x1, x2, y, y1, y2;
291 G4int count = 0, current = 2, start = 1;
292
293 // First element always goes and is never tested.
294 aBuff[0] = theData[0];
295
296 // Find the rest
297 while(current < GetVectorLength())
298 {
299 x1=aBuff[count].GetX();
300 y1=aBuff[count].GetY();
301 x2=theData[current].GetX();
302 y2=theData[current].GetY();
303 for(G4int j=start; j<current; j++)
304 {
305 x = theData[j].GetX();
306 if(x1-x2 == 0) y = (y2+y1)/2.;
307 else y = theInt.Lin(x, x1, x2, y1, y2);
308 if (std::abs(y-theData[j].GetY())>precision*y)
309 {
310 aBuff[++count] = theData[current-1]; // for this one, everything was fine
311 start = current; // the next candidate
312 break;
313 }
314 }
315 current++ ;
316 }
317 // The last one also always goes, and is never tested.
318 aBuff[++count] = theData[GetVectorLength()-1];
319 delete [] theData;
320 theData = aBuff;
321 nEntries = count+1;
322
323 // Rebuild the Hash;
324 if(theHash.Prepared())
325 {
326 ReHash();
327 }
328 }
329
330 G4bool G4NeutronHPVector::IsBlocked(G4double aX)
331 {
332 G4bool result = false;
333 std::vector<G4double>::iterator i;
334 for(i=theBlocked.begin(); i!=theBlocked.end(); i++)
335 {
336 G4double aBlock = *i;
337 if(std::abs(aX-aBlock) < 0.1*MeV)
338 {
339 result = true;
340 theBlocked.erase(i);
341 break;
342 }
343 }
344 return result;
345 }
346
347 G4double G4NeutronHPVector::Sample() // Samples X according to distribution Y
348 {
349 G4double result;
350 G4int j;
351 for(j=0; j<GetVectorLength(); j++)
352 {
353 if(GetY(j)<0) SetY(j, 0);
354 }
355
356 if(theBuffered.size() !=0 && G4UniformRand()<0.5)
357 {
358 result = theBuffered[0];
359 theBuffered.erase(theBuffered.begin());
360 if(result < GetX(GetVectorLength()-1) ) return result;
361 }
362 if(GetVectorLength()==1)
363 {
364 result = theData[0].GetX();
365 }
366 else
367 {
368 if(theIntegral==0) { IntegrateAndNormalise(); }
369 do
370 {
371//080808
372/*
373 G4double rand;
374 G4double value, test, baseline;
375 baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
376 do
377 {
378 value = baseline*G4UniformRand();
379 value += theData[0].GetX();
380 test = GetY(value)/maxValue;
381 rand = G4UniformRand();
382 }
383 //while(test<rand);
384 while( test < rand && test > 0 );
385 result = value;
386*/
387 G4double rand;
388 G4double value, test;
389 do
390 {
391 rand = G4UniformRand();
392 G4int ibin = -1;
393 for ( G4int i = 0 ; i < GetVectorLength() ; i++ )
394 {
395 if ( rand < theIntegral[i] )
396 {
397 ibin = i;
398 break;
399 }
400 }
401 if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl;
402 // result
403 rand = G4UniformRand();
404 G4double x1, x2;
405 if ( ibin == 0 )
406 {
407 x1 = theData[ ibin ].GetX();
408 value = x1;
409 break;
410 }
411 else
412 {
413 x1 = theData[ ibin-1 ].GetX();
414 }
415
416 x2 = theData[ ibin ].GetX();
417 value = rand * ( x2 - x1 ) + x1;
418 test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
419 }
420 while ( G4UniformRand() > test );
421 result = value;
422//080807
423 }
424 while(IsBlocked(result));
425 }
426 return result;
427 }
428
429 G4double G4NeutronHPVector::Get15percentBorder()
430 {
431 if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
432 G4double result;
433 if(GetVectorLength()==1)
434 {
435 result = theData[0].GetX();
436 the15percentBorderCash = result;
437 }
438 else
439 {
440 if(theIntegral==0) { IntegrateAndNormalise(); }
441 G4int i;
442 result = theData[GetVectorLength()-1].GetX();
443 for(i=0;i<GetVectorLength();i++)
444 {
445 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15)
446 {
447 result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
448 the15percentBorderCash = result;
449 break;
450 }
451 }
452 the15percentBorderCash = result;
453 }
454 return result;
455 }
456
457 G4double G4NeutronHPVector::Get50percentBorder()
458 {
459 if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
460 G4double result;
461 if(GetVectorLength()==1)
462 {
463 result = theData[0].GetX();
464 the50percentBorderCash = result;
465 }
466 else
467 {
468 if(theIntegral==0) { IntegrateAndNormalise(); }
469 G4int i;
470 G4double x = 0.5;
471 result = theData[GetVectorLength()-1].GetX();
472 for(i=0;i<GetVectorLength();i++)
473 {
474 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x)
475 {
476 G4int it;
477 it = i;
478 if(it == GetVectorLength()-1)
479 {
480 result = theData[GetVectorLength()-1].GetX();
481 }
482 else
483 {
484 G4double x1, x2, y1, y2;
485 x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
486 x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
487 y1 = theData[i-1].GetX();
488 y2 = theData[i].GetX();
489 result = theLin.Lin(x, x1, x2, y1, y2);
490 }
491 the50percentBorderCash = result;
492 break;
493 }
494 }
495 the50percentBorderCash = result;
496 }
497 return result;
498 }
Note: See TracBrowser for help on using the repository browser.