Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // neutron_hp -- source file 26 // neutron_hp -- source file 27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996 28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model. 29 // 29 // 30 // 080612 SampleDiscreteTwoBody contribution f << 30 // 080612 SampleDiscreteTwoBody contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3 31 // #3 << 32 // 31 // 33 // P. Arce, June-2014 Conversion neutron_hp to 32 // P. Arce, June-2014 Conversion neutron_hp to particle_hp 34 // 33 // 35 #include "G4ParticleHPLegendreStore.hh" 34 #include "G4ParticleHPLegendreStore.hh" 36 << 37 #include "G4ParticleHPFastLegendre.hh" << 38 #include "G4ParticleHPInterpolator.hh" << 39 #include "G4ParticleHPVector.hh" 35 #include "G4ParticleHPVector.hh" >> 36 #include "G4ParticleHPInterpolator.hh" >> 37 #include "G4ParticleHPFastLegendre.hh" 40 #include "Randomize.hh" 38 #include "Randomize.hh" 41 << 42 #include <iostream> 39 #include <iostream> 43 40 44 // 080612TK contribution from Benoit Pirard an << 45 G4double G4ParticleHPLegendreStore::SampleDisc << 46 { << 47 G4double result = 0.; << 48 41 >> 42 >> 43 //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3 >> 44 G4double G4ParticleHPLegendreStore::SampleDiscreteTwoBody (G4double anEnergy) >> 45 { >> 46 G4double result; >> 47 49 G4int i0; 48 G4int i0; 50 G4int low(0), high(0); 49 G4int low(0), high(0); 51 G4ParticleHPFastLegendre theLeg; 50 G4ParticleHPFastLegendre theLeg; 52 for (i0 = 0; i0 < nEnergy; i0++) { << 51 for (i0=0; i0<nEnergy; i0++) >> 52 { 53 high = i0; 53 high = i0; 54 if (theCoeff[i0].GetEnergy() > anEnergy) b << 54 if(theCoeff[i0].GetEnergy()>anEnergy) break; 55 } 55 } 56 low = std::max(0, high - 1); << 56 low = std::max(0, high-1); 57 G4ParticleHPInterpolator theInt; 57 G4ParticleHPInterpolator theInt; 58 G4double x, x1, x2; 58 G4double x, x1, x2; 59 x = anEnergy; 59 x = anEnergy; 60 x1 = theCoeff[low].GetEnergy(); 60 x1 = theCoeff[low].GetEnergy(); 61 x2 = theCoeff[high].GetEnergy(); 61 x2 = theCoeff[high].GetEnergy(); 62 G4double theNorm = 0; 62 G4double theNorm = 0; 63 G4double try01 = 0, try02 = 0; << 63 G4double try01=0, try02=0; 64 G4double max1, max2, costh; 64 G4double max1, max2, costh; 65 max1 = 0; << 65 max1 = 0; max2 = 0; 66 max2 = 0; << 66 G4int l,m_tmp; 67 G4int l, m_tmp; << 67 for(i0=0; i0<601; i0++) 68 for (i0 = 0; i0 < 601; i0++) { << 68 { 69 costh = G4double(i0 - 300) / 300.; << 69 costh = G4double(i0-300)/300.; 70 try01 = 0.5; << 70 try01 = 0.5; 71 for (m_tmp = 0; m_tmp < theCoeff[low].GetN << 71 for(m_tmp=0; m_tmp<theCoeff[low].GetNumberOfPoly() ; m_tmp++) 72 l = m_tmp + 1; << 72 { 73 try01 += (2. * l + 1) / 2. * theCoeff[lo << 73 l=m_tmp+1; 74 } << 74 try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(m_tmp)*theLeg.Evaluate(l, costh); 75 if (try01 > max1) max1 = try01; << 75 } 76 try02 = 0.5; << 76 if(try01>max1) max1=try01; 77 for (m_tmp = 0; m_tmp < theCoeff[high].Get << 77 try02 = 0.5; 78 l = m_tmp + 1; << 78 for(m_tmp=0; m_tmp<theCoeff[high].GetNumberOfPoly() ; m_tmp++) 79 try02 += (2. * l + 1) / 2. * theCoeff[hi << 79 { 80 } << 80 l=m_tmp+1; 81 if (try02 > max2) max2 = try02; << 81 try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*theLeg.Evaluate(l, costh); 82 } << 82 } >> 83 if(try02>max2) max2=try02; >> 84 } 83 theNorm = theInt.Interpolate(theManager.GetS 85 theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2); 84 << 86 85 G4double value, random; 87 G4double value, random; 86 G4double v1, v2; 88 G4double v1, v2; 87 G4int icounter = 0; << 89 do 88 G4int icounter_max = 1024; << 90 { 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; 91 v1 = 0.5; 97 v2 = 0.5; 92 v2 = 0.5; 98 result = 2. * G4UniformRand() - 1.; << 93 result = 2.*G4UniformRand()-1.; 99 for (m_tmp = 0; m_tmp < theCoeff[low].GetN << 94 for(m_tmp=0; m_tmp<theCoeff[low].GetNumberOfPoly() ; m_tmp++) 100 l = m_tmp + 1; << 95 { 101 G4double legend = theLeg.Evaluate(l, res << 96 l=m_tmp+1; 102 v1 += (2. * l + 1) / 2. * theCoeff[low]. << 97 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN 103 } << 98 v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(m_tmp)*legend; 104 for (m_tmp = 0; m_tmp < theCoeff[high].Get << 99 } 105 l = m_tmp + 1; << 100 for(m_tmp=0; m_tmp<theCoeff[high].GetNumberOfPoly() ; m_tmp++) 106 G4double legend = theLeg.Evaluate(l, res << 101 { 107 v2 += (2. * l + 1) / 2. * theCoeff[high] << 102 l=m_tmp+1; 108 } << 103 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN >> 104 v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(m_tmp)*legend; >> 105 } 109 // v1 = std::max(0.,v1); // Workaround in 106 // v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical. 110 // v2 = std::max(0.,v2); << 107 // v2 = std::max(0.,v2); 111 value = theInt.Interpolate(theManager.GetS 108 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2); 112 random = G4UniformRand(); 109 random = G4UniformRand(); 113 if (0 >= theNorm) break; // Workaround fo << 110 if(0>=theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000 114 } while (random > value / theNorm); // Loop << 111 } >> 112 while(random>value/theNorm); 115 113 116 return result; 114 return result; 117 } 115 } 118 116 119 G4double G4ParticleHPLegendreStore::SampleMax( << 120 { << 121 G4double result = 0.; << 122 117 >> 118 >> 119 G4double G4ParticleHPLegendreStore::SampleMax (G4double anEnergy) >> 120 { >> 121 G4double result; >> 122 123 G4int i0; 123 G4int i0; 124 G4int low(0), high(0); 124 G4int low(0), high(0); 125 G4ParticleHPFastLegendre theLeg; 125 G4ParticleHPFastLegendre theLeg; 126 for (i0 = 0; i0 < nEnergy; i0++) { << 126 for (i0=0; i0<nEnergy; i0++) >> 127 { 127 high = i0; 128 high = i0; 128 if (theCoeff[i0].GetEnergy() > anEnergy) b << 129 if(theCoeff[i0].GetEnergy()>anEnergy) break; 129 } 130 } 130 low = std::max(0, high - 1); << 131 low = std::max(0, high-1); 131 G4ParticleHPInterpolator theInt; 132 G4ParticleHPInterpolator theInt; 132 G4double x, x1, x2; 133 G4double x, x1, x2; 133 x = anEnergy; 134 x = anEnergy; 134 x1 = theCoeff[low].GetEnergy(); 135 x1 = theCoeff[low].GetEnergy(); 135 x2 = theCoeff[high].GetEnergy(); 136 x2 = theCoeff[high].GetEnergy(); 136 G4double theNorm = 0; 137 G4double theNorm = 0; 137 G4double try01 = 0, try02 = 0; << 138 G4double try01=0, try02=0; 138 G4double max1, max2, costh; 139 G4double max1, max2, costh; 139 max1 = 0; << 140 max1 = 0; max2 = 0; 140 max2 = 0; << 141 G4int l; 141 G4int l; 142 for (i0 = 0; i0 < 601; i0++) { << 142 for(i0=0; i0<601; i0++) 143 costh = G4double(i0 - 300) / 300.; << 143 { >> 144 costh = G4double(i0-300)/300.; 144 try01 = 0; 145 try01 = 0; 145 for (l = 0; l < theCoeff[low].GetNumberOfP << 146 for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++) 146 try01 += (2. * l + 1) / 2. * theCoeff[lo << 147 { 147 } << 148 try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, costh); 148 if (try01 > max1) max1 = try01; << 149 } >> 150 if(try01>max1) max1=try01; 149 try02 = 0; 151 try02 = 0; 150 for (l = 0; l < theCoeff[high].GetNumberOf << 152 for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++) 151 try02 += (2. * l + 1) / 2. * theCoeff[hi << 153 { >> 154 try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, costh); 152 } 155 } 153 if (try02 > max2) max2 = try02; << 156 if(try02>max2) max2=try02; 154 } << 157 } 155 theNorm = theInt.Interpolate(theManager.GetS 158 theNorm = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, max1, max2); 156 << 159 157 G4double value, random; 160 G4double value, random; 158 G4double v1, v2; 161 G4double v1, v2; 159 G4int icounter = 0; << 162 do 160 G4int icounter_max = 1024; << 163 { 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; 164 v1 = 0; 169 v2 = 0; 165 v2 = 0; 170 result = 2. * G4UniformRand() - 1.; << 166 result = 2.*G4UniformRand()-1.; 171 for (l = 0; l < theCoeff[low].GetNumberOfP << 167 for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++) 172 G4double legend = theLeg.Evaluate(l, res << 168 { 173 v1 += (2. * l + 1) / 2. * theCoeff[low]. << 169 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN 174 } << 170 v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend; 175 for (l = 0; l < theCoeff[high].GetNumberOf << 171 } 176 G4double legend = theLeg.Evaluate(l, res << 172 for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++) 177 v2 += (2. * l + 1) / 2. * theCoeff[high] << 173 { 178 } << 174 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN 179 v1 = std::max(0., v1); // Workaround in c << 175 v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend; 180 v2 = std::max(0., v2); << 176 } >> 177 v1 = std::max(0.,v1); // Workaround in case one of the distributions is fully non-physical. >> 178 v2 = std::max(0.,v2); 181 value = theInt.Interpolate(theManager.GetS 179 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2); 182 random = G4UniformRand(); 180 random = G4UniformRand(); 183 if (0 >= theNorm) break; // Workaround fo << 181 if(0>=theNorm) break; // Workaround for negative cross-section values. @@@@ 31 May 2000 184 } while (random > value / theNorm); // Loop << 182 } >> 183 while(random>value/theNorm); >> 184 185 return result; 185 return result; 186 } 186 } 187 187 188 G4double G4ParticleHPLegendreStore::SampleElas << 189 { << 190 G4double result = 0.; << 191 188 >> 189 G4double G4ParticleHPLegendreStore::SampleElastic (G4double anEnergy) >> 190 { >> 191 G4double result; >> 192 192 G4int i0; 193 G4int i0; 193 G4int low(0), high(0); 194 G4int low(0), high(0); 194 G4ParticleHPFastLegendre theLeg; 195 G4ParticleHPFastLegendre theLeg; 195 for (i0 = 0; i0 < nEnergy; i0++) { << 196 for (i0=0; i0<nEnergy; i0++) >> 197 { 196 high = i0; 198 high = i0; 197 if (theCoeff[i0].GetEnergy() > anEnergy) b << 199 if(theCoeff[i0].GetEnergy()>anEnergy) break; 198 } 200 } 199 low = std::max(0, high - 1); << 201 low = std::max(0, high-1); 200 G4ParticleHPInterpolator theInt; 202 G4ParticleHPInterpolator theInt; 201 G4double x, x1, x2; 203 G4double x, x1, x2; 202 x = anEnergy; 204 x = anEnergy; 203 x1 = theCoeff[low].GetEnergy(); 205 x1 = theCoeff[low].GetEnergy(); 204 x2 = theCoeff[high].GetEnergy(); 206 x2 = theCoeff[high].GetEnergy(); 205 G4double theNorm = 0; 207 G4double theNorm = 0; 206 G4double try01 = 0, try02 = 0, try11 = 0, tr << 208 G4double try01=0, try02=0, try11=0, try12=0; 207 G4double try1, try2; 209 G4double try1, try2; 208 G4int l; 210 G4int l; 209 for (l = 0; l < theCoeff[low].GetNumberOfPol << 211 for(l=0; l<theCoeff[low].GetNumberOfPoly(); l++) 210 try01 += (2. * l + 1) / 2. * theCoeff[low] << 212 { 211 try11 += (2. * l + 1) / 2. * theCoeff[low] << 213 try01 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, -1.); 212 } << 214 try11 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*theLeg.Evaluate(l, +1.); 213 for (l = 0; l < theCoeff[high].GetNumberOfPo << 215 } 214 try02 += (2. * l + 1) / 2. * theCoeff[high << 216 for(l=0; l<theCoeff[high].GetNumberOfPoly(); l++) 215 try12 += (2. * l + 1) / 2. * theCoeff[high << 217 { 216 } << 218 try02 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, -1.); >> 219 try12 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*theLeg.Evaluate(l, +1.); >> 220 } 217 try1 = theInt.Interpolate(theManager.GetSche 221 try1 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try01, try02); 218 try2 = theInt.Interpolate(theManager.GetSche 222 try2 = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, try11, try12); 219 theNorm = std::max(try1, try2); 223 theNorm = std::max(try1, try2); 220 << 224 221 G4double value, random; 225 G4double value, random; 222 G4double v1, v2; 226 G4double v1, v2; 223 G4int icounter = 0; << 227 do 224 G4int icounter_max = 1024; << 228 { 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; 229 v1 = 0; 233 v2 = 0; 230 v2 = 0; 234 result = 2. * G4UniformRand() - 1.; << 231 result = 2.*G4UniformRand()-1.; 235 for (l = 0; l < theCoeff[low].GetNumberOfP << 232 for(l=0; l<theCoeff[low].GetNumberOfPoly() ; l++) 236 G4double legend = theLeg.Evaluate(l, res << 233 { 237 v1 += (2. * l + 1) / 2. * theCoeff[low]. << 234 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN 238 } << 235 v1 += (2.*l+1)/2.*theCoeff[low].GetCoeff(l)*legend; 239 for (l = 0; l < theCoeff[high].GetNumberOf << 236 } 240 G4double legend = theLeg.Evaluate(l, res << 237 for(l=0; l<theCoeff[high].GetNumberOfPoly() ; l++) 241 v2 += (2. * l + 1) / 2. * theCoeff[high] << 238 { 242 } << 239 G4double legend = theLeg.Evaluate(l, result); // @@@ done to avoid optimization error on SUN >> 240 v2 += (2.*l+1)/2.*theCoeff[high].GetCoeff(l)*legend; >> 241 } 243 value = theInt.Interpolate(theManager.GetS 242 value = theInt.Interpolate(theManager.GetScheme(high), x, x1, x2, v1, v2); 244 random = G4UniformRand(); 243 random = G4UniformRand(); 245 } while (random > value / theNorm); // Loop << 244 } 246 << 245 while(random>value/theNorm); >> 246 247 return result; 247 return result; 248 } 248 } 249 249 250 G4double G4ParticleHPLegendreStore::Sample(G4d << 250 G4double G4ParticleHPLegendreStore::Sample (G4double energy) // still in interpolation; do not use 251 { 251 { 252 G4int i0; 252 G4int i0; 253 G4int low(0), high(0); 253 G4int low(0), high(0); 254 // G4cout << "G4ParticleHPLegendreStore::Sa << 254 // G4cout << "G4ParticleHPLegendreStore::Sample "<<energy<<" "<<energy<<" "<<nEnergy<<G4endl; 255 for (i0 = 0; i0 < nEnergy; i0++) { << 255 for (i0=0; i0<nEnergy; i0++) 256 // G4cout <<"theCoeff["<<i0<<"].GetEne << 256 { >> 257 // G4cout <<"theCoeff["<<i0<<"].GetEnergy() = "<<theCoeff[i0].GetEnergy()<<G4endl; 257 high = i0; 258 high = i0; 258 if (theCoeff[i0].GetEnergy() > energy) bre << 259 if(theCoeff[i0].GetEnergy()>energy) break; 259 } 260 } 260 low = std::max(0, high - 1); << 261 low = std::max(0, high-1); 261 // G4cout << "G4ParticleHPLegendreStore::Sa << 262 // G4cout << "G4ParticleHPLegendreStore::Sample high, low: "<<high<<", "<<low<<G4endl; 262 G4ParticleHPVector theBuffer; 263 G4ParticleHPVector theBuffer; 263 G4ParticleHPInterpolator theInt; 264 G4ParticleHPInterpolator theInt; 264 G4double x1, x2, y1, y2, y; 265 G4double x1, x2, y1, y2, y; 265 x1 = theCoeff[low].GetEnergy(); 266 x1 = theCoeff[low].GetEnergy(); 266 x2 = theCoeff[high].GetEnergy(); 267 x2 = theCoeff[high].GetEnergy(); 267 // G4cout << "the xes "<<x1<<" "<<x2<<G4end << 268 // G4cout << "the xes "<<x1<<" "<<x2<<G4endl; 268 G4double costh = 0; << 269 G4double costh=0; 269 for (i0 = 0; i0 < 601; i0++) { << 270 for(i0=0; i0<601; i0++) 270 costh = G4double(i0 - 300) / 300.; << 271 { >> 272 costh = G4double(i0-300)/300.; 271 y1 = Integrate(low, costh); 273 y1 = Integrate(low, costh); 272 y2 = Integrate(high, costh); 274 y2 = Integrate(high, costh); 273 y = theInt.Interpolate(theManager.GetSchem 275 y = theInt.Interpolate(theManager.GetScheme(high), energy, x1, x2, y1, y2); 274 theBuffer.SetData(i0, costh, y); 276 theBuffer.SetData(i0, costh, y); 275 // G4cout << "Integration "<<low<<" "< << 277 // G4cout << "Integration "<<low<<" "<<costh<<" "<<y1<<" "<<y2<<" "<<y<<G4endl; 276 } 278 } 277 G4double rand = G4UniformRand(); 279 G4double rand = G4UniformRand(); 278 G4int it; 280 G4int it; 279 for (i0 = 1; i0 < 601; i0++) { << 281 for (i0=1; i0<601; i0++) >> 282 { 280 it = i0; 283 it = i0; 281 if (rand < theBuffer.GetY(i0) / theBuffer. << 284 if(rand < theBuffer.GetY(i0)/theBuffer.GetY(600)) break; 282 // G4cout <<"sampling now "<<i0<<" " << 285 // G4cout <<"sampling now "<<i0<<" " 283 // << theBuffer.GetY(i0)<<" " << 286 // << theBuffer.GetY(i0)<<" " 284 // << theBuffer.GetY(600)<<" " << 287 // << theBuffer.GetY(600)<<" " 285 // << rand<<" " << 288 // << rand<<" " 286 // << theBuffer.GetY(i0)/theBuffer << 289 // << theBuffer.GetY(i0)/theBuffer.GetY(600)<<G4endl;; 287 } 290 } 288 if (it == 601) it = 600; << 291 if(it==601) it=600; 289 // G4cout << "G4ParticleHPLegendreStore::Sa << 292 // G4cout << "G4ParticleHPLegendreStore::Sample it "<<rand<<" "<<it<<G4endl; 290 G4double norm = theBuffer.GetY(600); 293 G4double norm = theBuffer.GetY(600); 291 if (norm == 0) return -DBL_MAX; << 294 if(norm==0) return -DBL_MAX; 292 x1 = theBuffer.GetY(it) / norm; << 295 x1 = theBuffer.GetY(it)/norm; 293 x2 = theBuffer.GetY(it - 1) / norm; << 296 x2 = theBuffer.GetY(it-1)/norm; 294 y1 = theBuffer.GetX(it); 297 y1 = theBuffer.GetX(it); 295 y2 = theBuffer.GetX(it - 1); << 298 y2 = theBuffer.GetX(it-1); 296 // G4cout << "G4ParticleHPLegendreStore::Sa << 299 // G4cout << "G4ParticleHPLegendreStore::Sample x y "<<x1<<" "<<y1<<" "<<x2<<" "<<y2<<G4endl; 297 return theInt.Interpolate(theManager.GetSche 300 return theInt.Interpolate(theManager.GetScheme(high), rand, x1, x2, y1, y2); 298 } 301 } 299 302 300 G4double << 303 G4double G4ParticleHPLegendreStore::Integrate(G4int k, G4double costh) // still in interpolation; not used anymore 301 G4ParticleHPLegendreStore::Integrate(G4int k, << 302 G4double << 303 { 304 { 304 G4double result = 0.; << 305 G4double result=0; 305 G4ParticleHPFastLegendre theLeg; 306 G4ParticleHPFastLegendre theLeg; 306 // G4cout <<"the COEFFS "<<k<<" "; << 307 // G4cout <<"the COEFFS "<<k<<" "; 307 // G4cout <<theCoeff[k].GetNumberOfPoly()<< << 308 // G4cout <<theCoeff[k].GetNumberOfPoly()<<" "; 308 for (G4int l = 0; l < theCoeff[k].GetNumberO << 309 for(G4int l=0; l<theCoeff[k].GetNumberOfPoly() ; l++) 309 result += theCoeff[k].GetCoeff(l) * theLeg << 310 { 310 // G4cout << theCoeff[k].GetCoeff(l)<<" << 311 result += theCoeff[k].GetCoeff(l)*theLeg.Integrate(l, costh); 311 } << 312 // G4cout << theCoeff[k].GetCoeff(l)<<" "; 312 // G4cout <<G4endl; << 313 } >> 314 // G4cout <<G4endl; 313 return result; 315 return result; 314 } 316 } 315 317