Geant4  10.01
G4NeutronHPPartial.cc
Go to the documentation of this file.
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 // 13-Jan-06 fix in Sample by T. Koi
31 //
32 #include "G4NeutronHPPartial.hh"
34 #include "Randomize.hh"
35 
37  {
38  G4NeutronHPVector * aBuffer = new G4NeutronHPVector();
39  G4int i;
40  if(nData==1)
41  {
42  for(i=0; i<data[0].GetVectorLength(); i++)
43  {
44  aBuffer->SetInterpolationManager(data[0].GetInterpolationManager());
45  aBuffer->SetData(i , data[0].GetX(i), data[0].GetY(i));
46  }
47  return aBuffer;
48  }
49  for (i=0; i<nData; i++)
50  {
51  if(X[i]>e1) break;
52  }
53  if(i==nData) i--;
54  if(0==i) i=1;
55  G4double x1,x2,y1,y2,y;
56  G4int i1=0, ib=0;
57  G4double E1 = X[i-1];
58  G4double E2 = X[i];
59  for(G4int ii=0; ii<data[i].GetVectorLength(); ii++)
60  {
61  x1 = data[i-1].GetX(std::min(i1, data[i-1].GetVectorLength()-1));
62  x2 = data[i].GetX(ii);
63  if(x1<x2&&i1<data[i-1].GetVectorLength())
64  {
65  y1 = data[i-1].GetY(i1);
66  y2 = data[i].GetY(x1);
67  if(E2-E1!=0)
68  {
69  y = theInt.Interpolate(theManager.GetScheme(i), e1, E1, E2, y1, y2);
70  }
71  else
72  {
73  y = 0.5*(y1+y2);
74  }
75  aBuffer->SetData(ib, x1, y);
76  aBuffer->SetScheme(ib++, data[i-1].GetScheme(i1));
77  i1++;
78  if(x2-x1>0.001*x1)
79  {
80  ii--;
81  }
82  }
83  else
84  {
85  y1 = data[i-1].GetY(x2);
86  y2 = data[i].GetY(ii);
87  if(E2-E1!=0)
88  {
89  y = theInt.Interpolate(theManager.GetScheme(i), e1, E1, E2, y1, y2);
90  }
91  else
92  {
93  y = 0.5*(y1+y2);
94  }
95  aBuffer->SetData(ib, x2, y);
96  aBuffer->SetScheme(ib++, data[i].GetScheme(ii));
97  if(x1-x2<0.001*x2) i1++;
98  }
99  }
100  return aBuffer;
101  }
102 
104  {
105  G4int i;
106  for (i=0; i<nData; i++)
107  {
108  if(x<X[i]) break;
109  }
110  G4NeutronHPVector theBuff;
111  if(i==0)
112  {
113  theBuff.SetInterpolationManager(data[0].GetInterpolationManager());
114  for(G4int ii=0;ii<GetNEntries(0);ii++)
115  {
116  theBuff.SetX(ii, GetX(0,ii));
117  theBuff.SetY(ii, GetY(0,ii));
118  }
119  }
120  //else if(i==nData-1) this line will be delete
121  else if ( i == nData )
122  {
123  for(i=0;i<GetNEntries(nData-1);i++)
124  {
125  theBuff.SetX(i, GetX(nData-1,i));
126  theBuff.SetY(i, GetY(nData-1,i));
127  theBuff.SetInterpolationManager(data[nData-1].GetInterpolationManager());
128  }
129  }
130  else
131  {
132  G4int low = i-1;
133  G4int high = low+1;
134  G4double x1,x2,y1,y2;
135  G4int i1=0, i2=0, ii=0;
136  x1 = X[low];
137  x2 = X[high];
138  while(i1<GetNEntries(low)||i2<GetNEntries(high))
139  {
140  if( (GetX(low,i1)<GetX(high,i2) && i1<GetNEntries(low))
141  ||(i2==GetNEntries(high)) )
142  {
143  theBuff.SetX(ii, GetX(low,i1));
144  y1 = GetY(low,i1);
145  y2 = GetY(high, GetX(low,i1)); //prob at ident theta
146  theBuff.SetY(ii, theInt.Interpolate(theManager.GetScheme(high),
147  x, x1, x2, y1, y2)); //energy interpol
148  theBuff.SetScheme(ii, data[low].GetScheme(i1));
149  if(std::abs(GetX(low,i1)-GetX(high,i2))<0.001) i2++;
150  i1++;
151  ii++;
152  }
153  else
154  {
155  theBuff.SetX(ii, GetX(high,i2));
156  //*******************************************************************************
157  //Change by E.Mendoza and D.Cano (CIEMAT):
158  //y1 = GetY(high,i2);
159  //y2 = GetY(low, GetX(high,i2)); //prob at ident theta
160  y2 = GetY(high,i2);
161  y1 = GetY(low, GetX(high,i2)); //prob at ident theta
162  //*******************************************************************************
163  theBuff.SetY(ii, theInt.Interpolate(theManager.GetScheme(high),
164  x, x1, x2, y1, y2)); //energy interpol
165  theBuff.SetScheme(ii, data[high].GetScheme(i2));
166  if(std::abs(GetX(low,i1)-GetX(high,i2))<0.001) i1++;
167  i2++;
168  ii++;
169  }
170  }
171  }
172  //buff is full, now sample.
173  return theBuff.Sample();
174  }
G4double GetY(G4double x)
G4int GetVectorLength() const
void SetInterpolationManager(const G4InterpolationManager &aManager)
G4double GetX(G4int i) const
void SetData(G4int i, G4double x, G4double y)
int G4int
Definition: G4Types.hh:78
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4NeutronHPInterpolator theInt
void SetY(G4int i, G4double x)
G4NeutronHPVector * data
G4int GetNEntries(G4int i)
G4InterpolationScheme GetScheme(G4int index) const
void SetX(G4int i, G4double e)
static const G4double e1
G4double GetY(G4int i, G4int j)
void SetScheme(G4int aPoint, const G4InterpolationScheme &aScheme)
G4double GetX(G4int i)
G4InterpolationManager theManager
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double G4double
Definition: G4Types.hh:76
G4double Sample(G4double x)