| Geant4 Cross Reference |
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 //
27 // FermiBreakUp de-excitation model
28 // by V. Ivanchenko (July 2016)
29 //
30
31 #include "G4FermiBreakUpVI.hh"
32 #include "G4FermiBreakUpUtil.hh"
33 #include "G4FermiFragmentsPoolVI.hh"
34 #include "G4FermiChannels.hh"
35 #include "G4FermiPair.hh"
36 #include "G4PhysicalConstants.hh"
37 #include "G4NuclearLevelData.hh"
38 #include "G4DeexPrecoParameters.hh"
39 #include "G4PhysicsModelCatalog.hh"
40 #include "Randomize.hh"
41 #include "G4RandomDirection.hh"
42
43 G4FermiFragmentsPoolVI* G4FermiBreakUpVI::fPool = nullptr;
44
45 G4FermiBreakUpVI::G4FermiBreakUpVI()
46 {
47 frag.reserve(10);
48 lvect.reserve(10);
49 secID = G4PhysicsModelCatalog::GetModelID("model_G4FermiBreakUpVI");
50 prob.resize(12,0.0);
51 if (nullptr == fPool) {
52 fPool = new G4FermiFragmentsPoolVI();
53 fPool->Initialise();
54 isFirst = true;
55 }
56 }
57
58 G4FermiBreakUpVI::~G4FermiBreakUpVI()
59 {
60 if (isFirst) {
61 delete fPool;
62 fPool = nullptr;
63 }
64 }
65
66 void G4FermiBreakUpVI::Initialise()
67 {
68 G4DeexPrecoParameters* param =
69 G4NuclearLevelData::GetInstance()->GetParameters();
70 fTolerance = param->GetMinExcitation();
71 fElim = param->GetFBUEnergyLimit();
72 fTimeLim = param->GetMaxLifeTime();
73 if (verbose > 1) {
74 G4cout << "### G4FermiBreakUpVI::Initialise(): the pool is initilized="
75 << fPool->IsInitialized() << " fTolerance(eV)=" << fTolerance/CLHEP::eV
76 << " Elim(MeV)=" << fElim/CLHEP::MeV << G4endl;
77 }
78 }
79
80 G4bool G4FermiBreakUpVI::IsApplicable(G4int Z, G4int A, G4double eexc) const
81 {
82 return (Z < maxZ && A < maxA && eexc <= fElim && fPool->HasDecay(Z, A, eexc));
83 }
84
85 void G4FermiBreakUpVI::BreakFragment(G4FragmentVector* theResult,
86 G4Fragment* theNucleus)
87 {
88 if (verbose > 1) {
89 G4cout << "### G4FermiBreakUpVI::BreakFragment start new fragment "
90 << G4endl;
91 G4cout << *theNucleus << G4endl;
92 }
93 if (!fPool->IsInitialized()) { fPool->Initialise(); }
94
95 // initial fragment
96 G4int Z = theNucleus->GetZ_asInt();
97 G4int A = theNucleus->GetA_asInt();
98 G4double excitation = theNucleus->GetExcitationEnergy();
99 if (!IsApplicable(Z, A, excitation)) { return; }
100 G4double mass = theNucleus->GetGroundStateMass() + excitation;
101 G4LorentzVector lv0 = theNucleus->GetMomentum();
102
103 // sample first decay of an initial state
104 // if not possible to decay - exit
105 if (!SampleDecay(Z, A, mass, excitation, lv0)) { return; }
106
107 G4double time = theNucleus->GetCreationTime();
108 delete theNucleus;
109
110 static const G4int imax = 100;
111
112 // loop over vector of Fermi fragments
113 // vector may grow at each iteraction
114 for (std::size_t i=0; i<frag.size(); ++i) {
115 Z = frag[i]->GetZ();
116 A = frag[i]->GetA();
117 excitation = frag[i]->GetExcitationEnergy();
118 lv0 = lvect[i];
119 G4bool unstable = (IsApplicable(Z, A, excitation) && frag[i]->GetLifeTime() < fTimeLim);
120 if (unstable) {
121 mass = frag[i]->GetTotalEnergy();
122 if (verbose > 1) {
123 G4cout << "# FermiFrag " << i << ". Z= " << Z << " A= " << A
124 << " mass= " << mass << " exc= "
125 << frag[i]->GetExcitationEnergy() << G4endl;
126 }
127 unstable = SampleDecay(Z, A, mass, excitation, lv0);
128 }
129 // stable fragment
130 if (!unstable) {
131 if(verbose > 1) { G4cout << " New G4Fragment" << G4endl; }
132 G4Fragment* f = new G4Fragment(A, Z, lv0);
133 f->SetCreationTime(time);
134 f->SetCreatorModelID(secID);
135 theResult->push_back(f);
136 }
137 // limit the loop
138 if (i == imax) { break; }
139 }
140 frag.clear();
141 lvect.clear();
142 }
143
144 G4bool G4FermiBreakUpVI::SampleDecay(const G4int Z, const G4int A, const G4double mass,
145 const G4double exc, G4LorentzVector& lv0)
146 {
147 const G4FermiChannels* chan = fPool->ClosestChannels(Z, A, mass);
148 if (nullptr == chan) { return false; }
149 std::size_t nn = chan->NumberPairs();
150 if (verbose > 1) {
151 G4cout << "G4FermiBreakUpVI::SampleDecay " << nn << " channels Eex= "
152 << chan->GetExcitation() << G4endl;
153 }
154 if (0 == nn) { return false; }
155 if (nn > prob.size()) { prob.resize(nn, 0.0); }
156
157 const G4FermiPair* fpair = nullptr;
158
159 // one unstable fragment
160 if (1 == nn) {
161 fpair = chan->GetPair(0);
162
163 // more pairs
164 } else {
165
166 G4double q = G4UniformRand();
167 const std::vector<G4FermiPair*>& pvect = chan->GetChannels();
168 std::size_t i{0};
169 G4bool pre = true;
170 if (std::abs(exc - chan->GetExcitation()) < fTolerance) {
171 // static probabilities may be used
172 for (; i<nn; ++i) {
173 if (q <= pvect[i]->Probability()) {
174 fpair = pvect[i];
175 break;
176 }
177 }
178 } else {
179 // recompute probabilities
180 pre = false;
181 G4double ptot = 0.0;
182 for (i=0; i<nn; ++i) {
183 ptot += G4FermiBreakUpUtil::Probability(A, pvect[i]->GetFragment1(),
184 pvect[i]->GetFragment2(),
185 mass, exc);
186 prob[i] = ptot;
187 }
188 ptot *= q;
189 for (i=0; i<nn; ++i) {
190 if(ptot <= prob[i]) {
191 fpair = pvect[i];
192 break;
193 }
194 }
195 }
196 if (verbose > 2) {
197 G4cout << "Probabilities of 2-body decay: Nchannels=" << nn
198 << " channels; i=" << i << " is selected; predefined="
199 << pre << G4endl;
200 for (std::size_t j=0; j<nn; ++j) {
201 G4cout << j << ". ";
202 if (pre) { G4cout << pvect[j]->Probability(); }
203 else { G4cout << prob[j]; }
204 G4cout << " Z1= " << pvect[j]->GetFragment1()->GetZ()
205 << " A1= " << pvect[j]->GetFragment1()->GetA()
206 << " Z2= " << pvect[j]->GetFragment2()->GetZ()
207 << " A2= " << pvect[j]->GetFragment2()->GetA()
208 << G4endl;
209 }
210 }
211 }
212 if (nullptr == fpair) { return false; }
213
214 auto frag1 = fpair->GetFragment1();
215 auto frag2 = fpair->GetFragment2();
216
217 G4double mass1 = frag1->GetTotalEnergy();
218 G4double mass2 = frag2->GetTotalEnergy();
219 if (verbose > 2) {
220 G4cout << " M= " << mass << " M1= " << mass1 << " M2= " << mass2
221 << " Exc1= " << frag1->GetExcitationEnergy()
222 << " Exc2= " << frag2->GetExcitationEnergy() << G4endl;
223 }
224 // sample fragment1
225 G4double e1 = 0.5*(mass*mass - mass2*mass2 + mass1*mass1)/mass;
226 //G4cout << "ekin1= " << e1 - mass1 << G4endl;
227 G4double p1(0.0);
228 if (e1 > mass1) {
229 p1 = std::sqrt((e1 - mass1)*(e1 + mass1));
230 } else {
231 e1 = mass1;
232 }
233 G4LorentzVector lv1 = G4LorentzVector(G4RandomDirection()*p1, e1);
234
235 // compute kinematics
236 auto boostVector = lv0.boostVector();
237 lv1.boost(boostVector);
238 G4LorentzVector lv2 = lv0 - lv1;
239
240 frag.push_back(frag1);
241 frag.push_back(frag2);
242 lvect.push_back(lv1);
243 lvect.push_back(lv2);
244
245 return true;
246 }
247