| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 13-Jan-06 fix in Sample by T. Koi
31 //
32 // P. Arce, June-2014 Conversion neutron_hp to
33 //
34 #include "G4ParticleHPPartial.hh"
35
36 #include "G4ParticleHPInterpolator.hh"
37 #include "Randomize.hh"
38
39 G4ParticleHPVector* G4ParticleHPPartial::GetY(
40 {
41 auto aBuffer = new G4ParticleHPVector();
42 G4int i;
43 if (nData == 1) {
44 for (i = 0; i < data[0].GetVectorLength();
45 aBuffer->SetInterpolationManager(data[0]
46 aBuffer->SetData(i, data[0].GetX(i), dat
47 }
48 return aBuffer;
49 }
50 for (i = 0; i < nData; i++) {
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].GetVectorLen
60 x1 = data[i - 1].GetX(std::min(i1, data[i
61 x2 = data[i].GetX(ii);
62 if (x1 < x2 && i1 < data[i - 1].GetVectorL
63 y1 = data[i - 1].GetY(i1);
64 y2 = data[i].GetY(x1);
65 if (E2 - E1 != 0) {
66 y = theInt.Interpolate(theManager.GetS
67 }
68 else {
69 y = 0.5 * (y1 + y2);
70 }
71 aBuffer->SetData(ib, x1, y);
72 aBuffer->SetScheme(ib++, data[i - 1].Get
73 i1++;
74 if (x2 - x1 > 0.001 * x1) {
75 ii--;
76 }
77 }
78 else {
79 y1 = data[i - 1].GetY(x2);
80 y2 = data[i].GetY(ii);
81 if (E2 - E1 != 0) {
82 y = theInt.Interpolate(theManager.GetS
83 }
84 else {
85 y = 0.5 * (y1 + y2);
86 }
87 aBuffer->SetData(ib, x2, y);
88 aBuffer->SetScheme(ib++, data[i].GetSche
89 if (x1 - x2 < 0.001 * x2) i1++;
90 }
91 }
92 return aBuffer;
93 }
94
95 G4double G4ParticleHPPartial::Sample(G4double
96 {
97 G4int i;
98 for (i = 0; i < nData; i++) {
99 if (x < X[i]) break;
100 }
101 G4ParticleHPVector theBuff;
102 if (i == 0) {
103 theBuff.SetInterpolationManager(data[0].Ge
104 for (G4int ii = 0; ii < GetNEntries(0); ii
105 theBuff.SetX(ii, GetX(0, ii));
106 theBuff.SetY(ii, GetY(0, ii));
107 }
108 }
109 // else if(i==nData-1) this line will be del
110 else if (i == nData) {
111 for (i = 0; i < GetNEntries(nData - 1); i+
112 theBuff.SetX(i, GetX(nData - 1, i));
113 theBuff.SetY(i, GetY(nData - 1, i));
114 theBuff.SetInterpolationManager(data[nDa
115 }
116 }
117 else {
118 G4int low = i - 1;
119 G4int high = low + 1;
120 G4double x1, x2, y1, y2;
121 G4int i1 = 0, i2 = 0, ii = 0;
122 x1 = X[low];
123 x2 = X[high];
124 while (i1 < GetNEntries(low) || i2 < GetNE
125 {
126 if ((GetX(low, i1) < GetX(high, i2) && i
127 theBuff.SetX(ii, GetX(low, i1));
128 y1 = GetY(low, i1);
129 y2 = GetY(high, GetX(low, i1)); // pr
130 theBuff.SetY(ii, theInt.Interpolate(th
131 y2
132 theBuff.SetScheme(ii, data[low].GetSch
133 if (std::abs(GetX(low, i1) - GetX(high
134 i1++;
135 ii++;
136 }
137 else {
138 theBuff.SetX(ii, GetX(high, i2));
139 //************************************
140 // Change by E.Mendoza and D.Cano (CIE
141 // y1 = GetY(high,i2);
142 // y2 = GetY(low, GetX(high,i2)); //p
143 y2 = GetY(high, i2);
144 y1 = GetY(low, GetX(high, i2)); // pr
145 //************************************
146 theBuff.SetY(ii, theInt.Interpolate(th
147 y2
148 theBuff.SetScheme(ii, data[high].GetSc
149 if (std::abs(GetX(low, i1) - GetX(high
150 i2++;
151 ii++;
152 }
153 }
154 }
155 // buff is full, now sample.
156 return theBuff.Sample();
157 }
158