| 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 // ABLAXX statistical de-excitation model
27 // Jose Luis Rodriguez, UDC (translation from ABLA07 and contact person)
28 // Pekka Kaitaniemi, HIP (initial translation of ablav3p)
29 // Aleksandra Kelic, GSI (ABLA07 code)
30 // Davide Mancusi, CEA (contact person INCL)
31 // Aatos Heikkinen, HIP (project coordination)
32 //
33
34 #include "globals.hh"
35 #include <cmath>
36 #include <iostream>
37
38 #include "G4AblaInterface.hh"
39 #include "G4DoubleHyperDoubleNeutron.hh"
40 #include "G4DoubleHyperH4.hh"
41 #include "G4DynamicParticle.hh"
42 #include "G4ExcitationHandler.hh"
43 #include "G4HyperAlpha.hh"
44 #include "G4HyperH4.hh"
45 #include "G4HyperHe5.hh"
46 #include "G4HyperTriton.hh"
47 #include "G4IonTable.hh"
48 #include "G4ParticleDefinition.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4PhysicsModelCatalog.hh"
51 #include "G4ReactionProduct.hh"
52 #include "G4ReactionProductVector.hh"
53 #include "G4SystemOfUnits.hh"
54
55 G4AblaInterface::G4AblaInterface(G4ExcitationHandler* ptr)
56 : G4VPreCompoundModel(ptr, "ABLAXX")
57 , ablaResult(new G4VarNtp)
58 , theABLAModel(new G4Abla(ablaResult))
59 , eventNumber(0)
60 , secID(-1)
61 , isInitialised(false)
62 {
63 secID = G4PhysicsModelCatalog::GetModelID("model_" + GetModelName());
64 // G4cout << "### NEW PrecompoundModel " << this << G4endl;
65 if (!ptr)
66 SetExcitationHandler(new G4ExcitationHandler);
67 InitialiseModel();
68 G4cout << G4endl << "G4AblaInterface::InitialiseModel() was right." << G4endl;
69 }
70
71 G4AblaInterface::~G4AblaInterface()
72 {
73 applyYourselfResult.Clear();
74 delete ablaResult;
75 delete theABLAModel;
76 delete GetExcitationHandler();
77 }
78
79 void G4AblaInterface::BuildPhysicsTable(const G4ParticleDefinition&) { InitialiseModel(); }
80
81 void G4AblaInterface::InitialiseModel()
82 {
83 if (isInitialised)
84 return;
85 isInitialised = true;
86 theABLAModel->initEvapora();
87 theABLAModel->SetParameters();
88 GetExcitationHandler()->Initialise();
89 }
90
91 G4HadFinalState* G4AblaInterface::ApplyYourself(const G4HadProjectile& thePrimary, G4Nucleus& theNucleus)
92 {
93 // This method is adapted from G4PreCompoundModel::ApplyYourself,
94 // and it is used only by Binary Cascade (BIC) when the latter is coupled with
95 // Abla for nuclear de-excitation. This method allows BIC+ABLA to be used also
96 // for proton and neutron projectile with kinetic energies below 45 MeV, by
97 // creating a "compound" nucleus made by the system "target nucleus +
98 // projectile", before calling the DeExcite method.
99 const G4ParticleDefinition* primary = thePrimary.GetDefinition();
100 if (primary != G4Neutron::Definition() && primary != G4Proton::Definition())
101 {
102 G4ExceptionDescription ed;
103 ed << "G4AblaModel is used for ";
104 if (primary)
105 ed << primary->GetParticleName();
106 G4Exception("G4AblaInterface::ApplyYourself()", "had040", FatalException, ed, "");
107 return nullptr;
108 }
109
110 G4int Zp = 0;
111 G4int Ap = 1;
112 if (primary == G4Proton::Definition())
113 Zp = 1;
114 G4double timePrimary = thePrimary.GetGlobalTime();
115 G4int A = theNucleus.GetA_asInt();
116 G4int Z = theNucleus.GetZ_asInt();
117 G4LorentzVector p = thePrimary.Get4Momentum();
118 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
119 p += G4LorentzVector(0.0, 0.0, 0.0, mass);
120
121 G4Fragment anInitialState(A + Ap, Z + Zp, p);
122 anInitialState.SetNumberOfExcitedParticle(1, Zp);
123 anInitialState.SetNumberOfHoles(1, Zp);
124 anInitialState.SetCreationTime(thePrimary.GetGlobalTime());
125 anInitialState.SetCreatorModelID(secID);
126
127 G4ReactionProductVector* deExciteResult = DeExcite(anInitialState);
128
129 applyYourselfResult.Clear();
130 applyYourselfResult.SetStatusChange(stopAndKill);
131 for (auto const& prod : *deExciteResult)
132 {
133 G4DynamicParticle* aNewDP =
134 new G4DynamicParticle(prod->GetDefinition(), prod->GetTotalEnergy(), prod->GetMomentum());
135 G4HadSecondary aNew = G4HadSecondary(aNewDP);
136 G4double time = std::max(prod->GetFormationTime(), 0.0);
137 aNew.SetTime(timePrimary + time);
138 aNew.SetCreatorModelID(prod->GetCreatorModelID());
139 delete prod;
140 applyYourselfResult.AddSecondary(aNew);
141 }
142 delete deExciteResult;
143 return &applyYourselfResult;
144 }
145
146 G4ReactionProductVector* G4AblaInterface::DeExcite(G4Fragment& aFragment)
147 {
148 if (!isInitialised)
149 InitialiseModel();
150
151 ablaResult->clear();
152
153 const G4int ARem = aFragment.GetA_asInt();
154 const G4int ZRem = aFragment.GetZ_asInt();
155 const G4int SRem = -aFragment.GetNumberOfLambdas(); // Strangeness = - (Number of lambdas)
156 const G4double eStarRem = aFragment.GetExcitationEnergy() / MeV;
157 const G4double jRem = aFragment.GetAngularMomentum().mag() / hbar_Planck;
158 const G4LorentzVector& pRem = aFragment.GetMomentum();
159 const G4double pxRem = pRem.x() / MeV;
160 const G4double pyRem = pRem.y() / MeV;
161 const G4double pzRem = pRem.z() / MeV;
162
163 ++eventNumber;
164
165 theABLAModel->DeexcitationAblaxx(ARem, ZRem, eStarRem, jRem, pxRem, pyRem, pzRem, (G4int)eventNumber, SRem);
166
167 G4ReactionProductVector* result = new G4ReactionProductVector;
168
169 for (G4int j = 0; j < ablaResult->ntrack; ++j)
170 { // Copy ABLA result to the EventInfo
171 G4ReactionProduct* product = toG4Particle(ablaResult->avv[j],
172 ablaResult->zvv[j],
173 ablaResult->svv[j],
174 ablaResult->enerj[j],
175 ablaResult->pxlab[j],
176 ablaResult->pylab[j],
177 ablaResult->pzlab[j]);
178 if (product)
179 {
180 product->SetCreatorModelID(secID);
181 result->push_back(product);
182 }
183 }
184 return result;
185 }
186
187 G4ParticleDefinition* G4AblaInterface::toG4ParticleDefinition(G4int A, G4int Z, G4int S) const
188 {
189 if (A == 1 && Z == 1 && S == 0)
190 return G4Proton::Proton();
191 else if (A == 1 && Z == 0 && S == 0)
192 return G4Neutron::Neutron();
193 else if (A == 1 && Z == 0 && S == -1)
194 return G4Lambda::Lambda();
195 else if (A == -1 && Z == 1 && S == 0)
196 return G4PionPlus::PionPlus();
197 else if (A == -1 && Z == -1 && S == 0)
198 return G4PionMinus::PionMinus();
199 else if (A == -1 && Z == 0 && S == 0)
200 return G4PionZero::PionZero();
201 else if (A == 0 && Z == 0 && S == 0)
202 return G4Gamma::Gamma();
203 else if (A == 2 && Z == 1 && S == 0)
204 return G4Deuteron::Deuteron();
205 else if (A == 3 && Z == 1 && S == 0)
206 return G4Triton::Triton();
207 else if (A == 3 && Z == 2 && S == 0)
208 return G4He3::He3();
209 else if (A == 3 && Z == 1 && S == -1)
210 return G4HyperTriton::Definition();
211 else if (A == 4 && Z == 2 && S == 0)
212 return G4Alpha::Alpha();
213 else if (A == 4 && Z == 1 && S == -1)
214 return G4HyperH4::Definition();
215 else if (A == 4 && Z == 2 && S == -1)
216 return G4HyperAlpha::Definition();
217 else if (A == 4 && Z == 1 && S == -2)
218 return G4DoubleHyperH4::Definition();
219 else if (A == 4 && Z == 0 && S == -2)
220 return G4DoubleHyperDoubleNeutron::Definition();
221 else if (A == 5 && Z == 2 && S == -1)
222 return G4HyperHe5::Definition();
223 else if (A > 0 && Z > 0 && A > Z)
224 { // Returns ground state ion definition.
225 auto ionfromtable = G4IonTable::GetIonTable()->GetIon(Z, A, std::abs(S), 0); // S is the number of lambdas
226 if (ionfromtable)
227 return ionfromtable;
228 else
229 {
230 G4cout << "Can't convert particle with A=" << A << ", Z=" << Z << ", S=" << S
231 << " to G4ParticleDefinition, trouble ahead" << G4endl;
232 return 0;
233 }
234 }
235 else
236 { // Error, unrecognized particle
237 G4cout << "Can't convert particle with A=" << A << ", Z=" << Z << ", S=" << S
238 << " to G4ParticleDefinition, trouble ahead" << G4endl;
239 return 0;
240 }
241 }
242
243 G4ReactionProduct*
244 G4AblaInterface::toG4Particle(G4int A, G4int Z, G4int S, G4double kinE, G4double px, G4double py, G4double pz) const
245 {
246 G4ParticleDefinition* def = toG4ParticleDefinition(A, Z, S);
247 if (def == 0)
248 { // Check if we have a valid particle definition
249 return 0;
250 }
251
252 const G4double energy = kinE * MeV;
253 const G4ThreeVector momentum(px, py, pz);
254 const G4ThreeVector momentumDirection = momentum.unit();
255 G4DynamicParticle p(def, momentumDirection, energy);
256 G4ReactionProduct* r = new G4ReactionProduct(def);
257 (*r) = p;
258 return r;
259 }
260
261 void G4AblaInterface::ModelDescription(std::ostream& outFile) const
262 {
263 outFile << "ABLA++ does not provide an implementation of the ApplyYourself "
264 "method!\n\n";
265 }
266
267 void G4AblaInterface::DeExciteModelDescription(std::ostream& outFile) const
268 {
269 outFile << "ABLA++ is a statistical model for nuclear de-excitation. It simulates\n"
270 << "the gamma emission and the evaporation of neutrons, light charged\n"
271 << "particles and IMFs, as well as fission where applicable. The code\n"
272 << "included in Geant4 is a C++ translation of the original Fortran\n"
273 << "code ABLA07. Although the model has been recently extended to\n"
274 << "hypernuclei by including the evaporation of lambda particles.\n"
275 << "More details about the physics are available in the Geant4\n"
276 << "Physics Reference Manual and in the reference articles.\n\n"
277 << "References:\n"
278 << "(1) A. Kelic, M. V. Ricciardi, and K. H. Schmidt, in Proceedings of Joint\n"
279 << "ICTP-IAEA Advanced Workshop on Model Codes for Spallation Reactions,\n"
280 << "ICTP Trieste, Italy, 4–8 February 2008, edited by D. Filges, S. "
281 "Leray, Y. Yariv, A. Mengoni, A. Stanculescu, and G. Mank (IAEA "
282 "INDC(NDS)-530, Vienna, 2008), pp. 181–221.\n\n"
283 << "(2) J.L. Rodriguez-Sanchez, J.-C. David et al., Phys. Rev. C 98, 021602R (2018)\n"
284 << "(3) J.L. Rodriguez-Sanchez et al., Phys. Rev. C 105, 014623 (2022)\n"
285 << "(4) J.L. Rodriguez-Sanchez et al., Phys. Rev. Lett. 130, 132501 (2023)\n\n";
286 }
287