Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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  }
void SetData(G4int i, G4double x, G4double y)
void SetInterpolationManager(const G4InterpolationManager &aManager)
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
const XML_Char const XML_Char * data
Definition: expat.h:268
G4double GetX(G4int i)
void SetY(G4int i, G4double x)
G4double Sample(G4double x)
G4InterpolationScheme GetScheme(G4int index) const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
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)