Geant4  10.01.p02
G4FermiPhaseSpaceDecay.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4ExcitationHandler.hh,v 1.13 2010-11-17 16:20:31 vnivanch Exp $
27 //
28 // Hadronic Process: Phase space decay for the Fermi BreakUp model
29 // by V. Lara
30 //
31 // Modifications:
32 // 01.04.2011 General cleanup by V.Ivanchenko:
33 // - IsotropicVector is inlined
34 // - Momentum computation return zero or positive value
35 // - DumpProblem method is added providing more information
36 // - Reduced usage of exotic std functions
37 
38 #include <numeric>
39 
41 #include "G4SystemOfUnits.hh"
42 #include "G4HadronicException.hh"
43 
45 {
47 }
48 
50 {
51 }
52 
53 std::vector<G4LorentzVector*> *
55  const std::vector<G4double>& mr) const
56  // Calculates momentum for N fragments (Kopylov's method of sampling is used)
57 {
58  size_t N = mr.size();
59 
60  std::vector<G4LorentzVector*>* P = new std::vector<G4LorentzVector*>(N, 0);
61 
62  G4double mtot = std::accumulate( mr.begin(), mr.end(), 0.0);
63  G4double mu = mtot;
64  G4double PFragMagCM = 0.0;
65  G4double Mass = M;
66  G4double T = Mass-mtot;
67  G4LorentzVector PFragCM(0.0,0.0,0.0,0.0);
68  G4LorentzVector PRestCM(0.0,0.0,0.0,0.0);
69  G4LorentzVector PRestLab(0.0,0.0,0.0,Mass);
70 
71  for (size_t k = N-1; k>0; --k)
72  {
73  mu -= mr[k];
74  if (k>1) { T *= BetaKopylov(k); }
75  else { T = 0.0; }
76 
77  G4double RestMass = mu + T;
78 
79  PFragMagCM = PtwoBody(Mass,mr[k],RestMass);
80 
81  // Create a unit vector with a random direction isotropically distributed
82  G4ThreeVector RandVector(IsotropicVector(PFragMagCM));
83 
84  PFragCM.setVect(RandVector);
85  PFragCM.setE(std::sqrt(PFragMagCM*PFragMagCM + mr[k]*mr[k]));
86 
87  PRestCM.setVect(-RandVector);
88  PRestCM.setE(std::sqrt(PFragMagCM*PFragMagCM + RestMass*RestMass));
89 
90 
91  G4ThreeVector BoostV = PRestLab.boostVector();
92 
93  PFragCM.boost(BoostV);
94  PRestCM.boost(BoostV);
95  PRestLab = PRestCM;
96 
97  (*P)[k] = new G4LorentzVector(PFragCM);
98 
99  Mass = RestMass;
100  }
101 
102  (*P)[0] = new G4LorentzVector(PRestLab);
103 
104  return P;
105 }
106 
107 
108 std::vector<G4LorentzVector*> *
109 G4FermiPhaseSpaceDecay::NBodyDecay(G4double M, const std::vector<G4double>& mr) const
110 {
111  // Number of fragments
112  size_t N = mr.size();
113  size_t i, j;
114  // Total Daughters Mass
115  G4double mtot = std::accumulate( mr.begin(), mr.end(), 0.0);
116  G4double Emax = M - mtot + mr[0];
117  G4double Emin = 0.0;
118  G4double Wmax = 1.0;
119  for (i = 1; i < N; i++)
120  {
121  Emax += mr[i];
122  Emin += mr[i-1];
123  Wmax *= PtwoBody(Emax, Emin, mr[i]);
124  }
125 
126  G4int ntries = 0;
127  G4double weight = 1.0;
128  std::vector<G4double> p(N, 0.0);
129  std::vector<G4double> r(N,0.0);
130  std::vector<G4double> vm(N, 0.0);
131  r[N-1] = 1.0;
132 
133  do
134  {
135  // Sample uniform random numbers in increasing order
136  for (i = 1; i < N-1; i++) { r[i] = G4UniformRand(); }
137  std::sort(r.begin(),r.end(), std::less<G4double>());
138 
139  // Calculate virtual masses
140  std::partial_sum(mr.begin(), mr.end(), vm.begin());
141  std::transform(r.begin(), r.end(), r.begin(),
142  std::bind2nd(std::multiplies<G4double>(), M-mtot));
143  std::transform(r.begin(), r.end(), vm.begin(), vm.begin(), std::plus<G4double>());
144 
145  // Calcualte daughter momenta
146  weight = 1.0;
147  for (j = 0; j < N-1; j++)
148  {
149  p[j] = PtwoBody(vm[j+1],vm[j],mr[j+1]);
150  weight *= p[j];
151  }
152  p[N-1] = PtwoBody(vm[N-2],mr[N-2],mr[N-1]);
153 
154 
155  if (ntries++ > 1000000)
156  {
157  throw G4HadronicException(__FILE__, __LINE__, "Failed to decay");
158  }
159  }
160  while ( weight < G4UniformRand()*Wmax );
161 
162  std::vector<G4LorentzVector*> * P = new std::vector<G4LorentzVector*>(N, 0);
163 
164  G4ThreeVector a3P = IsotropicVector(p[0]);
165 
166  (*P)[0] = new G4LorentzVector( a3P, std::sqrt(a3P.mag2()+mr[0]*mr[0]) );
167  (*P)[1] = new G4LorentzVector(-a3P, std::sqrt(a3P.mag2()+mr[1]*mr[1]) );
168  for (i = 2; i < N; i++)
169  {
170  a3P = IsotropicVector(p[i-1]);
171  (*P)[i] = new G4LorentzVector(a3P, std::sqrt(a3P.mag2() + mr[i]*mr[i]));
172  G4ThreeVector Beta = -((*P)[i]->boostVector());
173  // boost already created particles
174  for (j = 0; j < i; j++)
175  {
176  (*P)[j]->boost(Beta);
177  }
178  }
179 
180  return P;
181 }
182 
183 std::vector<G4LorentzVector*> *
185  const std::vector<G4double>& mass) const
186 {
187  G4double m0 = mass.front();
188  G4double m1 = mass.back();
189  G4double mom = PtwoBody(M,m0,m1);
191 
192  G4LorentzVector * P41 = new G4LorentzVector;
193  P41->setVect(p);
194  P41->setE(std::sqrt(mom*mom + m0*m0));
195 
196  G4LorentzVector * P42 = new G4LorentzVector;
197  P42->setVect(-p);
198  P42->setE(std::sqrt(mom*mom + m1*m1));
199 
200  std::vector<G4LorentzVector*> * result = new std::vector<G4LorentzVector*>;
201  result->push_back(P41);
202  result->push_back(P42);
203  return result;
204 }
205 
206 void
208  G4double P) const
209 {
210  G4cout << "G4FermiPhaseSpaceDecay: problem of decay of M(GeV)= " << E/GeV
211  << " on M1(GeV)= " << P1/GeV << " and M2(GeV)= " << P2/GeV
212  << " P(MeV)= " << P/MeV << " < 0" << G4endl;
213  // exception only if the problem is numerically significant
214  if(P < -CLHEP::eV) {
215  throw G4HadronicException(__FILE__, __LINE__,"Error in decay kinematics");
216  }
217 }
218 
219 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static const double MeV
Definition: G4SIunits.hh:193
std::vector< G4LorentzVector * > * NBodyDecay(G4double, const std::vector< G4double > &) const
static const G4double * P1[nN]
CLHEP::Hep3Vector G4ThreeVector
std::vector< G4LorentzVector * > * TwoBodyDecay(G4double, const std::vector< G4double > &) const
int G4int
Definition: G4Types.hh:78
G4double PtwoBody(G4double E, G4double P1, G4double P2) const
void DumpProblem(G4double E, G4double P1, G4double P2, G4double P) const
#define G4UniformRand()
Definition: Randomize.hh:93
G4GLOB_DLL std::ostream G4cout
static const double GeV
Definition: G4SIunits.hh:196
static const double eV
Definition: G4SIunits.hh:194
G4ThreeVector IsotropicVector(const G4double Magnitude=1.0) const
static const G4double Emin
static const G4double * P2[nN]
static const G4double Emax
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
std::vector< G4LorentzVector * > * KopylovNBodyDecay(G4double, const std::vector< G4double > &) const
CLHEP::HepLorentzVector G4LorentzVector