Geant4_10
G4FermiConfigurationList.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 - more clean usage of static
33 // 23.04.2011 V.Ivanchenko: make this class to be responsible for
34 // selection of decay channel and decay
35 
36 #include <set>
37 
39 #include "G4FermiFragmentsPool.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "Randomize.hh"
42 #include "G4Pow.hh"
43 
44 const G4double G4FermiConfigurationList::Kappa = 6.0;
45 const G4double G4FermiConfigurationList::r0 = 1.3*CLHEP::fermi;
46 
48 {
49  thePool = G4FermiFragmentsPool::Instance();
50  g4pow = G4Pow::GetInstance();
51  Coef = 0.6*(CLHEP::elm_coupling/r0)/g4pow->Z13(1+G4int(Kappa));
52  ConstCoeff = g4pow->powN(r0/hbarc,3)*Kappa*std::sqrt(2.0/pi)/3.0;
53 
54  // 16 is the max number
55  nmax = 50;
56  NormalizedWeights.resize(nmax,0.0);
57 }
58 
60 {}
61 
62 G4double
63 G4FermiConfigurationList::CoulombBarrier(
64  const std::vector<const G4VFermiFragment*>& conf)
65 {
66  // Calculates Coulomb Barrier (MeV) for given channel with K fragments.
67  G4int SumA = 0;
68  G4int SumZ = 0;
69  G4double CoulombEnergy = 0.;
70  size_t nn = conf.size();
71  for (size_t i=0; i<nn; ++i) {
72  G4int z = conf[i]->GetZ();
73  G4int a = conf[i]->GetA();
74  CoulombEnergy += G4double(z*z)/g4pow->Z13(a);
75  SumA += a;
76  SumZ += z;
77  }
78  CoulombEnergy -= SumZ*SumZ/g4pow->Z13(SumA);
79  return -Coef * CoulombEnergy;
80 }
81 
82 G4double
83 G4FermiConfigurationList::DecayProbability(G4int A, G4double TotalE,
85  // Decay probability for a given channel with K fragments
86 {
87  // A: Atomic Weight
88  // TotalE: Total energy of nucleus
89 
90  G4double KineticEnergy = TotalE; // MeV
91  G4double ProdMass = 1.0;
92  G4double SumMass = 0.0;
93  G4double S_n = 1.0;
94  std::set<G4int> combSet;
95  std::multiset<G4int> combmSet;
96 
97  const std::vector<const G4VFermiFragment*> flist =
98  conf->GetFragmentList();
99 
100  // Number of fragments
101  G4int K = flist.size();
102 
103  for (G4int i=0; i<K; ++i) {
104  G4int a = flist[i]->GetA();
105  combSet.insert(a);
106  combmSet.insert(a);
107  G4double mass = flist[i]->GetFragmentMass();
108  ProdMass *= mass;
109  SumMass += mass;
110  // Spin factor S_n
111  S_n *= flist[i]->GetPolarization();
112  KineticEnergy -= mass + flist[i]->GetExcitationEnergy();
113  }
114 
115  // Check that there is enough energy to produce K fragments
116  KineticEnergy -= CoulombBarrier(flist);
117  if (KineticEnergy <= 0.0) { return 0.0; }
118 
119  G4double MassFactor = ProdMass/SumMass;
120  MassFactor *= std::sqrt(MassFactor);
121 
122  // This is the constant (doesn't depend on nucleus) part
123  G4double Coeff = g4pow->powN(ConstCoeff*A, K-1);
124 
125  //JMQ 111009 Bug fixed: gamma function for odd K was wrong by a factor 2
126  // new implementation explicitely according to standard properties of Gamma function
127  // Calculation of 1/Gamma(3(k-1)/2)
128  G4double Gamma = 1.0;
129  // G4double arg = 3.0*(K-1)/2.0 - 1.0;
130  // while (arg > 1.1)
131  // {
132  // Gamma *= arg;
133  // arg--;
134  // }
135  // if ((K-1)%2 == 1) Gamma *= std::sqrt(pi);
136 
137  if ((K-1)%2 != 1)
138 
139  {
140  G4double arg = 3.0*(K-1)/2.0 - 1.0;
141  while (arg > 1.1)
142  {
143  Gamma *= arg;
144  arg--;
145  }
146  }
147  else {
148  G4double n= 3.0*K/2.0-2.0;
149  G4double arg2=2*n-1;
150  while (arg2>1.1)
151  {
152  Gamma *= arg2;
153  arg2 -= 2;
154  }
155  Gamma *= std::sqrt(pi)/g4pow->powZ(2,n);
156  }
157  // end of new implementation of Gamma function
158 
159 
160  // Permutation Factor G_n
161  G4double G_n = 1.0;
162  for (std::set<G4int>::iterator itr = combSet.begin();
163  itr != combSet.end(); ++itr)
164  {
165  for (G4int ni = combmSet.count(*itr); ni > 1; ni--) { G_n *= ni; }
166  }
167 
168  G4double Weight = Coeff * MassFactor * (S_n / G_n) / Gamma;
169  Weight *= std::sqrt(g4pow->powN(KineticEnergy,3*(K-1)))/KineticEnergy;
170 
171  return Weight;
172 }
173 
176 {
177  // Calculate Momenta of K fragments
178  G4double M = theNucleus.GetMomentum().m();
179  const std::vector<const G4VFermiFragment*>* conf =
180  SelectConfiguration(theNucleus.GetZ_asInt(),
181  theNucleus.GetA_asInt(), M);
182 
183 
184  G4FragmentVector* theResult = new G4FragmentVector();
185  size_t nn = conf->size();
186  if(1 >= nn) {
187  theResult->push_back(new G4Fragment(theNucleus));
188  delete conf;
189  return theResult;
190  }
191 
192  G4ThreeVector boostVector = theNucleus.GetMomentum().boostVector();
193  std::vector<G4double> mr;
194  mr.reserve(nn);
195  for(size_t i=0; i<nn; ++i) {
196  mr.push_back( (*conf)[i]->GetTotalEnergy() );
197  }
198  std::vector<G4LorentzVector*>* mom = thePhaseSpace.Decay(M,mr);
199  if(!mom) {
200  delete conf;
201  return theResult;
202  }
203 
204  size_t nmom = mom->size();
205 
206  // Go back to the Lab Frame
207  if(0 < nmom) {
208  for (size_t j=0; j<nmom; ++j) {
209  G4LorentzVector* FourMomentum = (*mom)[j];
210 
211  // Lorentz boost
212  FourMomentum->boost(boostVector);
213 
214  G4FragmentVector* fragment = (*conf)[j]->GetFragment(*FourMomentum);
215 
216  size_t nfrag = fragment->size();
217  for (size_t k=0; k<nfrag; ++k) { theResult->push_back((*fragment)[k]); }
218  delete fragment;
219  delete (*mom)[j];
220  }
221  }
222 
223  delete mom;
224  delete conf;
225  return theResult;
226 }
227 
228 const std::vector<const G4VFermiFragment*>*
229 G4FermiConfigurationList::SelectConfiguration(G4int Z, G4int A, G4double mass)
230 {
231  std::vector<const G4VFermiFragment*>* res =
232  new std::vector<const G4VFermiFragment*>;
233  const std::vector<G4FermiConfiguration*>* conflist =
234  thePool->GetConfigurationList(Z, A, mass);
235  if(!conflist) { return res; }
236  size_t nn = conflist->size();
237  if(0 < nn) {
238  size_t idx = 0;
239  if(1 < nn) {
240  if(nn > nmax) {
241  nmax = nn;
242  NormalizedWeights.resize(nmax,0.0);
243  }
244  G4double prob = 0.0;
245  for(size_t i=0; i<nn; ++i) {
246  prob += DecayProbability(A, mass, (*conflist)[i]);
247  NormalizedWeights[i] = prob;
248  }
249  prob *= G4UniformRand();
250  for(idx=0; idx<nn; ++idx) {
251  if(NormalizedWeights[idx] >= prob) { break; }
252  }
253  }
254  const std::vector<const G4VFermiFragment*> flist =
255  (*conflist)[idx]->GetFragmentList();
256  size_t nf = flist.size();
257  for(size_t i=0; i<nf; ++i) { res->push_back(flist[i]); }
258  //G4cout << "FermiBreakUp: " << nn << " conf; selected idx= "
259  // << idx << " Nprod= " << nf << G4endl;
260  }
261  delete conflist;
262  return res;
263 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
Hep3Vector boostVector() const
G4double powN(G4double x, G4int n) const
Definition: G4Pow.cc:125
tuple a
Definition: test.py:11
const std::vector< const G4VFermiFragment * > & GetFragmentList()
G4FragmentVector * GetFragments(const G4Fragment &theNucleus)
int G4int
Definition: G4Types.hh:78
std::vector< G4LorentzVector * > * Decay(const G4double, const std::vector< G4double > &) const
Char_t n[5]
#define G4UniformRand()
Definition: Randomize.hh:87
G4double Z13(G4int Z) const
Definition: G4Pow.hh:129
G4int GetA_asInt() const
Definition: G4Fragment.hh:218
Float_t Z
Definition: plot.C:39
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:251
HepLorentzVector & boost(double, double, double)
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:64
const std::vector< G4FermiConfiguration * > * GetConfigurationList(G4int Z, G4int A, G4double mass)
G4int GetZ_asInt() const
Definition: G4Fragment.hh:223
tuple z
Definition: test.py:28
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:258
double G4double
Definition: G4Types.hh:76
static G4FermiFragmentsPool * Instance()