Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
G4FermiBreakUpVI.cc
이 파일의 문서화 페이지로 가기
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 
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 = 0.0;
56  frag1 = frag2 = nullptr;
57  prob.resize(12,0.0);
58  Initialise();
59 }
60 
62 {
64  delete thePool;
65  thePool = nullptr;
66  }
67 }
68 
70 {
71  if(verbose > 0) {
72  G4cout << "### G4FermiBreakUpVI::Initialise(): " << thePool << G4endl;
73  }
74  if(thePool == nullptr) { InitialisePool(); }
77  //tolerance = thePool->GetTolerance();
78 }
79 
81 {
82 #ifdef G4MULTITHREADED
83  G4MUTEXLOCK(&G4FermiBreakUpVI::FermiBreakUpVIMutex);
84 #endif
85  if(thePool == nullptr) {
87  }
88 #ifdef G4MULTITHREADED
89  G4MUTEXUNLOCK(&G4FermiBreakUpVI::FermiBreakUpVIMutex);
90 #endif
91 }
92 
94 {
95  return (ZZ < maxZ && AA < maxA && AA > 0 && eexc <= elim) ? true : false;
96 }
97 
99  G4Fragment* theNucleus)
100 {
101  if(verbose > 0) {
102  G4cout << "### G4FermiBreakUpVI::BreakFragment start new fragment " << G4endl;
103  G4cout << *theNucleus << G4endl;
104  }
105  frag.clear();
106  lvect.clear();
107 
108  // initial fragment
109  Z = theNucleus->GetZ_asInt();
110  A = theNucleus->GetA_asInt();
111  excitation = theNucleus->GetExcitationEnergy();
112  mass = theNucleus->GetGroundStateMass() + excitation;
113  spin = -1;
114  G4double time = theNucleus->GetCreationTime();
115 
116  lv0 = theNucleus->GetMomentum();
117  rndmEngine = G4Random::getTheEngine();
118 
119  // sample first decay of an initial state
120  if(!SampleDecay()) {
121  theResult->push_back(theNucleus);
122  return;
123  }
124  delete theNucleus;
125 
126  static const G4int imax = 100;
127 
128  // loop over vector of Fermi fragments
129  // vector may grow at each iteraction
130  for(size_t i=0; i<frag.size(); ++i) {
131  Z = frag[i]->GetZ();
132  A = frag[i]->GetA();
133  spin = frag[i]->GetSpin();
134  mass = frag[i]->GetTotalEnergy();
135  excitation = 0.0;
136  if(thePool->IsPhysical(Z, A)) {
137  excitation = frag[i]->GetExcitationEnergy();
138  }
139  lv0 = lvect[i];
140  if(verbose > 0) {
141  G4cout << "# FermiFrag " << i << ". Z= " << Z << " A= " << A
142  << " mass= " << mass << " exc= " << excitation << G4endl;
143  }
144  // stable fragment
145  if(!SampleDecay()) {
146  if(verbose > 0) { G4cout << " New G4Fragment" << G4endl; }
147  G4Fragment* f = new G4Fragment(A, Z, lv0);
148  f->SetSpin(0.5*spin);
149  f->SetCreationTime(time);
150  theResult->push_back(f);
151  }
152  // limit the loop
153  if(i == imax) {
154  break;
155  }
156  }
157 }
158 
160 {
161  const G4FermiChannels* chan = thePool->ClosestChannels(Z, A, mass);
162  if(!chan) { return false; }
163  size_t nn = chan->GetNumberOfChannels();
164  if(verbose > 0) {
165  G4cout << "== SampleDecay " << nn << " channels Eex= "
166  << chan->GetExcitation() << G4endl;
167  }
168  if(0 == nn) { return false; }
169 
170  const G4FermiPair* fpair = nullptr;
171 
172  // one unstable fragment
173  if(1 == nn) {
174  fpair = chan->GetPair(0);
175 
176  // more pairs
177  } else {
178 
179  // in static probabilities may be used
180  if(std::abs(excitation - chan->GetExcitation()) < tolerance) {
181  fpair = chan->SamplePair(rndmEngine->flat());
182 
183  } else {
184 
185  // recompute probabilities
186  const std::vector<const G4FermiPair*>& pvect = chan->GetChannels();
187  if(nn > 12) { prob.resize(nn, 0.0); }
188  G4double ptot = 0.0;
189  if(verbose > 1) {
190  G4cout << "Start recompute probabilities" << G4endl;
191  }
192  for(size_t i=0; i<nn; ++i) {
193  ptot += theDecay->ComputeProbability(Z, A, -1, mass,
194  pvect[i]->GetFragment1(),
195  pvect[i]->GetFragment2());
196  prob[i] = ptot;
197  if(verbose > 1) {
198  G4cout << i << ". " << prob[i]
199  << " Z1= " << pvect[i]->GetFragment1()->GetZ()
200  << " A1= " << pvect[i]->GetFragment1()->GetA()
201  << " Z2= " << pvect[i]->GetFragment2()->GetZ()
202  << " A2= " << pvect[i]->GetFragment2()->GetA()
203  << G4endl;
204  }
205  }
206  ptot *= rndmEngine->flat();
207  for(size_t i=0; i<nn; ++i) {
208  if(ptot <= prob[i] || i+1 == nn) {
209  fpair = pvect[i];
210  break;
211  }
212  }
213  }
214  }
215  if(!fpair) { return false; }
216 
217  frag1 = fpair->GetFragment1();
218  frag2 = fpair->GetFragment2();
219 
220  G4double mass1 = frag1->GetTotalEnergy();
221  G4double mass2 = frag2->GetTotalEnergy();
222  if(verbose > 1) {
223  G4cout << " M= " << mass << " M1= " << mass1 << " M2= " << mass2
224  << " Exc1= " << frag1->GetExcitationEnergy()
225  << " Exc2= " << frag2->GetExcitationEnergy() << G4endl;
226  }
227  // sample fragment1
228  G4double e1 = 0.5*(mass*mass - mass2*mass2 + mass1*mass1)/mass;
229  //G4cout << "ekin1= " << e1 - mass1 << G4endl;
230  G4double p1(0.0);
231  if(e1 > mass1) {
232  p1 = std::sqrt((e1 - mass1)*(e1 + mass1));
233  } else {
234  e1 = mass1;
235  }
237  G4LorentzVector lv1 = G4LorentzVector(v*p1, e1);
238 
239  // compute kinematics
241  lv1.boost(boostVector);
242 
243  G4double e2 = mass - e1;
244  if(e2 < mass2) {
245  e2 = mass2;
246  p1 = 0.0;
247  }
248  lv0.set(-v*p1, e2);
250 
251  frag.push_back(frag1);
252  frag.push_back(frag2);
253  lvect.push_back(lv1);
254  lvect.push_back(lv0);
255 
256  return true;
257 }
258 
G4double GetCreationTime() const
Definition: G4Fragment.hh:441
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:300
std::vector< G4LorentzVector > lvect
std::vector< G4double > prob
void SetCreationTime(G4double time)
Definition: G4Fragment.hh:446
virtual void Initialise() final
const G4FermiFragment * GetFragment1() const
Definition: G4FermiPair.hh:97
#define G4endl
Definition: G4ios.hh:61
const G4FermiPair * SamplePair(G4double rand) const
G4int GetA_asInt() const
Definition: G4Fragment.hh:259
G4bool IsMasterThread()
Definition: G4Threading.cc:130
G4double ComputeProbability(G4int Z, G4int A, G4int spin, G4double TotalE, const G4FermiFragment *f1, const G4FermiFragment *f2) const
virtual ~G4FermiBreakUpVI()
G4ThreeVector G4RandomDirection()
G4double GetGroundStateMass() const
Definition: G4Fragment.hh:281
virtual void BreakFragment(G4FragmentVector *, G4Fragment *theNucleus) final
std::vector< const G4FermiFragment * > frag
static const G4int maxZ
const std::vector< const G4FermiPair * > & GetChannels() const
const G4FermiDecayProbability * FermiDecayProbability() const
virtual void InitialisePool() final
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:88
const G4FermiDecayProbability * theDecay
G4LorentzVector lv0
virtual double flat()=0
G4bool IsPhysical(G4int Z, G4int A) const
CLHEP::HepRandomEngine * rndmEngine
void SetSpin(G4double value)
Definition: G4Fragment.hh:415
static constexpr double MeV
static const G4int maxA
G4ThreeVector boostVector
void set(double x, double y, double z, double t)
G4double GetTotalEnergy(void) const
G4int GetZ_asInt() const
Definition: G4Fragment.hh:264
#define G4MUTEXUNLOCK(mutex)
Definition: G4Threading.hh:234
const G4FermiPair * GetPair(size_t idx) const
size_t GetNumberOfChannels() const
static const G4int imax
int G4int
Definition: G4Types.hh:78
const G4FermiFragment * frag2
const G4FermiFragment * frag1
#define G4MUTEXLOCK(mutex)
Definition: G4Threading.hh:233
const G4FermiChannels * ClosestChannels(G4int Z, G4int A, G4double mass) const
G4GLOB_DLL std::ostream G4cout
virtual G4bool IsApplicable(G4int ZZ, G4int AA, G4double etot) const final
Hep3Vector boostVector() const
G4double GetExcitationEnergy(void) const
static G4FermiFragmentsPoolVI * thePool
CLHEP::HepLorentzVector G4LorentzVector
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:276
std::vector< G4Fragment * > G4FragmentVector
Definition: G4Fragment.hh:63
G4double GetExcitation() const
HepLorentzVector & boost(double, double, double)
const G4FermiFragment * GetFragment2() const
Definition: G4FermiPair.hh:102
std::mutex G4Mutex
Definition: G4Threading.hh:84