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 #include "G4AngularDistributionNP.hh"
00027 #include "G4PhysicalConstants.hh"
00028 #include "G4SystemOfUnits.hh"
00029 #include "Randomize.hh"
00030 #include "G4ios.hh"
00031
00032
00033 #include "G4AngularDistributionNPData.hh"
00034 #include "Randomize.hh"
00035
00036
00037 G4double G4AngularDistributionNP::CosTheta(G4double S, G4double m_1, G4double m_2) const
00038 {
00039 G4int verboseLevel=1;
00040
00041 G4double ek= ((S - sqr(m_1) -sqr(m_2) )/(2*m_1) - m_1 )/GeV ;
00042
00043
00044
00045 G4int je1 = 0;
00046 G4int je2 = NENERGY - 1;
00047 do {
00048 G4int midBin = (je1 + je2)/2;
00049 if (ek < elab[midBin])
00050 je2 = midBin;
00051 else
00052 je1 = midBin;
00053 } while (je2 - je1 > 1);
00054
00055
00056 G4double delab = elab[je2] - elab[je1];
00057
00058
00059
00060 G4float sample = G4UniformRand();
00061 G4int ke1 = 0;
00062 G4int ke2 = NANGLE - 1;
00063 G4double dsig = sig[je2][0] - sig[je1][0];
00064 G4double rc = dsig/delab;
00065 G4double b = sig[je1][0] - rc*elab[je1];
00066 G4double sigint1 = rc*ek + b;
00067 G4double sigint2 = 0.;
00068
00069 if (verboseLevel > 1) G4cout << "sample=" << sample << G4endl
00070 << ek << " " << ke1 << " " << ke2 << " "
00071 << sigint1 << " " << sigint2 << G4endl;
00072
00073 do {
00074 G4int midBin = (ke1 + ke2)/2;
00075 dsig = sig[je2][midBin] - sig[je1][midBin];
00076 rc = dsig/delab;
00077 b = sig[je1][midBin] - rc*elab[je1];
00078 G4double sigint = rc*ek + b;
00079 if (sample < sigint) {
00080 ke2 = midBin;
00081 sigint2 = sigint;
00082 }
00083 else {
00084 ke1 = midBin;
00085 sigint1 = sigint;
00086 }
00087 if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " "
00088 << sigint1 << " " << sigint2 << G4endl;
00089 } while (ke2 - ke1 > 1);
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103 dsig = sigint2 - sigint1;
00104 rc = 1./dsig;
00105 b = ke1 - rc*sigint1;
00106 G4double kint = rc*sample + b;
00107 G4double theta = (0.5 + kint)*pi/180.;
00108
00109
00110
00111
00112
00113 if (verboseLevel > 1) {
00114 G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl;
00115 G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl;
00116 }
00117 G4double costh= std::cos(theta);
00118 return costh;
00119 }
00120
00121 G4double G4AngularDistributionNP::Phi() const
00122 {
00123 return twopi * G4UniformRand();
00124 }