Changeset 1315 for trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include

Timestamp:

Jun 18, 2010, 11:42:07 AM (14 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:

• trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include/G4FermiPhaseSpaceDecay.hh (modified) (4 diffs)

trunk/source/processes/hadronic/models/de_excitation/fermi_breakup/include/G4FermiPhaseSpaceDecay.hh

@@ -33 +33 @@
 #include "G4LorentzVector.hh"
 #include "G4ParticleMomentum.hh"
+#include "Randomize.hh"
+#include "G4Pow.hh"
 #include <CLHEP/Random/RandGamma.h>
 
@@ -41 +43 @@
 {
 public:
-  inline G4FermiPhaseSpaceDecay();
-  inline ~G4FermiPhaseSpaceDecay();
+
+  G4FermiPhaseSpaceDecay();
+  ~G4FermiPhaseSpaceDecay();
   
   inline std::vector<G4LorentzVector*> * 
@@ -48 +51 @@
 
 private:
-  inline G4FermiPhaseSpaceDecay(const G4FermiPhaseSpaceDecay&);
-  inline const G4FermiPhaseSpaceDecay & operator=(const G4FermiPhaseSpaceDecay &); 
-  inline G4bool operator==(const G4FermiPhaseSpaceDecay&);
-  inline G4bool operator!=(const G4FermiPhaseSpaceDecay&);
+
+  G4FermiPhaseSpaceDecay(const G4FermiPhaseSpaceDecay&);
+  const G4FermiPhaseSpaceDecay & operator=(const G4FermiPhaseSpaceDecay &); 
+  G4bool operator==(const G4FermiPhaseSpaceDecay&);
+  G4bool operator!=(const G4FermiPhaseSpaceDecay&);
 
   inline G4double PtwoBody(G4double E, G4double P1, G4double P2) const;
@@ -69 +73 @@
 };
 
-#include "G4FermiPhaseSpaceDecay.icc"
+inline G4double 
+G4FermiPhaseSpaceDecay::PtwoBody(G4double E, G4double P1, G4double P2) const
+{
+  G4double res = -1.0;
+  G4double P = (E+P1+P2)*(E+P1-P2)*(E-P1+P2)*(E-P1-P2)/(4.0*E*E);
+  if (P>0.0) { res = std::sqrt(P); }
+  return res;
+}
+
+inline std::vector<G4LorentzVector*> * G4FermiPhaseSpaceDecay::
+Decay(const G4double parent_mass, const std::vector<G4double>& fragment_masses) const
+{
+  return KopylovNBodyDecay(parent_mass,fragment_masses);
+}
+
+inline G4double 
+G4FermiPhaseSpaceDecay::BetaKopylov(const G4int K) const
+{
+  //JMQ 250410 old algorithm has been commented
+  // Notice that alpha > beta always
+  // const G4double beta = 1.5;
+  // G4double alpha = 1.5*(K-1);
+  // G4double Y1 = CLHEP::RandGamma::shoot(alpha,1);
+  // G4double Y2 = CLHEP::RandGamma::shoot(beta,1);
+  
+  // return Y1/(Y1+Y2);
+
+  G4Pow* g4pow = G4Pow::GetInstance();
+  G4int N = 3*K - 5;
+  G4double xN = G4double(N);
+  G4double F;
+  //G4double Fmax = std::pow((3.*K-5.)/(3.*K-4.),(3.*K-5.)/2.)*std::sqrt(1-((3.*K-5.)/(3.*K-4.))); 
+  // VI variant
+  G4double Fmax = std::sqrt(g4pow->powZ(N, xN/(xN + 1))/(xN + 1)); 
+  G4double chi;
+  do 
+    {
+      chi = G4UniformRand();
+      F = std::sqrt(g4pow->powZ(N, chi)*(1-chi));      
+    } while ( Fmax*G4UniformRand() > F);  
+  return chi;
+
+}
 
 #endif