| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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::GetInstanc
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
71 theTableOfIons = G4ParticleTable::GetParticl
72 }
73
74 G4NeutronRadCaptureHP::~G4NeutronRadCaptureHP(
75 {
76 if (fLocalPE) { delete photonEvaporation; }
77 }
78
79 void G4NeutronRadCaptureHP::BuildPhysicsTable(
80 {
81 if (photonEvaporation != nullptr) { return;
82 G4HadronicInteraction* p =
83 G4HadronicInteractionRegistry::Instance()-
84 if (nullptr != p) {
85 auto handler =
86 (static_cast<G4VPreCompoundModel*>(p))->
87 if (nullptr != handler)
88 photonEvaporation = handler->GetPhotonEv
89 }
90 G4DeexPrecoParameters* param =
91 G4NuclearLevelData::GetInstance()->GetPara
92 minExcitation = param->GetMinExcitation();
93 icID = G4PhysicsModelCatalog::GetModelID("mo
94 secID = G4PhysicsModelCatalog::GetModelID("m
95 if (nullptr == photonEvaporation) {
96 photonEvaporation = new G4PhotonEvaporatio
97 fLocalPE = true;
98 }
99 photonEvaporation->Initialise();
100 photonEvaporation->SetICM(true);
101 }
102
103 G4HadFinalState* G4NeutronRadCaptureHP::ApplyY
104 const G4HadProjectile& aTrack, G4Nucleus&
105 {
106 theParticleChange.Clear();
107 G4double ekin = aTrack.GetKineticEnergy();
108 if (ekin > emax) {
109 return &theParticleChange;
110 }
111
112 theParticleChange.SetStatusChange(stopAndKil
113 G4double T = aTrack.GetMaterial()->GetTemper
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::GetNucle
122
123 // no Doppler broading
124 G4double factT = T/CLHEP::STP_Temperature;
125
126 if (ekin >= emaxT*factT || fManagerHP->GetNe
127 lab4mom.set(0.,0.,0.,mass);
128 } else {
129 G4double lambda = 1.0/(CLHEP::k_Boltzmann*
130 G4double erand = G4RandGamma::shoot(2.0, l
131 auto mom = G4RandomDirection()*std::sqrt(2
132 lab4mom.set(mom.x(), mom.y(), mom.z(), mas
133 }
134
135 lab4mom += aTrack.Get4Momentum();
136
137 G4double M = lab4mom.mag();
138 ++A;
139 mass = G4NucleiProperties::GetNuclearMass(A,
140 //G4cout << "Capture start: Z= " << Z << " A
141 // << " LabM= " << M << " Mcompound= " <<
142
143 // simplified method of 1 gamma emission
144 if (A <= 4) {
145
146 if (verboseLevel > 1) {
147 G4cout << "G4NeutronRadCaptureHP::DoIt:
148 << ekin/MeV << " Eexc(MeV)= "
149 << (M - mass)/MeV
150 << " Z= " << Z << " A= " << A << G4en
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
156 lv2.boost(bst);
157 if (verboseLevel > 1) {
158 G4cout << "Gamma 4-mom: " << lv2 << "
159 }
160 lab4mom -= lv2;
161 G4HadSecondary* news =
162 new G4HadSecondary(new G4DynamicParticle(G4G
163 news->SetTime(time);
164 news->SetCreatorModelID(secID);
165 theParticleChange.AddSecondary(*news);
166 delete news;
167 }
168
169 const G4ParticleDefinition* theDef = nullp
170
171 if (Z == 1 && A == 2) {theDef = G4Deu
172 else if (Z == 1 && A == 3) {theDef = G4Tri
173 else if (Z == 2 && A == 3) {theDef = G4He3
174 else if (Z == 2 && A == 4) {theDef = G4Alp
175 else { theDef = theTableOfIons->GetIon(Z,
176
177 if (nullptr != theDef) {
178 G4HadSecondary* news =
179 new G4HadSecondary(new G4DynamicPartic
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)*
193 G4ThreeVector v = lab4mom.vect().unit();
194 lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*
195 }
196
197 G4Fragment* aFragment = new G4Fragment(A,
198
199 if (verboseLevel > 1) {
200 G4cout << "G4NeutronRadCaptureHP::ApplyY
201 << G4endl;
202 G4cout << aFragment << G4endl;
203 }
204
205 //
206 // Sample final state
207 //
208 G4FragmentVector* fv = photonEvaporation->
209 if (nullptr == fv) { fv = new G4FragmentVe
210 fv->push_back(aFragment);
211
212 if (verboseLevel > 1) {
213 G4cout << "G4NeutronRadCaptureHP: " << f
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->Get
224 else if (Z == 1 && A == 2) { theDef = G4
225 else if (Z == 1 && A == 3) { theDef = G4
226 else if (Z == 2 && A == 3) { theDef = G4
227 else if (Z == 2 && A == 4) { theDef = G4
228 else {
229 G4double eexc = f->GetExcitationEnergy
230 if (eexc <= minExcitation) { eexc = 0.0; }
231 theDef = theTableOfIons->GetIon(Z, A, eexc,
232 /*
233 G4cout << "### NC Find ion Z= " << Z << " A=
234 << " Eexc(MeV)= " << eexc/MeV << " "
235 << theDef << G4endl;
236 */
237 }
238 ekin = std::max(0.0, etot - theDef->GetP
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 G4DynamicPartic
247
248
249 G4double timeF = std::max(f->GetCreation
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