| 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 // by V. Lara
27 //
28 // Modified:
29 // 01.04.2008 J.M.Quesada Several changes. Soft cut-off switched off.
30 // 01.05.2008 J.M.Quesada Protection against non-physical preeq.
31 // transitional regime has been set
32 // 03.09.2008 J.M.Quesada for external choice of inverse cross section option
33 // 06.09.2008 J.M.Quesada Also external choices have been added for:
34 // - superimposed Coulomb barrier (useSICB=true)
35 // - "never go back" hipothesis (useNGB=true)
36 // - soft cutoff from preeq. to equlibrium (useSCO=true)
37 // - CEM transition probabilities (useCEMtr=true)
38 // 20.08.2010 V.Ivanchenko Cleanup of the code:
39 // - integer Z and A;
40 // - emission and transition classes created at
41 // initialisation
42 // - options are set at initialisation
43 // - do not use copy-constructors for G4Fragment
44 // 03.01.2012 V.Ivanchenko Added pointer to G4ExcitationHandler to the
45 // constructor
46
47 #include "G4PreCompoundModel.hh"
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4PreCompoundEmission.hh"
51 #include "G4PreCompoundTransitions.hh"
52 #include "G4GNASHTransitions.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4Proton.hh"
55 #include "G4Neutron.hh"
56
57 #include "G4NucleiProperties.hh"
58 #include "G4NuclearLevelData.hh"
59 #include "G4DeexPrecoParameters.hh"
60 #include "Randomize.hh"
61 #include "G4DynamicParticle.hh"
62 #include "G4ParticleTypes.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4LorentzVector.hh"
65 #include "G4Exp.hh"
66 #include "G4PhysicsModelCatalog.hh"
67
68 ////////////////////////////////////////////////////////////////////////////////
69
70 G4PreCompoundModel::G4PreCompoundModel(G4ExcitationHandler* ptr)
71 : G4VPreCompoundModel(ptr,"PRECO")
72 {
73 //G4cout << "### NEW PrecompoundModel " << this << G4endl;
74 if(nullptr == ptr) { SetExcitationHandler(new G4ExcitationHandler()); }
75
76 fNuclData = G4NuclearLevelData::GetInstance();
77 proton = G4Proton::Proton();
78 neutron = G4Neutron::Neutron();
79 modelID = G4PhysicsModelCatalog::GetModelID("model_PRECO");
80 }
81
82 ////////////////////////////////////////////////////////////////////////////////
83
84 G4PreCompoundModel::~G4PreCompoundModel()
85 {
86 delete theEmission;
87 delete theTransition;
88 delete GetExcitationHandler();
89 theResult.Clear();
90 }
91
92 ////////////////////////////////////////////////////////////////////////////////
93
94 void G4PreCompoundModel::BuildPhysicsTable(const G4ParticleDefinition&)
95 {
96 InitialiseModel();
97 }
98
99 ////////////////////////////////////////////////////////////////////////////////
100
101 void G4PreCompoundModel::InitialiseModel()
102 {
103 if(isInitialised) { return; }
104 isInitialised = true;
105
106 //G4cout << "G4PreCompoundModel::InitialiseModel() started" << G4endl;
107
108 G4DeexPrecoParameters* param = fNuclData->GetParameters();
109
110 fLowLimitExc = param->GetPrecoLowEnergy();
111 fHighLimitExc = param->GetPrecoHighEnergy();
112
113 useSCO = param->UseSoftCutoff();
114
115 minZ = param->GetMinZForPreco();
116 minA = param->GetMinAForPreco();
117
118 theEmission = new G4PreCompoundEmission();
119 if(param->UseHETC()) { theEmission->SetHETCModel(); }
120 theEmission->SetOPTxs(param->GetPrecoModelType());
121
122 if(param->UseGNASH()) { theTransition = new G4GNASHTransitions; }
123 else { theTransition = new G4PreCompoundTransitions(); }
124 theTransition->UseNGB(param->NeverGoBack());
125 theTransition->UseCEMtr(param->UseCEM());
126
127 if(param->PrecoDummy()) { isActive = false; }
128
129 GetExcitationHandler()->Initialise();
130 }
131
132 ////////////////////////////////////////////////////////////////////////////////
133
134 G4HadFinalState*
135 G4PreCompoundModel::ApplyYourself(const G4HadProjectile & thePrimary,
136 G4Nucleus & theNucleus)
137 {
138 const G4ParticleDefinition* primary = thePrimary.GetDefinition();
139 if(primary != neutron && primary != proton) {
140 G4ExceptionDescription ed;
141 ed << "G4PreCompoundModel is used for ";
142 if(primary) { ed << primary->GetParticleName(); }
143 G4Exception("G4PreCompoundModel::ApplyYourself()","had0033",FatalException,
144 ed,"");
145 return nullptr;
146 }
147
148 G4int Zp = 0;
149 G4int Ap = 1;
150 if(primary == proton) { Zp = 1; }
151
152 G4double timePrimary=thePrimary.GetGlobalTime();
153
154 G4int A = theNucleus.GetA_asInt();
155 G4int Z = theNucleus.GetZ_asInt();
156
157 //G4cout << "### G4PreCompoundModel::ApplyYourself: A= " << A << " Z= " << Z
158 // << " Ap= " << Ap << " Zp= " << Zp << G4endl;
159 // 4-Momentum
160 G4LorentzVector p = thePrimary.Get4Momentum();
161 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
162 p += G4LorentzVector(0.0,0.0,0.0,mass);
163 //G4cout << "Primary 4-mom " << p << " mass= " << mass << G4endl;
164
165 // prepare fragment
166 G4Fragment anInitialState(A + Ap, Z + Zp, p);
167 anInitialState.SetNumberOfExcitedParticle(2, 1);
168 anInitialState.SetNumberOfHoles(1,0);
169 anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
170 anInitialState.SetCreatorModelID(modelID);
171
172 // call excitation handler
173 G4ReactionProductVector* result = DeExcite(anInitialState);
174
175 // fill particle change
176 theResult.Clear();
177 theResult.SetStatusChange(stopAndKill);
178 for(auto const & prod : *result) {
179 G4DynamicParticle * aNewDP = new G4DynamicParticle(prod->GetDefinition(),
180 prod->GetTotalEnergy(),
181 prod->GetMomentum());
182 G4HadSecondary aNew = G4HadSecondary(aNewDP);
183 G4double time = std::max(prod->GetFormationTime(), 0.0);
184 aNew.SetTime(timePrimary + time);
185 aNew.SetCreatorModelID(prod->GetCreatorModelID());
186 delete prod;
187 theResult.AddSecondary(aNew);
188 }
189 delete result;
190
191 //return the filled particle change
192 return &theResult;
193 }
194
195 ////////////////////////////////////////////////////////////////////////////////
196
197 G4ReactionProductVector* G4PreCompoundModel::DeExcite(G4Fragment& aFragment)
198 {
199 if(!isInitialised) { InitialiseModel(); }
200
201 G4ReactionProductVector * Result = new G4ReactionProductVector;
202 G4double U = aFragment.GetExcitationEnergy();
203 G4int Z = aFragment.GetZ_asInt();
204 G4int A = aFragment.GetA_asInt();
205
206 //G4cout << "### G4PreCompoundModel::DeExcite" << G4endl;
207 //G4cout << aFragment << G4endl;
208
209 // Conditions to skip pre-compound and perform equilibrium emission
210 if (!isActive || (Z < minZ && A < minA) ||
211 U < fLowLimitExc*A || U > A*fHighLimitExc ||
212 0 < aFragment.GetNumberOfLambdas()) {
213 PerformEquilibriumEmission(aFragment, Result);
214 return Result;
215 }
216
217 // main loop
218 G4int count = 0;
219 const G4double ldfact = 12.0/CLHEP::pi2;
220 const G4int countmax = 1000;
221 for (;;) {
222 //G4cout << "### PreCompound loop over fragment" << G4endl;
223 //G4cout << aFragment << G4endl;
224 U = aFragment.GetExcitationEnergy();
225 Z = aFragment.GetZ_asInt();
226 A = aFragment.GetA_asInt();
227 G4int eqExcitonNumber =
228 G4lrint(std::sqrt(ldfact*U*fNuclData->GetLevelDensity(Z, A, U)));
229 //
230 // G4cout<<"Neq="<<EquilibriumExcitonNumber<<G4endl;
231 //
232 // J. M. Quesada (Jan. 08) equilibrium hole number could be used as preeq.
233 // evap. delimiter (IAEA report)
234
235 // Loop for transitions, it is performed while there are
236 // preequilibrium transitions.
237 G4bool isTransition = false;
238
239 // G4cout<<"----------------------------------------"<<G4endl;
240 // G4double NP=aFragment.GetNumberOfParticles();
241 // G4double NH=aFragment.GetNumberOfHoles();
242 // G4double NE=aFragment.GetNumberOfExcitons();
243 // G4cout<<" Ex. Energy="<<aFragment.GetExcitationEnergy()<<G4endl;
244 // G4cout<<"N. excitons="<<NE<<" N. Part="<<NP<<"N. Holes ="<<NH<<G4endl;
245 do {
246 ++count;
247 //G4cout<<"transition number .."<<count
248 // <<" n ="<<aFragment.GetNumberOfExcitons()<<G4endl;
249 // soft cutoff criterium as an "ad-hoc" solution to force
250 // increase in evaporation
251 G4int ne = aFragment.GetNumberOfExcitons();
252 G4bool go_ahead = (ne <= eqExcitonNumber);
253
254 //J. M. Quesada (Apr. 08): soft-cutoff switched off by default
255 if (useSCO && go_ahead) {
256 G4double x = (G4double)(ne - eqExcitonNumber)/(G4double)eqExcitonNumber;
257 if( G4UniformRand() < 1.0 - G4Exp(-x*x/0.32) ) { go_ahead = false; }
258 }
259
260 // JMQ: WARNING: CalculateProbability MUST be called prior to Get!!
261 // (O values would be returned otherwise)
262 G4double transProbability =
263 theTransition->CalculateProbability(aFragment);
264 G4double P1 = theTransition->GetTransitionProb1();
265 G4double P2 = theTransition->GetTransitionProb2();
266 G4double P3 = theTransition->GetTransitionProb3();
267 //G4cout<<"#0 P1="<<P1<<" P2="<<P2<<" P3="<<P3<<G4endl;
268
269 //J.M. Quesada (May 2008) Physical criterium (lamdas) PREVAILS over
270 // approximation (critical exciton number)
271 //V.Ivanchenko (May 2011) added check on number of nucleons
272 // to send a fragment to FermiBreakUp
273 // or check on limits of excitation
274 if(!go_ahead || P1 <= P2+P3 || Z < minZ || A < minA ||
275 U <= fLowLimitExc*A || U > A*fHighLimitExc ||
276 aFragment.GetNumberOfExcitons() <= 0) {
277 // G4cout<<"#4 EquilibriumEmission"<<G4endl;
278 PerformEquilibriumEmission(aFragment,Result);
279 return Result;
280 }
281 G4double emissionProbability =
282 theEmission->GetTotalProbability(aFragment);
283 //G4cout<<"#1 TotalEmissionProbability="<<emissionProbability
284 //<<" Nex= " <<aFragment.GetNumberOfExcitons()<<G4endl;
285 //J.M.Quesada (May 08) this has already been done in order to decide
286 // what to do (preeq-eq)
287 // Sum of all probabilities
288 G4double TotalProbability = emissionProbability + transProbability;
289
290 // Select subprocess
291 if (TotalProbability*G4UniformRand() > emissionProbability) {
292 //G4cout<<"#2 Transition"<<G4endl;
293 // It will be transition to state with a new number of excitons
294 isTransition = true;
295 // Perform the transition
296 theTransition->PerformTransition(aFragment);
297 } else {
298 //G4cout<<"#3 Emission"<<G4endl;
299 // It will be fragment emission
300 isTransition = false;
301 Result->push_back(theEmission->PerformEmission(aFragment));
302 }
303 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
304 } while (isTransition); // end of do loop
305
306 // stop if too many iterations
307 if(count >= countmax) {
308 G4ExceptionDescription ed;
309 ed << "G4PreCompoundModel loop over " << countmax << " iterations; "
310 << "current G4Fragment: \n" << aFragment;
311 G4Exception("G4PreCompoundModel::DeExcite()","had0034",JustWarning,
312 ed,"");
313 PerformEquilibriumEmission(aFragment, Result);
314 return Result;
315 }
316 } // end of for (;;) loop
317 return Result;
318 }
319
320 ////////////////////////////////////////////////////////////////////////////////
321 // Documentation
322 ////////////////////////////////////////////////////////////////////////////////
323
324 void G4PreCompoundModel::ModelDescription(std::ostream& outFile) const
325 {
326 outFile
327 << "The GEANT4 precompound model is considered as an extension of the\n"
328 << "hadron kinetic model. It gives a possibility to extend the low energy range\n"
329 << "of the hadron kinetic model for nucleon-nucleus inelastic collision and it \n"
330 << "provides a ”smooth” transition from kinetic stage of reaction described by the\n"
331 << "hadron kinetic model to the equilibrium stage of reaction described by the\n"
332 << "equilibrium deexcitation models.\n"
333 << "The initial information for calculation of pre-compound nuclear stage\n"
334 << "consists of the atomic mass number A, charge Z of residual nucleus, its\n"
335 << "four momentum P0 , excitation energy U and number of excitons n, which equals\n"
336 << "the sum of the number of particles p (from them p_Z are charged) and the number of\n"
337 << "holes h.\n"
338 << "At the preequilibrium stage of reaction, we follow the exciton model approach in ref. [1],\n"
339 << "taking into account the competition among all possible nuclear transitions\n"
340 << "with ∆n = +2, −2, 0 (which are defined by their associated transition probabilities) and\n"
341 << "the emission of neutrons, protons, deuterons, thritium and helium nuclei (also defined by\n"
342 << "their associated emission probabilities according to exciton model)\n"
343 << "\n"
344 << "[1] K.K. Gudima, S.G. Mashnik, V.D. Toneev, Nucl. Phys. A401 329 (1983)\n"
345 << "\n";
346 }
347
348 ////////////////////////////////////////////////////////////////////////////////
349
350 void G4PreCompoundModel::DeExciteModelDescription(std::ostream& outFile) const
351 {
352 outFile << "description of precompound model as used with DeExcite()" << "\n";
353 }
354
355 ////////////////////////////////////////////////////////////////////////////////
356