00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef G4CASCADE_SAMPLER_ICC
00027 #define G4CASCADE_SAMPLER_ICC
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "Randomize.hh"
00041 #include <iostream>
00042 #include <vector>
00043
00044
00045 template <int NBINS, int NMULT> inline
00046 G4double G4CascadeSampler<NBINS,NMULT>::
00047 findCrossSection(double ke,
00048 const G4double (&xsec)[energyBins]) const {
00049 return interpolator.interpolate(ke, xsec);
00050 }
00051
00052 template <int NBINS, int NMULT> inline
00053 G4int G4CascadeSampler<NBINS,NMULT>::
00054 findMultiplicity(G4double ke,
00055 const G4double xmult[][energyBins]) const {
00056 fillSigmaBuffer(ke, xmult);
00057 return sampleFlat() + 2;
00058 }
00059
00060 template <int NBINS, int NMULT> inline
00061 G4int G4CascadeSampler<NBINS,NMULT>::
00062 findFinalStateIndex(G4int mult, G4double ke, const G4int index[],
00063 const G4double xsec[][energyBins]) const {
00064 G4int start = index[mult-2];
00065 G4int stop = index[mult-1];
00066 if (stop-start <= 1) return start;
00067
00068 fillSigmaBuffer(ke, xsec, start, stop);
00069 return sampleFlat();
00070 }
00071
00072
00073 template <int NBINS, int NMULT> inline
00074 void G4CascadeSampler<NBINS,NMULT>::
00075 fillSigmaBuffer(G4double ke, const G4double x[][energyBins],
00076 G4int startBin, G4int stopBin) const {
00077 sigmaBuf.clear();
00078 if (stopBin-startBin <= 1) return;
00079
00080
00081 sigmaBuf.reserve(stopBin-startBin);
00082 for(G4int i = startBin; i < stopBin; i++)
00083 sigmaBuf.push_back(interpolator.interpolate(ke, x[i]));
00084 }
00085
00086
00087 template <int NBINS, int NMULT> inline
00088 G4int G4CascadeSampler<NBINS,NMULT>::sampleFlat() const {
00089 G4int nbins = sigmaBuf.size();
00090 if (nbins <= 1) return 0;
00091
00092 #ifdef G4CASCADE_DEBUG_SAMPLER
00093 G4cout << "G4CascadeSampler::sampleFlat() has " << nbins << "bins:" << G4endl;
00094 for (G4int sbi=0; sbi<nbins; sbi++) G4cout << " " << sigmaBuf[sbi];
00095 G4cout << G4endl;
00096 #endif
00097
00098 G4int i;
00099 G4double fsum = 0.;
00100 for (i = 0; i < nbins; i++) fsum += sigmaBuf[i];
00101 #ifdef G4CASCADE_DEBUG_SAMPLER
00102 G4cout << " buffer total (fsum) " << fsum;
00103 #endif
00104 fsum *= G4UniformRand();
00105 #ifdef G4CASCADE_DEBUG_SAMPLER
00106 G4cout << " *random-scale " << fsum << G4endl;
00107 #endif
00108
00109 G4double partialSum = 0.0;
00110 for (i = 0; i < nbins; i++) {
00111 partialSum += sigmaBuf[i];
00112 if (fsum < partialSum) return i;
00113 }
00114
00115 return 0;
00116 }
00117
00118
00119 template <int NBINS, int NMULT> inline
00120 void G4CascadeSampler<NBINS,NMULT>::print(std::ostream& os) const {
00121 interpolator.printBins(os);
00122 }
00123
00124 #endif