Geant4  10.01.p03
G4FermiBreakUp.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: G4VFermiBreakUp.cc,v 1.5 2006-06-29 20:13:13 gunter Exp $
27 //
28 // Hadronic Process: Nuclear De-excitations
29 // by V. Lara (Nov 1998)
30 //
31 // Modifications:
32 // 01.04.2011 General cleanup by V.Ivanchenko - use only one object
33 // theConfigurationList and do not instantiate it at each decay
34 
35 #include "G4FermiBreakUp.hh"
36 #include "G4FermiFragmentsPool.hh"
37 #include "G4PhysicalConstants.hh"
38 #include "Randomize.hh"
39 #include "G4Pow.hh"
40 
41 // Kappa = V/V_0 it is used in calculation of Coulomb energy
42 const G4double Kappa = 6.0;
43 // Nuclear radius r0 (is a model parameter)
44 const G4double r0 = 1.3*CLHEP::fermi;
45 const G4double sqrtpi = std::sqrt(CLHEP::pi);
46 
48 {
51  Coef = 0.6*(CLHEP::elm_coupling/r0)/g4pow->Z13(1+G4int(Kappa));
52  ConstCoeff = g4pow->powN(r0/hbarc,3)*Kappa/(6.0*pi*pi);
53 
55 
56  nmax = 16;
57  NormalizedWeights.resize(nmax,0.0);
58  massRes.reserve(4);
59  frag.resize(4, 0);
60 }
61 
63 {}
64 
66  const std::vector<const G4VFermiFragment*>* conf)
67 {
68  // Calculates Coulomb Barrier (MeV) for given channel with K fragments.
69  G4int SumA = 0;
70  G4int SumZ = 0;
71  G4double CoulombEnergy = 0.;
72  size_t nn = conf->size();
73  for (size_t i=0; i<nn; ++i) {
74  G4int z = (*conf)[i]->GetZ();
75  G4int a = (*conf)[i]->GetA();
76  CoulombEnergy += G4double(z*z)/g4pow->Z13(a);
77  SumA += a;
78  SumZ += z;
79  }
80  CoulombEnergy -= SumZ*SumZ/g4pow->Z13(SumA);
81  return -Coef * CoulombEnergy;
82 }
83 
85 {
86  G4FragmentVector* result = new G4FragmentVector();
87  G4Fragment* nucleus = new G4Fragment(theNucleus);
88  BreakFragment(result, nucleus);
89  return result;
90 }
91 
93  G4Fragment* theNucleus)
94 {
95  // Calculate Momenta of K fragments
96  G4double M = theNucleus->GetMomentum().m();
97  const std::vector<const G4VFermiFragment*>* conf =
98  SelectConfiguration(theNucleus->GetZ_asInt(),
99  theNucleus->GetA_asInt(), M);
100 
101  // should never happen
102  if(!conf) {
103  theResult->push_back(theNucleus);
104  return;
105  }
106 
107  size_t nn = conf->size();
108 
109  // should never happen
110  if(0 == nn) {
111  theResult->push_back(theNucleus);
112  return;
113  }
114 
115  G4LorentzVector fourMomentum = theNucleus->GetMomentum();
116 
117  // one unstable fragment
118  if(1 == nn) {
119  (*conf)[0]->FillFragment(theResult, fourMomentum);
120 
121  // normal case
122  } else {
123  G4ThreeVector boostVector = fourMomentum.boostVector();
124  massRes.clear();
125  for(size_t i=0; i<nn; ++i) {
126  massRes.push_back( (*conf)[i]->GetTotalEnergy() );
127  }
128  std::vector<G4LorentzVector*>* mom = thePhaseSpace->Decay(M, massRes);
129 
130  // size_t nmom = mom->size();
131  // G4cout << "nmom= " << nmom << G4endl;
132 
133  // Go back to the Lab Frame
134  for (size_t j=0; j<nn; ++j) {
135  (*mom)[j]->boost(boostVector);
136  (*conf)[j]->FillFragment(theResult, *((*mom)[j]));
137  delete (*mom)[j];
138  }
139  delete mom;
140  }
141  delete theNucleus;
142 }
143 
144 const std::vector<const G4VFermiFragment*>*
146 {
147  const std::vector<const G4VFermiFragment*>* res = 0;
148  // new std::vector<const G4VFermiFragment*>;
149  const std::vector<const G4FermiConfiguration*>* conflist =
150  thePool->GetConfigurationList(Z, A, mass);
151  if(!conflist) { return res; }
152  size_t nn = conflist->size();
153  if(0 < nn) {
154  size_t idx = 0;
155  if(1 < nn) {
156  if(nn > nmax) {
157  nmax = nn;
158  NormalizedWeights.resize(nmax,0.0);
159  }
160  G4double prob = 0.0;
161  for(size_t i=0; i<nn; ++i) {
162  prob += DecayProbability(A, mass, (*conflist)[i]);
163  NormalizedWeights[i] = prob;
164  }
165  prob *= G4UniformRand();
166  for(idx=0; idx<nn; ++idx) {
167  if(NormalizedWeights[idx] >= prob) { break; }
168  }
169  }
170  //const std::vector<const G4VFermiFragment*> flist =
171  res = (*conflist)[idx]->GetFragmentList();
172  //size_t nf = flist.size();
173  //for(size_t i=0; i<nf; ++i) { res->push_back(flist[i]); }
174  //G4cout << "FermiBreakUp: " << nn << " conf; selected idx= "
175  // << idx << " Nprod= " << nf << G4endl;
176  }
177  delete conflist;
178  return res;
179 }
180 
182  const G4FermiConfiguration* conf)
183  // Decay probability for a given channel with K fragments
184 {
185  // A: Atomic Weight
186  // TotalE: Total energy of nucleus
187 
188  G4double KineticEnergy = TotalE;
189  const std::vector<const G4VFermiFragment*>* flist =
190  conf->GetFragmentList();
191 
192  // Number of fragments
193  size_t K = flist->size();
194  if(K > frag.size()) { frag.resize(K, 0); }
195 
196  for (size_t i=0; i<K; ++i) {
197  frag[i] = (*flist)[i];
198  KineticEnergy -=
199  ((*flist)[i]->GetFragmentMass() + (*flist)[i]->GetExcitationEnergy());
200  }
201 
202  // Check that there is enough energy to produce K fragments
203  KineticEnergy -= CoulombBarrier(flist);
204  if (KineticEnergy <= 0.0) { return 0.0; }
205 
206  // Spin factor S_n
207  G4double S_n = 1.0;
208 
209  // mass factors
210  G4double ProdMass = 1.0;
211  G4double SumMass = 0.0;
212 
213  for (size_t i=0; i<K; ++i) {
214  G4double mass = (*flist)[i]->GetFragmentMass();
215  ProdMass *= mass;
216  SumMass += mass;
217  S_n *= (*flist)[i]->GetPolarization();
218  }
219 
220  G4double MassFactor = ProdMass/SumMass;
221  MassFactor *= std::sqrt(MassFactor);
222 
223  G4double Coeff = g4pow->powN(ConstCoeff*A, K-1);
224 
225  //JMQ 111009 Bug fixed: gamma function for odd K was wrong by a factor 2
226  //VI 251014 Use G4Pow
227  G4double Gamma = 1.0;
228  G4double Energ = twopi*KineticEnergy;
229 
230  // integer argument of Gamma function
231  if ((K-1)%2 != 1) {
232  G4int N = 3*(K-1)/2;
233  Gamma = g4pow->factorial(N - 1);
234  Energ = g4pow->powN(Energ, N);
235 
236  // n + 1/2 argument of Gamma function
237  // http://en.wikipedia.org/wiki/Gamma_function
238  } else {
239  G4int n2 = 3*K - 4;
240  G4int n1 = n2/2;
241  Gamma = sqrtpi*g4pow->factorial(n2)/
242  (g4pow->powN(4.0,n1)*g4pow->factorial(n1));
243  Energ = g4pow->powN(Energ, n1)*std::sqrt(Energ);
244  }
245 
246  // Permutation Factor G_n
247  // search for identical fragments
248  G4double G_n = 1.0;
249  for(size_t i=0; i<K-1; ++i) {
250  if(frag[i]) {
251  G4int nf = 1;
252  for(size_t j=i+1; j<K; ++j) {
253  if(frag[i] == frag[j]) {
254  ++nf;
255  frag[j] = 0;
256  }
257  }
258  if(1 < nf) { G_n *= g4pow->factorial(nf); }
259  }
260  }
261 
262  G4double Weight = Coeff*MassFactor*S_n*Energ/(G_n*Gamma*KineticEnergy);
263 
264  return Weight;
265 }
266 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
CLHEP::Hep3Vector G4ThreeVector
G4double z
Definition: TRTMaterials.hh:39
const G4double pi
G4double a
Definition: TRTMaterials.hh:39
const G4double Kappa
int G4int
Definition: G4Types.hh:78
std::vector< G4LorentzVector * > * Decay(const G4double, const std::vector< G4double > &) const
const G4double sqrtpi
G4double CoulombBarrier(const std::vector< const G4VFermiFragment * > *v)
#define G4UniformRand()
Definition: Randomize.hh:93
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
G4int GetA_asInt() const
Definition: G4Fragment.hh:243
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:276
const G4FermiPhaseSpaceDecay * thePhaseSpace
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:65
G4double factorial(G4int Z) const
Definition: G4Pow.hh:266
std::vector< G4double > NormalizedWeights
virtual ~G4FermiBreakUp()
void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus)
G4FermiFragmentsPool * thePool
static const G4double A[nN]
const std::vector< const G4VFermiFragment * > * SelectConfiguration(G4int Z, G4int A, G4double mass)
std::vector< const G4VFermiFragment * > frag
std::vector< G4double > massRes
G4double ConstCoeff
G4int GetZ_asInt() const
Definition: G4Fragment.hh:248
double G4double
Definition: G4Types.hh:76
static G4FermiFragmentsPool * Instance()
const std::vector< const G4VFermiFragment * > * GetFragmentList() const
const std::vector< const G4FermiConfiguration * > * GetConfigurationList(G4int Z, G4int A, G4double mass) const
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)
static const double fermi
Definition: G4SIunits.hh:93
const G4FermiPhaseSpaceDecay * GetFermiPhaseSpaceDecay() const
const G4double r0
G4double DecayProbability(G4int A, G4double TotalE, const G4FermiConfiguration *)
CLHEP::HepLorentzVector G4LorentzVector