| Geant4 Cross Reference |
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"
35
36 #include "G4ParticleHPInterpolator.hh"
37 #include "Randomize.hh"
38
39 G4ParticleHPVector* G4ParticleHPPartial::GetY(G4double e1)
40 {
41 auto aBuffer = new G4ParticleHPVector();
42 G4int i;
43 if (nData == 1) {
44 for (i = 0; i < data[0].GetVectorLength(); i++) {
45 aBuffer->SetInterpolationManager(data[0].GetInterpolationManager());
46 aBuffer->SetData(i, data[0].GetX(i), data[0].GetY(i));
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].GetVectorLength(); ii++) {
60 x1 = data[i - 1].GetX(std::min(i1, data[i - 1].GetVectorLength() - 1));
61 x2 = data[i].GetX(ii);
62 if (x1 < x2 && i1 < data[i - 1].GetVectorLength()) {
63 y1 = data[i - 1].GetY(i1);
64 y2 = data[i].GetY(x1);
65 if (E2 - E1 != 0) {
66 y = theInt.Interpolate(theManager.GetScheme(i), e1, E1, E2, y1, y2);
67 }
68 else {
69 y = 0.5 * (y1 + y2);
70 }
71 aBuffer->SetData(ib, x1, y);
72 aBuffer->SetScheme(ib++, data[i - 1].GetScheme(i1));
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.GetScheme(i), e1, E1, E2, y1, y2);
83 }
84 else {
85 y = 0.5 * (y1 + y2);
86 }
87 aBuffer->SetData(ib, x2, y);
88 aBuffer->SetScheme(ib++, data[i].GetScheme(ii));
89 if (x1 - x2 < 0.001 * x2) i1++;
90 }
91 }
92 return aBuffer;
93 }
94
95 G4double G4ParticleHPPartial::Sample(G4double x)
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].GetInterpolationManager());
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 delete
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[nData - 1].GetInterpolationManager());
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 < GetNEntries(high)) // Loop checking, 11.05.2015, T. Koi
125 {
126 if ((GetX(low, i1) < GetX(high, i2) && i1 < GetNEntries(low)) || (i2 == GetNEntries(high))) {
127 theBuff.SetX(ii, GetX(low, i1));
128 y1 = GetY(low, i1);
129 y2 = GetY(high, GetX(low, i1)); // prob at ident theta
130 theBuff.SetY(ii, theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, y1,
131 y2)); // energy interpol
132 theBuff.SetScheme(ii, data[low].GetScheme(i1));
133 if (std::abs(GetX(low, i1) - GetX(high, i2)) < 0.001) i2++;
134 i1++;
135 ii++;
136 }
137 else {
138 theBuff.SetX(ii, GetX(high, i2));
139 //*******************************************************************************
140 // Change by E.Mendoza and D.Cano (CIEMAT):
141 // y1 = GetY(high,i2);
142 // y2 = GetY(low, GetX(high,i2)); //prob at ident theta
143 y2 = GetY(high, i2);
144 y1 = GetY(low, GetX(high, i2)); // prob at ident theta
145 //*******************************************************************************
146 theBuff.SetY(ii, theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, y1,
147 y2)); // energy interpol
148 theBuff.SetScheme(ii, data[high].GetScheme(i2));
149 if (std::abs(GetX(low, i1) - GetX(high, i2)) < 0.001) i1++;
150 i2++;
151 ii++;
152 }
153 }
154 }
155 // buff is full, now sample.
156 return theBuff.Sample();
157 }
158