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 // 080612 SampleDiscreteTwoBody contribution f 31 // #3 32 // 33 // P. Arce, June-2014 Conversion neutron_hp to 34 // 35 #include "G4ParticleHPLegendreStore.hh" 36 37 #include "G4ParticleHPFastLegendre.hh" 38 #include "G4ParticleHPInterpolator.hh" 39 #include "G4ParticleHPVector.hh" 40 #include "Randomize.hh" 41 42 #include <iostream> 43 44 // 080612TK contribution from Benoit Pirard an 45 G4double G4ParticleHPLegendreStore::SampleDisc 46 { 47 G4double result = 0.; 48 49 G4int i0; 50 G4int low(0), high(0); 51 G4ParticleHPFastLegendre theLeg; 52 for (i0 = 0; i0 < nEnergy; i0++) { 53 high = i0; 54 if (theCoeff[i0].GetEnergy() > anEnergy) b 55 } 56 low = std::max(0, high - 1); 57 G4ParticleHPInterpolator theInt; 58 G4double x, x1, x2; 59 x = anEnergy; 60 x1 = theCoeff[low].GetEnergy(); 61 x2 = theCoeff[high].GetEnergy(); 62 G4double theNorm = 0; 63 G4double try01 = 0, try02 = 0; 64 G4double max1, max2, costh; 65 max1 = 0; 66 max2 = 0; 67 G4int l, m_tmp; 68 for (i0 = 0; i0 < 601; i0++) { 69 costh = G4double(i0 - 300) / 300.; 70 try01 = 0.5; 71 for (m_tmp = 0; m_tmp < theCoeff[low].GetN 72 l = m_tmp + 1; 73 try01 += (2. * l + 1) / 2. * theCoeff[lo 74 } 75 if (try01 > max1) max1 = try01; 76 try02 = 0.5; 77 for (m_tmp = 0; m_tmp < theCoeff[high].Get 78 l = m_tmp + 1; 79 try02 += (2. * l + 1) / 2. * theCoeff[hi 80 } 81 if (try02 > max2) max2 = try02; 82 } 83 theNorm = theInt.Interpolate(theManager.GetS 84 85 G4double value, random; 86 G4double v1, v2; 87 G4int icounter = 0; 88 G4int icounter_max = 1024; 89 do { 90 icounter++; 91 if (icounter > icounter_max) { 92 G4cout << "Loop-counter exceeded the thr 93 << __FILE__ << "." << G4endl; 94 break; 95 } 96 v1 = 0.5; 97 v2 = 0.5; 98 result = 2. * G4UniformRand() - 1.; 99 for (m_tmp = 0; m_tmp < theCoeff[low].GetN 100 l = m_tmp + 1; 101 G4double legend = theLeg.Evaluate(l, res 102 v1 += (2. * l + 1) / 2. * theCoeff[low]. 103 } 104 for (m_tmp = 0; m_tmp < theCoeff[high].Get 105 l = m_tmp + 1; 106 G4double legend = theLeg.Evaluate(l, res 107 v2 += (2. * l + 1) / 2. * theCoeff[high] 108 } 109 // v1 = std::max(0.,v1); // Workaround in 110 // v2 = std::max(0.,v2); 111 value = theInt.Interpolate(theManager.GetS 112 random = G4UniformRand(); 113 if (0 >= theNorm) break; // Workaround fo 114 } while (random > value / theNorm); // Loop 115 116 return result; 117 } 118 119 G4double G4ParticleHPLegendreStore::SampleMax( 120 { 121 G4double result = 0.; 122 123 G4int i0; 124 G4int low(0), high(0); 125 G4ParticleHPFastLegendre theLeg; 126 for (i0 = 0; i0 < nEnergy; i0++) { 127 high = i0; 128 if (theCoeff[i0].GetEnergy() > anEnergy) b 129 } 130 low = std::max(0, high - 1); 131 G4ParticleHPInterpolator theInt; 132 G4double x, x1, x2; 133 x = anEnergy; 134 x1 = theCoeff[low].GetEnergy(); 135 x2 = theCoeff[high].GetEnergy(); 136 G4double theNorm = 0; 137 G4double try01 = 0, try02 = 0; 138 G4double max1, max2, costh; 139 max1 = 0; 140 max2 = 0; 141 G4int l; 142 for (i0 = 0; i0 < 601; i0++) { 143 costh = G4double(i0 - 300) / 300.; 144 try01 = 0; 145 for (l = 0; l < theCoeff[low].GetNumberOfP 146 try01 += (2. * l + 1) / 2. * theCoeff[lo 147 } 148 if (try01 > max1) max1 = try01; 149 try02 = 0; 150 for (l = 0; l < theCoeff[high].GetNumberOf 151 try02 += (2. * l + 1) / 2. * theCoeff[hi 152 } 153 if (try02 > max2) max2 = try02; 154 } 155 theNorm = theInt.Interpolate(theManager.GetS 156 157 G4double value, random; 158 G4double v1, v2; 159 G4int icounter = 0; 160 G4int icounter_max = 1024; 161 do { 162 icounter++; 163 if (icounter > icounter_max) { 164 G4cout << "Loop-counter exceeded the thr 165 << __FILE__ << "." << G4endl; 166 break; 167 } 168 v1 = 0; 169 v2 = 0; 170 result = 2. * G4UniformRand() - 1.; 171 for (l = 0; l < theCoeff[low].GetNumberOfP 172 G4double legend = theLeg.Evaluate(l, res 173 v1 += (2. * l + 1) / 2. * theCoeff[low]. 174 } 175 for (l = 0; l < theCoeff[high].GetNumberOf 176 G4double legend = theLeg.Evaluate(l, res 177 v2 += (2. * l + 1) / 2. * theCoeff[high] 178 } 179 v1 = std::max(0., v1); // Workaround in c 180 v2 = std::max(0., v2); 181 value = theInt.Interpolate(theManager.GetS 182 random = G4UniformRand(); 183 if (0 >= theNorm) break; // Workaround fo 184 } while (random > value / theNorm); // Loop 185 return result; 186 } 187 188 G4double G4ParticleHPLegendreStore::SampleElas 189 { 190 G4double result = 0.; 191 192 G4int i0; 193 G4int low(0), high(0); 194 G4ParticleHPFastLegendre theLeg; 195 for (i0 = 0; i0 < nEnergy; i0++) { 196 high = i0; 197 if (theCoeff[i0].GetEnergy() > anEnergy) b 198 } 199 low = std::max(0, high - 1); 200 G4ParticleHPInterpolator theInt; 201 G4double x, x1, x2; 202 x = anEnergy; 203 x1 = theCoeff[low].GetEnergy(); 204 x2 = theCoeff[high].GetEnergy(); 205 G4double theNorm = 0; 206 G4double try01 = 0, try02 = 0, try11 = 0, tr 207 G4double try1, try2; 208 G4int l; 209 for (l = 0; l < theCoeff[low].GetNumberOfPol 210 try01 += (2. * l + 1) / 2. * theCoeff[low] 211 try11 += (2. * l + 1) / 2. * theCoeff[low] 212 } 213 for (l = 0; l < theCoeff[high].GetNumberOfPo 214 try02 += (2. * l + 1) / 2. * theCoeff[high 215 try12 += (2. * l + 1) / 2. * theCoeff[high 216 } 217 try1 = theInt.Interpolate(theManager.GetSche 218 try2 = theInt.Interpolate(theManager.GetSche 219 theNorm = std::max(try1, try2); 220 221 G4double value, random; 222 G4double v1, v2; 223 G4int icounter = 0; 224 G4int icounter_max = 1024; 225 do { 226 icounter++; 227 if (icounter > icounter_max) { 228 G4cout << "Loop-counter exceeded the thr 229 << __FILE__ << "." << G4endl; 230 break; 231 } 232 v1 = 0; 233 v2 = 0; 234 result = 2. * G4UniformRand() - 1.; 235 for (l = 0; l < theCoeff[low].GetNumberOfP 236 G4double legend = theLeg.Evaluate(l, res 237 v1 += (2. * l + 1) / 2. * theCoeff[low]. 238 } 239 for (l = 0; l < theCoeff[high].GetNumberOf 240 G4double legend = theLeg.Evaluate(l, res 241 v2 += (2. * l + 1) / 2. * theCoeff[high] 242 } 243 value = theInt.Interpolate(theManager.GetS 244 random = G4UniformRand(); 245 } while (random > value / theNorm); // Loop 246 247 return result; 248 } 249 250 G4double G4ParticleHPLegendreStore::Sample(G4d 251 { 252 G4int i0; 253 G4int low(0), high(0); 254 // G4cout << "G4ParticleHPLegendreStore::Sa 255 for (i0 = 0; i0 < nEnergy; i0++) { 256 // G4cout <<"theCoeff["<<i0<<"].GetEne 257 high = i0; 258 if (theCoeff[i0].GetEnergy() > energy) bre 259 } 260 low = std::max(0, high - 1); 261 // G4cout << "G4ParticleHPLegendreStore::Sa 262 G4ParticleHPVector theBuffer; 263 G4ParticleHPInterpolator theInt; 264 G4double x1, x2, y1, y2, y; 265 x1 = theCoeff[low].GetEnergy(); 266 x2 = theCoeff[high].GetEnergy(); 267 // G4cout << "the xes "<<x1<<" "<<x2<<G4end 268 G4double costh = 0; 269 for (i0 = 0; i0 < 601; i0++) { 270 costh = G4double(i0 - 300) / 300.; 271 y1 = Integrate(low, costh); 272 y2 = Integrate(high, costh); 273 y = theInt.Interpolate(theManager.GetSchem 274 theBuffer.SetData(i0, costh, y); 275 // G4cout << "Integration "<<low<<" "< 276 } 277 G4double rand = G4UniformRand(); 278 G4int it; 279 for (i0 = 1; i0 < 601; i0++) { 280 it = i0; 281 if (rand < theBuffer.GetY(i0) / theBuffer. 282 // G4cout <<"sampling now "<<i0<<" " 283 // << theBuffer.GetY(i0)<<" " 284 // << theBuffer.GetY(600)<<" " 285 // << rand<<" " 286 // << theBuffer.GetY(i0)/theBuffer 287 } 288 if (it == 601) it = 600; 289 // G4cout << "G4ParticleHPLegendreStore::Sa 290 G4double norm = theBuffer.GetY(600); 291 if (norm == 0) return -DBL_MAX; 292 x1 = theBuffer.GetY(it) / norm; 293 x2 = theBuffer.GetY(it - 1) / norm; 294 y1 = theBuffer.GetX(it); 295 y2 = theBuffer.GetX(it - 1); 296 // G4cout << "G4ParticleHPLegendreStore::Sa 297 return theInt.Interpolate(theManager.GetSche 298 } 299 300 G4double 301 G4ParticleHPLegendreStore::Integrate(G4int k, 302 G4double 303 { 304 G4double result = 0.; 305 G4ParticleHPFastLegendre theLeg; 306 // G4cout <<"the COEFFS "<<k<<" "; 307 // G4cout <<theCoeff[k].GetNumberOfPoly()<< 308 for (G4int l = 0; l < theCoeff[k].GetNumberO 309 result += theCoeff[k].GetCoeff(l) * theLeg 310 // G4cout << theCoeff[k].GetCoeff(l)<<" 311 } 312 // G4cout <<G4endl; 313 return result; 314 } 315