Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 G4FermiBreakUp::G4FermiBreakUp() : thePhaseSpace(nullptr)
42 {
43  g4calc = G4Pow::GetInstance();
44  Coef = ConstCoeff = 0.0;
45 
46  nmax = 16;
47  NormalizedWeights.resize(nmax,0.0);
48  massRes.reserve(4);
49  frag.resize(4, 0);
50  thePool = G4FermiFragmentsPool::Instance();
51  Initialise();
52 }
53 
55 {}
56 
58 {
59  if(thePhaseSpace != nullptr) { return; }
60  // Kappa = V/V_0 it is used in calculation of Coulomb energy
61  // Nuclear radius r0 (is a model parameter)
62  G4double Kappa = 6.0;
63  G4double r0 = 1.3*CLHEP::fermi;
64 
65  Coef = 0.6*(CLHEP::elm_coupling/r0)/g4calc->Z13(1+G4int(Kappa));
66  ConstCoeff = g4calc->powN(r0/hbarc,3)*Kappa/(6.0*pi*pi);
67 
68  thePhaseSpace = thePool->GetFermiPhaseSpaceDecay();
69 }
70 
72 {
73  return thePool->IsApplicable(Z, A, mass);
74 }
75 
76 G4double G4FermiBreakUp::CoulombBarrier(
77  const std::vector<const G4VFermiFragment*>* conf)
78 {
79  // Calculates Coulomb Barrier (MeV) for given channel with K fragments.
80  G4int SumA = 0;
81  G4int SumZ = 0;
82  G4double CoulombEnergy = 0.;
83  size_t nn = conf->size();
84  for (size_t i=0; i<nn; ++i) {
85  G4int z = (*conf)[i]->GetZ();
86  G4int a = (*conf)[i]->GetA();
87  CoulombEnergy += G4double(z*z)/g4calc->Z13(a);
88  SumA += a;
89  SumZ += z;
90  }
91  CoulombEnergy -= SumZ*SumZ/g4calc->Z13(SumA);
92  return -Coef * CoulombEnergy;
93 }
94 
96  G4Fragment* theNucleus)
97 {
98  if(thePool == nullptr) { Initialise(); }
99  // Calculate Momenta of K fragments
100  G4double M = theNucleus->GetMomentum().m();
101  const std::vector<const G4VFermiFragment*>* conf =
102  SelectConfiguration(theNucleus->GetZ_asInt(),
103  theNucleus->GetA_asInt(), M);
104 
105  // should never happen
106  if(!conf) {
107  theResult->push_back(theNucleus);
108  return;
109  }
110 
111  size_t nn = conf->size();
112 
113  // should never happen
114  if(0 == nn) {
115  theResult->push_back(theNucleus);
116  return;
117  }
118 
119  G4LorentzVector fourMomentum = theNucleus->GetMomentum();
120 
121  // one unstable fragment
122  if(1 == nn) {
123  (*conf)[0]->FillFragment(theResult, fourMomentum);
124 
125  // normal case
126  } else {
127  G4ThreeVector boostVector = fourMomentum.boostVector();
128  massRes.clear();
129  for(size_t i=0; i<nn; ++i) {
130  massRes.push_back( (*conf)[i]->GetTotalEnergy() );
131  }
132  std::vector<G4LorentzVector*>* mom = thePhaseSpace->Decay(M, massRes);
133 
134  // size_t nmom = mom->size();
135  // G4cout << "nmom= " << nmom << G4endl;
136 
137  // Go back to the Lab Frame
138  for (size_t j=0; j<nn; ++j) {
139  (*mom)[j]->boost(boostVector);
140  (*conf)[j]->FillFragment(theResult, *((*mom)[j]));
141  delete (*mom)[j];
142  }
143  delete mom;
144  }
145  delete theNucleus;
146 }
147 
148 const std::vector<const G4VFermiFragment*>*
149 G4FermiBreakUp::SelectConfiguration(G4int Z, G4int A, G4double mass)
150 {
151  const std::vector<const G4VFermiFragment*>* res = 0;
152  // new std::vector<const G4VFermiFragment*>;
153  const std::vector<const G4FermiConfiguration*>* conflist =
154  thePool->GetConfigurationList(Z, A, mass);
155  if(!conflist) { return res; }
156  size_t nn = conflist->size();
157  if(0 < nn) {
158  size_t idx = 0;
159  if(1 < nn) {
160  if(nn > nmax) {
161  nmax = nn;
162  NormalizedWeights.resize(nmax,0.0);
163  }
164  G4double prob = 0.0;
165  for(size_t i=0; i<nn; ++i) {
166  prob += DecayProbability(A, mass, (*conflist)[i]);
167  NormalizedWeights[i] = prob;
168  }
169  prob *= G4UniformRand();
170  for(idx=0; idx<nn; ++idx) {
171  if(NormalizedWeights[idx] >= prob) { break; }
172  }
173  }
174  //const std::vector<const G4VFermiFragment*> flist =
175  res = (*conflist)[idx]->GetFragmentList();
176  //size_t nf = flist.size();
177  //for(size_t i=0; i<nf; ++i) { res->push_back(flist[i]); }
178  //G4cout << "FermiBreakUp: " << nn << " conf; selected idx= "
179  // << idx << " Nprod= " << nf << G4endl;
180  }
181  delete conflist;
182  return res;
183 }
184 
185 G4double G4FermiBreakUp::DecayProbability(G4int A, G4double TotalE,
186  const G4FermiConfiguration* conf)
187  // Decay probability for a given channel with K fragments
188 {
189  // A: Atomic Weight
190  // TotalE: Total energy of nucleus
191 
192  G4double KineticEnergy = TotalE;
193  const std::vector<const G4VFermiFragment*>* flist =
194  conf->GetFragmentList();
195 
196  // Number of fragments
197  size_t K = flist->size();
198  if(K > frag.size()) { frag.resize(K, 0); }
199 
200  for (size_t i=0; i<K; ++i) {
201  frag[i] = (*flist)[i];
202  KineticEnergy -=
203  ((*flist)[i]->GetFragmentMass() + (*flist)[i]->GetExcitationEnergy());
204  }
205 
206  // Check that there is enough energy to produce K fragments
207  KineticEnergy -= CoulombBarrier(flist);
208  if (KineticEnergy <= 0.0) { return 0.0; }
209 
210  // Spin factor S_n
211  G4double S_n = 1.0;
212 
213  // mass factors
214  G4double ProdMass = 1.0;
215  G4double SumMass = 0.0;
216 
217  for (size_t i=0; i<K; ++i) {
218  G4double mass = (*flist)[i]->GetFragmentMass();
219  ProdMass *= mass;
220  SumMass += mass;
221  S_n *= (*flist)[i]->GetPolarization();
222  }
223 
224  G4double MassFactor = ProdMass/SumMass;
225  MassFactor *= std::sqrt(MassFactor);
226 
227  G4double Coeff = g4calc->powN(ConstCoeff*A, K-1);
228 
229  //JMQ 111009 Bug fixed: gamma function for odd K was wrong by a factor 2
230  //VI 251014 Use G4Pow
231  G4double Gamma = 1.0;
232  G4double Energ = twopi*KineticEnergy;
233 
234  // integer argument of Gamma function
235  if ((K-1)%2 != 1) {
236  G4int N = 3*(K-1)/2;
237  Gamma = g4calc->factorial(N - 1);
238  Energ = g4calc->powN(Energ, N);
239 
240  // n + 1/2 argument of Gamma function
241  // http://en.wikipedia.org/wiki/Gamma_function
242  } else {
243  G4int n2 = 3*K - 4;
244  G4int n1 = n2/2;
245 
246  static const G4double sqrtpi = std::sqrt(CLHEP::pi);
247  Gamma = sqrtpi*g4calc->factorial(n2)/
248  (g4calc->powN(4.0,n1)*g4calc->factorial(n1));
249  Energ = g4calc->powN(Energ, n1)*std::sqrt(Energ);
250  }
251 
252  // Permutation Factor G_n
253  // search for identical fragments
254  G4double G_n = 1.0;
255  for(size_t i=0; i<K-1; ++i) {
256  if(frag[i]) {
257  G4int nf = 1;
258  for(size_t j=i+1; j<K; ++j) {
259  if(frag[i] == frag[j]) {
260  ++nf;
261  frag[j] = 0;
262  }
263  }
264  if(1 < nf) { G_n *= g4calc->factorial(nf); }
265  }
266  }
267 
268  G4double Weight = Coeff*MassFactor*S_n*Energ/(G_n*Gamma*KineticEnergy);
269 
270  return Weight;
271 }
272 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
Hep3Vector boostVector() const
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:128
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
int G4int
Definition: G4Types.hh:78
static constexpr double twopi
Definition: G4SIunits.hh:76
G4bool IsApplicable(G4int Z, G4int A, G4double mass) const
#define G4UniformRand()
Definition: Randomize.hh:97
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
double A(double temperature)
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
bool G4bool
Definition: G4Types.hh:79
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:307
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4double factorial(G4int Z) const
Definition: G4Pow.hh:264
virtual ~G4FermiBreakUp()
virtual G4bool IsApplicable(G4int Z, G4int A, G4double mass) const final
std::vector< G4LorentzVector * > * Decay(G4double parent_mass, const std::vector< G4double > &fragment_masses) const
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus) final
virtual void Initialise() final
static constexpr double fermi
Definition: SystemOfUnits.h:83
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
tuple z
Definition: test.py:28
static constexpr double elm_coupling
**D E S C R I P T I O N
Definition: HEPEvtcom.cc:77
static constexpr double pi
Definition: G4SIunits.hh:75
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
static constexpr double pi
Definition: SystemOfUnits.h:54
const G4FermiPhaseSpaceDecay * GetFermiPhaseSpaceDecay() const