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