| 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 // Physics model class G4NeutronRadCaptureHP
28 // derived from G4NeutronRadCapture
29 //
30 // Created: 02 October 2023
31 // Author V.Ivanchenko
32 //
33 //
34
35 #include "G4NeutronRadCaptureHP.hh"
36 #include "G4HadronicInteractionRegistry.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4Fragment.hh"
40 #include "G4FragmentVector.hh"
41 #include "G4NucleiProperties.hh"
42 #include "G4VPreCompoundModel.hh"
43 #include "G4ExcitationHandler.hh"
44 #include "G4VEvaporationChannel.hh"
45 #include "G4PhotonEvaporation.hh"
46 #include "G4DynamicParticle.hh"
47 #include "G4ParticleTable.hh"
48 #include "G4ParticleHPManager.hh"
49 #include "G4IonTable.hh"
50 #include "G4Electron.hh"
51 #include "G4Deuteron.hh"
52 #include "G4Triton.hh"
53 #include "G4He3.hh"
54 #include "G4Alpha.hh"
55 #include "Randomize.hh"
56 #include "G4RandomDirection.hh"
57 #include "G4HadronicParameters.hh"
58 #include "G4PhysicsModelCatalog.hh"
59
60 G4NeutronRadCaptureHP::G4NeutronRadCaptureHP()
61 : G4HadronicInteraction("nRadCaptureHP"),
62 electron(G4Electron::Electron()),
63 fManagerHP(G4ParticleHPManager::GetInstance()),
64 lowestEnergyLimit(1.0e-11*CLHEP::eV),
65 minExcitation(0.1*CLHEP::keV),
66 emax(20*CLHEP::MeV),
67 emaxT(fManagerHP->GetMaxEnergyDoppler()),
68 lab4mom(0.,0.,0.,0.)
69 {
70 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
71 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
72 }
73
74 G4NeutronRadCaptureHP::~G4NeutronRadCaptureHP()
75 {
76 if (fLocalPE) { delete photonEvaporation; }
77 }
78
79 void G4NeutronRadCaptureHP::BuildPhysicsTable(const G4ParticleDefinition&)
80 {
81 if (photonEvaporation != nullptr) { return; }
82 G4HadronicInteraction* p =
83 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
84 if (nullptr != p) {
85 auto handler =
86 (static_cast<G4VPreCompoundModel*>(p))->GetExcitationHandler();
87 if (nullptr != handler)
88 photonEvaporation = handler->GetPhotonEvaporation();
89 }
90 G4DeexPrecoParameters* param =
91 G4NuclearLevelData::GetInstance()->GetParameters();
92 minExcitation = param->GetMinExcitation();
93 icID = G4PhysicsModelCatalog::GetModelID("model_e-InternalConversion");
94 secID = G4PhysicsModelCatalog::GetModelID("model_" + GetModelName());
95 if (nullptr == photonEvaporation) {
96 photonEvaporation = new G4PhotonEvaporation();
97 fLocalPE = true;
98 }
99 photonEvaporation->Initialise();
100 photonEvaporation->SetICM(true);
101 }
102
103 G4HadFinalState* G4NeutronRadCaptureHP::ApplyYourself(
104 const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
105 {
106 theParticleChange.Clear();
107 G4double ekin = aTrack.GetKineticEnergy();
108 if (ekin > emax) {
109 return &theParticleChange;
110 }
111
112 theParticleChange.SetStatusChange(stopAndKill);
113 G4double T = aTrack.GetMaterial()->GetTemperature();
114
115 G4int A = theNucleus.GetA_asInt();
116 G4int Z = theNucleus.GetZ_asInt();
117
118 G4double time = aTrack.GetGlobalTime();
119
120 // Create initial state
121 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
122
123 // no Doppler broading
124 G4double factT = T/CLHEP::STP_Temperature;
125
126 if (ekin >= emaxT*factT || fManagerHP->GetNeglectDoppler()) {
127 lab4mom.set(0.,0.,0.,mass);
128 } else {
129 G4double lambda = 1.0/(CLHEP::k_Boltzmann*T);
130 G4double erand = G4RandGamma::shoot(2.0, lambda);
131 auto mom = G4RandomDirection()*std::sqrt(2*mass*erand);
132 lab4mom.set(mom.x(), mom.y(), mom.z(), mass + erand);
133 }
134
135 lab4mom += aTrack.Get4Momentum();
136
137 G4double M = lab4mom.mag();
138 ++A;
139 mass = G4NucleiProperties::GetNuclearMass(A, Z);
140 //G4cout << "Capture start: Z= " << Z << " A= " << A
141 // << " LabM= " << M << " Mcompound= " << mass << G4endl;
142
143 // simplified method of 1 gamma emission
144 if (A <= 4) {
145
146 if (verboseLevel > 1) {
147 G4cout << "G4NeutronRadCaptureHP::DoIt: Eini(MeV)="
148 << ekin/MeV << " Eexc(MeV)= "
149 << (M - mass)/MeV
150 << " Z= " << Z << " A= " << A << G4endl;
151 }
152 if (M - mass > lowestEnergyLimit) {
153 G4ThreeVector bst = lab4mom.boostVector();
154 G4double e1 = (M - mass)*(M + mass)/(2*M);
155 G4LorentzVector lv2(e1*G4RandomDirection(),e1);
156 lv2.boost(bst);
157 if (verboseLevel > 1) {
158 G4cout << "Gamma 4-mom: " << lv2 << " Escm(MeV)=" << e1/CLHEP::MeV << G4endl;
159 }
160 lab4mom -= lv2;
161 G4HadSecondary* news =
162 new G4HadSecondary(new G4DynamicParticle(G4Gamma::Gamma(), lv2));
163 news->SetTime(time);
164 news->SetCreatorModelID(secID);
165 theParticleChange.AddSecondary(*news);
166 delete news;
167 }
168
169 const G4ParticleDefinition* theDef = nullptr;
170
171 if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron(); }
172 else if (Z == 1 && A == 3) {theDef = G4Triton::Triton(); }
173 else if (Z == 2 && A == 3) {theDef = G4He3::He3(); }
174 else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha(); }
175 else { theDef = theTableOfIons->GetIon(Z, A, 0.0, noFloat, 0); }
176
177 if (nullptr != theDef) {
178 G4HadSecondary* news =
179 new G4HadSecondary(new G4DynamicParticle(theDef, lab4mom));
180 news->SetTime(time);
181 news->SetCreatorModelID(secID);
182 theParticleChange.AddSecondary(*news);
183 delete news;
184 }
185
186 // Use photon evaporation
187 } else {
188
189 // protection against wrong kinematic
190 if (M < mass) {
191 G4double etot = std::max(mass, lab4mom.e());
192 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
193 G4ThreeVector v = lab4mom.vect().unit();
194 lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*ptot,etot);
195 }
196
197 G4Fragment* aFragment = new G4Fragment(A, Z, lab4mom);
198
199 if (verboseLevel > 1) {
200 G4cout << "G4NeutronRadCaptureHP::ApplyYourself initial G4Fragmet:"
201 << G4endl;
202 G4cout << aFragment << G4endl;
203 }
204
205 //
206 // Sample final state
207 //
208 G4FragmentVector* fv = photonEvaporation->BreakUpFragment(aFragment);
209 if (nullptr == fv) { fv = new G4FragmentVector(); }
210 fv->push_back(aFragment);
211
212 if (verboseLevel > 1) {
213 G4cout << "G4NeutronRadCaptureHP: " << fv->size() << " final particle icID= "
214 << icID << G4endl;
215 }
216 for (auto const & f : *fv) {
217 G4double etot = f->GetMomentum().e();
218
219 Z = f->GetZ_asInt();
220 A = f->GetA_asInt();
221
222 const G4ParticleDefinition* theDef;
223 if (0 == Z && 0 == A) { theDef = f->GetParticleDefinition(); }
224 else if (Z == 1 && A == 2) { theDef = G4Deuteron::Deuteron(); }
225 else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
226 else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
227 else if (Z == 2 && A == 4) { theDef = G4Alpha::Alpha(); }
228 else {
229 G4double eexc = f->GetExcitationEnergy();
230 if (eexc <= minExcitation) { eexc = 0.0; }
231 theDef = theTableOfIons->GetIon(Z, A, eexc, noFloat, 0);
232 /*
233 G4cout << "### NC Find ion Z= " << Z << " A= " << A
234 << " Eexc(MeV)= " << eexc/MeV << " "
235 << theDef << G4endl;
236 */
237 }
238 ekin = std::max(0.0, etot - theDef->GetPDGMass());
239 if (verboseLevel > 1) {
240 G4cout << theDef->GetParticleName()
241 << " Ekin(MeV)= " << ekin/MeV
242 << " p: " << f->GetMomentum().vect()
243 << G4endl;
244 }
245 G4HadSecondary* news =
246 new G4HadSecondary(new G4DynamicParticle(theDef,
247 f->GetMomentum().vect().unit(),
248 ekin));
249 G4double timeF = std::max(f->GetCreationTime(), 0.0);
250 news->SetTime(time + timeF);
251 if (theDef == electron) {
252 news->SetCreatorModelID(icID);
253 } else {
254 news->SetCreatorModelID(secID);
255 }
256 theParticleChange.AddSecondary(*news);
257 delete news;
258 delete f;
259 }
260 delete fv;
261 }
262 //G4cout << "Capture done" << G4endl;
263 return &theParticleChange;
264 }
265
266