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
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #ifndef G4FermiPhaseSpaceDecay_hh
00041 #define G4FermiPhaseSpaceDecay_hh 1
00042
00043 #include <vector>
00044 #include <CLHEP/Units/PhysicalConstants.h>
00045
00046 #include "G4LorentzVector.hh"
00047 #include "G4ThreeVector.hh"
00048 #include "Randomize.hh"
00049 #include "G4Pow.hh"
00050
00051 class G4FermiPhaseSpaceDecay
00052 {
00053 public:
00054
00055 G4FermiPhaseSpaceDecay();
00056 ~G4FermiPhaseSpaceDecay();
00057
00058 inline std::vector<G4LorentzVector*> *
00059 Decay(const G4double, const std::vector<G4double>&) const;
00060
00061 private:
00062
00063 inline G4double PtwoBody(G4double E, G4double P1, G4double P2) const;
00064
00065 inline G4ThreeVector IsotropicVector(const G4double Magnitude = 1.0) const;
00066
00067 inline G4double BetaKopylov(G4int) const;
00068
00069 std::vector<G4LorentzVector*> *
00070 TwoBodyDecay(G4double, const std::vector<G4double>&) const;
00071
00072 std::vector<G4LorentzVector*> *
00073 NBodyDecay(G4double, const std::vector<G4double>&) const;
00074
00075 std::vector<G4LorentzVector*> *
00076 KopylovNBodyDecay(G4double, const std::vector<G4double>&) const;
00077
00078 void DumpProblem(G4double E, G4double P1, G4double P2, G4double P) const;
00079
00080 G4FermiPhaseSpaceDecay(const G4FermiPhaseSpaceDecay&);
00081 const G4FermiPhaseSpaceDecay & operator=(const G4FermiPhaseSpaceDecay &);
00082 G4bool operator==(const G4FermiPhaseSpaceDecay&);
00083 G4bool operator!=(const G4FermiPhaseSpaceDecay&);
00084
00085 G4Pow* g4pow;
00086 };
00087
00088 inline G4double
00089 G4FermiPhaseSpaceDecay::PtwoBody(G4double E, G4double P1, G4double P2) const
00090 {
00091 G4double res = 0.0;
00092 G4double P = (E+P1+P2)*(E+P1-P2)*(E-P1+P2)*(E-P1-P2)/(4.0*E*E);
00093 if (P>0.0) { res = std::sqrt(P); }
00094 else { DumpProblem(E,P1,P2,P); }
00095 return res;
00096 }
00097
00098 inline std::vector<G4LorentzVector*> * G4FermiPhaseSpaceDecay::
00099 Decay(G4double parent_mass, const std::vector<G4double>& fragment_masses) const
00100 {
00101 return KopylovNBodyDecay(parent_mass,fragment_masses);
00102 }
00103
00104 inline G4double G4FermiPhaseSpaceDecay::BetaKopylov(G4int K) const
00105 {
00106 G4int N = 3*K - 5;
00107 G4double xN = G4double(N);
00108 G4double F;
00109
00110
00111 G4double Fmax = std::sqrt(g4pow->powN(xN/(xN + 1),N)/(xN + 1));
00112 G4double chi;
00113 do {
00114 chi = G4UniformRand();
00115 F = std::sqrt(g4pow->powN(chi,N)*(1-chi));
00116 } while ( Fmax*G4UniformRand() > F);
00117 return chi;
00118 }
00119
00120 inline G4ThreeVector
00121 G4FermiPhaseSpaceDecay::IsotropicVector(G4double Magnitude) const
00122
00123
00124 {
00125 G4double CosTheta = 2.0*G4UniformRand() - 1.0;
00126 G4double SinTheta = std::sqrt((1. - CosTheta)*(1. + CosTheta));
00127 G4double Phi = CLHEP::twopi*G4UniformRand();
00128 G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
00129 Magnitude*std::sin(Phi)*SinTheta,
00130 Magnitude*CosTheta);
00131 return Vector;
00132 }
00133
00134 #endif