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 // 27 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 28 // 29 #ifndef G4ParticleHPFastLegendre_h 30 #define G4ParticleHPFastLegendre_h 1 31 32 #include "globals.hh" 33 34 class G4ParticleHPFastLegendre 35 { 36 public: 37 G4ParticleHPFastLegendre() 38 { 39 value = new const G4double*[31]; 40 value[0] = l0; 41 value[1] = l1; 42 value[2] = l2; 43 value[3] = l3; 44 value[4] = l4; 45 value[5] = l5; 46 value[6] = l6; 47 value[7] = l7; 48 value[8] = l8; 49 value[9] = l9; 50 value[10] = l10; 51 value[11] = l11; 52 value[12] = l12; 53 value[13] = l13; 54 value[14] = l14; 55 value[15] = l15; 56 value[16] = l16; 57 value[17] = l17; 58 value[18] = l18; 59 value[19] = l19; 60 value[20] = l20; 61 value[21] = l21; 62 value[22] = l22; 63 value[23] = l23; 64 value[24] = l24; 65 value[25] = l25; 66 value[26] = l26; 67 value[27] = l27; 68 value[28] = l28; 69 value[29] = l29; 70 value[30] = l30; 71 integral = new const G4double*[31]; 72 integral[0] = i0; 73 integral[1] = i1; 74 integral[2] = i2; 75 integral[3] = i3; 76 integral[4] = i4; 77 integral[5] = i5; 78 integral[6] = i6; 79 integral[7] = i7; 80 integral[8] = i8; 81 integral[9] = i9; 82 integral[10] = i10; 83 integral[11] = i11; 84 integral[12] = i12; 85 integral[13] = i13; 86 integral[14] = i14; 87 integral[15] = i15; 88 integral[16] = i16; 89 integral[17] = i17; 90 integral[18] = i18; 91 integral[19] = i19; 92 integral[20] = i20; 93 integral[21] = i21; 94 integral[22] = i22; 95 integral[23] = i23; 96 integral[24] = i24; 97 integral[25] = i25; 98 integral[26] = i26; 99 integral[27] = i27; 100 integral[28] = i28; 101 integral[29] = i29; 102 integral[30] = i30; 103 104 G4int i; 105 for (i = 0; i < 31; i++) 106 theNbin[i] = 1 + 200 * (i + 1); 107 } 108 109 ~G4ParticleHPFastLegendre() 110 { 111 delete[] value; 112 delete[] integral; 113 } 114 115 G4double Integrate(G4int l, G4double costh) 116 { 117 if (l > 30) return regularIntegrate(l, costh); 118 G4int bin = GetBin(l, costh); 119 G4double y1, y2; 120 // G4cout <<"Testhpw G4ParticleHPFastLegendre::Integrate "<<l<<" "<<bin<<G4endl; 121 y1 = integral[l][bin]; 122 y2 = integral[l][bin + 1]; 123 // G4cout <<"Testhpw G4ParticleHPFastLegendre::Integrate exit"<<G4endl; 124 return Interpolate(bin, l, y1, y2, costh); 125 } 126 127 inline G4double Evaluate(G4int l, G4double costh) 128 { 129 if (l > 30) return regularEvaluate(l, costh); 130 G4double result; 131 G4int bin = GetBin(l, costh); 132 if (bin != theNbin[l] - 1) { 133 G4double y1, y2; 134 y1 = value[l][bin]; 135 y2 = value[l][bin + 1]; 136 result = Interpolate(bin, l, y1, y2, costh); 137 } 138 else { 139 result = value[l][bin]; 140 } 141 return result; 142 } 143 144 private: 145 G4double regularEvaluate(int l, double x); 146 G4double regularIntegrate(int l, double x); 147 148 inline G4int GetBin(G4int l, G4double costh) 149 { 150 G4int bin = 0; 151 bin = G4int((theNbin[l] - 1) * (costh + 1) / 2.); 152 if (bin == theNbin[l] - 1) bin--; 153 return bin; 154 } 155 156 inline G4double Interpolate(G4int bin, G4int l, G4double y1, G4double y2, G4double x) 157 { 158 G4double slope = 0, off = 0, x2 = 0, x1mx2; 159 G4int half = (theNbin[l] - 1) / 2; 160 // x1 = (bin-half)/G4double(half); 161 x2 = (bin + 1 - half) / G4double(half); 162 x1mx2 = 1. / G4double((theNbin[l] - 1) / 2); 163 // slope = (y2-y1)/(x2-x1); 164 slope = (y2 - y1) / x1mx2; 165 off = y2 - x2 * slope; 166 return x * slope + off; 167 } 168 169 const G4double** value; 170 const G4double** integral; 171 G4int theNbin[31]; 172 static const G4double l0[201]; 173 static const G4double i0[201]; 174 static const G4double l1[401]; 175 static const G4double i1[401]; 176 static const G4double l2[601]; 177 static const G4double i2[601]; 178 static const G4double l3[801]; 179 static const G4double i3[801]; 180 static const G4double l4[1001]; 181 static const G4double i4[1001]; 182 static const G4double l5[1201]; 183 static const G4double i5[1201]; 184 static const G4double l6[1401]; 185 static const G4double i6[1401]; 186 static const G4double l7[1601]; 187 static const G4double i7[1601]; 188 static const G4double l8[1801]; 189 static const G4double i8[1801]; 190 static const G4double l9[2001]; 191 static const G4double i9[2001]; 192 static const G4double l10[2201]; 193 static const G4double i10[2201]; 194 static const G4double l11[2401]; 195 static const G4double i11[2401]; 196 static const G4double l12[2601]; 197 static const G4double i12[2601]; 198 static const G4double l13[2801]; 199 static const G4double i13[2801]; 200 static const G4double l14[3001]; 201 static const G4double i14[3001]; 202 static const G4double l15[3201]; 203 static const G4double i15[3201]; 204 static const G4double l16[3401]; 205 static const G4double i16[3401]; 206 static const G4double l17[3601]; 207 static const G4double i17[3601]; 208 static const G4double l18[3801]; 209 static const G4double i18[3801]; 210 static const G4double l19[4001]; 211 static const G4double i19[4001]; 212 static const G4double l20[4201]; 213 static const G4double i20[4201]; 214 static const G4double l21[4401]; 215 static const G4double i21[4401]; 216 static const G4double l22[4601]; 217 static const G4double i22[4601]; 218 static const G4double l23[4801]; 219 static const G4double i23[4801]; 220 static const G4double l24[5001]; 221 static const G4double i24[5001]; 222 static const G4double l25[5201]; 223 static const G4double i25[5201]; 224 static const G4double l26[5401]; 225 static const G4double i26[5401]; 226 static const G4double l27[5601]; 227 static const G4double i27[5601]; 228 static const G4double l28[5801]; 229 static const G4double i28[5801]; 230 static const G4double l29[6001]; 231 static const G4double i29[6001]; 232 static const G4double l30[6201]; 233 static const G4double i30[6201]; 234 }; 235 #endif 236