Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4FermiBreakUpVI.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: G4VFermiBreakUpVI.cc,v 1.5 2006-06-29 20:13:13 gunter Exp $
27 //
28 // FermiBreakUp de-excitation model
29 // by V. Ivanchenko (July 2016)
30 //
31 
32 #include "G4FermiBreakUpVI.hh"
35 #include "G4FermiChannels.hh"
36 #include "G4FermiPair.hh"
37 #include "G4RandomDirection.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "Randomize.hh"
40 
41 G4FermiFragmentsPoolVI* G4FermiBreakUpVI::thePool = nullptr;
42 
43 #ifdef G4MULTITHREADED
44 G4Mutex G4FermiBreakUpVI::FermiBreakUpVIMutex = G4MUTEX_INITIALIZER;
45 #endif
46 
48  : theDecay(nullptr), rndmEngine(nullptr), verbose(0), maxZ(9), maxA(17)
49 {
50  prob.reserve(10);
51  frag.reserve(10);
52  lvect.reserve(10);
53  Z = A = spin = 0;
54  mass = elim = excitation = tolerance = 0.0;
55  frag1 = frag2 = nullptr;
56  Initialise();
57 }
58 
60 {
62  delete thePool;
63  thePool = nullptr;
64  }
65 }
66 
68 {
69  if(thePool == nullptr) { InitialisePool(); }
70  theDecay = thePool->FermiDecayProbability();
71  elim = thePool->GetEnergyLimit();
72  tolerance = thePool->GetTolerance();
73 }
74 
75 void G4FermiBreakUpVI::InitialisePool()
76 {
77 #ifdef G4MULTITHREADED
78  G4MUTEXLOCK(&G4FermiBreakUpVI::FermiBreakUpVIMutex);
79 #endif
80  if(thePool == nullptr) {
81  thePool = new G4FermiFragmentsPoolVI();
82  }
83 #ifdef G4MULTITHREADED
84  G4MUTEXUNLOCK(&G4FermiBreakUpVI::FermiBreakUpVIMutex);
85 #endif
86 }
87 
89 {
90  return (ZZ < maxZ && AA < maxA && AA > 0 && eexc <= elim) ? true : false;
91 }
92 
94  G4Fragment* theNucleus)
95 {
96  if(verbose > 0) {
97  G4cout << "### G4FermiBreakUpVI::BreakFragment start new fragment " << G4endl;
98  G4cout << *theNucleus << G4endl;
99  }
100  frag.clear();
101  lvect.clear();
102 
103  // initial fragment
104  Z = theNucleus->GetZ_asInt();
105  A = theNucleus->GetA_asInt();
106  excitation = theNucleus->GetExcitationEnergy();
107  mass = theNucleus->GetGroundStateMass() + excitation;
108  spin = -1;
109  G4double time = theNucleus->GetCreationTime();
110 
111  lv0 = theNucleus->GetMomentum();
112  rndmEngine = G4Random::getTheEngine();
113 
114  // sample first decay of an initial state
115  if(!SampleDecay()) {
116  theResult->push_back(theNucleus);
117  return;
118  }
119  delete theNucleus;
120 
121  static const G4int imax = 100;
122 
123  // loop over vector of Fermi fragments
124  // vector may grouw at each iteraction
125  for(size_t i=0; i<frag.size(); ++i) {
126  Z = frag[i]->GetZ();
127  A = frag[i]->GetA();
128  spin = frag[i]->GetSpin();
129  mass = frag[i]->GetTotalEnergy();
130  excitation = frag[i]->GetExcitationEnergy();
131  lv0 = lvect[i];
132  if(verbose > 0) {
133  G4cout << "# FermiFrag " << i << ". Z= " << Z << " A= " << A
134  << " mass= " << mass << " exc= " << excitation << G4endl;
135  }
136  // stable fragment
137  if(!SampleDecay()) {
138  if(verbose > 0) { G4cout << " New G4Fragment" << G4endl; }
139  G4Fragment* f = new G4Fragment(A, Z, lv0);
140  f->SetSpin(0.5*spin);
141  f->SetCreationTime(time);
142  theResult->push_back(f);
143  }
144  // limit the loop
145  if(i == imax) {
146  break;
147  }
148  }
149 }
150 
151 G4bool G4FermiBreakUpVI::SampleDecay()
152 {
153  const G4FermiChannels* chan = thePool->ClosestChannels(Z, A, mass);
154  if(!chan) { return false; }
155  size_t nn = chan->GetNumberOfChannels();
156  if(verbose > 0) {
157  G4cout << "== SampleDecay " << nn << " channels Eex= "
158  << chan->GetExcitation() << G4endl;
159  }
160  if(0 == nn) { return false; }
161 
162  const G4FermiPair* fpair = nullptr;
163 
164  // one unstable fragment
165  if(1 == nn) {
166  fpair = chan->GetPair(0);
167 
168  // more pairs
169  } else {
170 
171  // in static probabilities may be used
172  if(std::abs(excitation - chan->GetExcitation()) < tolerance) {
173  fpair = chan->SamplePair(rndmEngine->flat());
174 
175  } else {
176 
177  // recompute probabilities
178  const std::vector<const G4FermiPair*>& pvect = chan->GetChannels();
179  prob.resize(nn, 0.0);
180  G4double ptot = 0.0;
181  if(verbose > 1) {
182  G4cout << "Start recompute probabilities" << G4endl;
183  }
184  for(size_t i=0; i<nn; ++i) {
185  ptot += theDecay->ComputeProbability(Z, A, -1, mass,
186  pvect[i]->GetFragment1(),
187  pvect[i]->GetFragment2());
188  prob[i] = ptot;
189  if(verbose > 1) {
190  G4cout << i << ". " << prob[i]
191  << " Z1= " << pvect[i]->GetFragment1()->GetZ()
192  << " A1= " << pvect[i]->GetFragment1()->GetA()
193  << " Z2= " << pvect[i]->GetFragment2()->GetZ()
194  << " A2= " << pvect[i]->GetFragment2()->GetA()
195  << G4endl;
196  }
197  }
198  ptot *= rndmEngine->flat();
199  for(size_t i=0; i<nn; ++i) {
200  if(ptot <= prob[i] || i+1 == nn) {
201  fpair = pvect[i];
202  break;
203  }
204  }
205  }
206  }
207  if(!fpair) { return false; }
208 
209  frag1 = fpair->GetFragment1();
210  frag2 = fpair->GetFragment2();
211 
212  G4double mass1 = frag1->GetTotalEnergy();
213  G4double mass2 = frag2->GetTotalEnergy();
214  if(verbose > 1) {
215  G4cout << " M= " << mass << " M1= " << mass1 << " M2= " << mass2
216  << " Exc1= " << frag1->GetExcitationEnergy()
217  << " Exc2= " << frag2->GetExcitationEnergy() << G4endl;
218  }
219  // sample fragment1
220  G4double e1 = 0.5*(mass*mass - mass2*mass2 + mass1*mass1)/mass;
221  //G4cout << "ekin1= " << e1 - mass1 << G4endl;
222  G4double p1(0.0);
223  if(e1 > mass1) {
224  p1 = std::sqrt((e1 - mass1)*(e1 + mass1));
225  } else {
226  e1 = mass1;
227  }
229  G4LorentzVector lv1 = G4LorentzVector(v*p1, e1);
230 
231  // compute kinematics
232  boostVector = lv0.boostVector();
233  lv1.boost(boostVector);
234 
235  G4double e2 = mass - e1;
236  if(e2 < mass2) {
237  e2 = mass2;
238  p1 = 0.0;
239  }
240  lv0.set(-v*p1, e2);
241  lv0.boost(boostVector);
242 
243  frag.push_back(frag1);
244  frag.push_back(frag2);
245  lvect.push_back(lv1);
246  lvect.push_back(lv0);
247 
248  return true;
249 }
250 
#define G4MUTEXUNLOCK
Definition: G4Threading.hh:180
G4double ComputeProbability(G4int Z, G4int A, G4int spin, G4double TotalE, const G4FermiFragment *f1, const G4FermiFragment *f2) const
Hep3Vector boostVector() const
const G4FermiFragment * GetFragment2() const
Definition: G4FermiPair.hh:103
virtual G4bool IsApplicable(G4int ZZ, G4int AA, G4double etot) const final
G4double GetExcitation() const
G4ThreeVector G4RandomDirection()
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus) final
virtual double flat()=0
const G4FermiPair * GetPair(size_t idx) const
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
const G4FermiChannels * ClosestChannels(G4int Z, G4int A, G4double mass) const
const G4FermiFragment * GetFragment1() const
Definition: G4FermiPair.hh:98
const std::vector< const G4FermiPair * > & GetChannels() const
G4GLOB_DLL std::ostream G4cout
G4int GetA_asInt() const
Definition: G4Fragment.hh:266
virtual ~G4FermiBreakUpVI()
G4double GetCreationTime() const
Definition: G4Fragment.hh:448
bool G4bool
Definition: G4Types.hh:79
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:307
HepLorentzVector & boost(double, double, double)
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
const G4FermiDecayProbability * FermiDecayProbability() const
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:288
virtual void Initialise() final
#define G4MUTEXLOCK
Definition: G4Threading.hh:179
tuple v
Definition: test.py:18
G4int G4Mutex
Definition: G4Threading.hh:173
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:453
void set(double x, double y, double z, double t)
G4int GetZ_asInt() const
Definition: G4Fragment.hh:271
const G4FermiPair * SamplePair(G4double rand) const
#define G4endl
Definition: G4ios.hh:61
G4bool IsMasterThread()
Definition: G4Threading.cc:146
void SetSpin(G4double value)
Definition: G4Fragment.hh:422
static const G4int imax
double G4double
Definition: G4Types.hh:76
G4double GetTotalEnergy(void) const
size_t GetNumberOfChannels() const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:283
CLHEP::HepLorentzVector G4LorentzVector
G4double GetExcitationEnergy(void) const