| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // Physics model class G4NeutronRadCaptureHP 27 // Physics model class G4NeutronRadCaptureHP
28 // derived from G4NeutronRadCapture 28 // derived from G4NeutronRadCapture
29 // 29 //
30 // Created: 02 October 2023 30 // Created: 02 October 2023
31 // Author V.Ivanchenko 31 // Author V.Ivanchenko
32 // 32 //
33 // 33 //
34 34
35 #include "G4NeutronRadCaptureHP.hh" 35 #include "G4NeutronRadCaptureHP.hh"
36 #include "G4HadronicInteractionRegistry.hh" 36 #include "G4HadronicInteractionRegistry.hh"
37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
38 #include "G4ParticleDefinition.hh" 38 #include "G4ParticleDefinition.hh"
39 #include "G4Fragment.hh" 39 #include "G4Fragment.hh"
40 #include "G4FragmentVector.hh" 40 #include "G4FragmentVector.hh"
41 #include "G4NucleiProperties.hh" 41 #include "G4NucleiProperties.hh"
42 #include "G4VPreCompoundModel.hh" 42 #include "G4VPreCompoundModel.hh"
43 #include "G4ExcitationHandler.hh" 43 #include "G4ExcitationHandler.hh"
44 #include "G4VEvaporationChannel.hh" 44 #include "G4VEvaporationChannel.hh"
45 #include "G4PhotonEvaporation.hh" 45 #include "G4PhotonEvaporation.hh"
46 #include "G4DynamicParticle.hh" 46 #include "G4DynamicParticle.hh"
47 #include "G4ParticleTable.hh" 47 #include "G4ParticleTable.hh"
48 #include "G4ParticleHPManager.hh" 48 #include "G4ParticleHPManager.hh"
49 #include "G4IonTable.hh" 49 #include "G4IonTable.hh"
50 #include "G4Electron.hh" 50 #include "G4Electron.hh"
51 #include "G4Deuteron.hh" 51 #include "G4Deuteron.hh"
52 #include "G4Triton.hh" 52 #include "G4Triton.hh"
53 #include "G4He3.hh" 53 #include "G4He3.hh"
54 #include "G4Alpha.hh" 54 #include "G4Alpha.hh"
55 #include "Randomize.hh" 55 #include "Randomize.hh"
56 #include "G4RandomDirection.hh" 56 #include "G4RandomDirection.hh"
57 #include "G4HadronicParameters.hh" 57 #include "G4HadronicParameters.hh"
58 #include "G4PhysicsModelCatalog.hh" 58 #include "G4PhysicsModelCatalog.hh"
59 59
60 G4NeutronRadCaptureHP::G4NeutronRadCaptureHP() 60 G4NeutronRadCaptureHP::G4NeutronRadCaptureHP()
61 : G4HadronicInteraction("nRadCaptureHP"), 61 : G4HadronicInteraction("nRadCaptureHP"),
62 electron(G4Electron::Electron()), 62 electron(G4Electron::Electron()),
63 fManagerHP(G4ParticleHPManager::GetInstanc 63 fManagerHP(G4ParticleHPManager::GetInstance()),
64 lowestEnergyLimit(1.0e-11*CLHEP::eV), 64 lowestEnergyLimit(1.0e-11*CLHEP::eV),
65 minExcitation(0.1*CLHEP::keV), 65 minExcitation(0.1*CLHEP::keV),
66 emax(20*CLHEP::MeV), 66 emax(20*CLHEP::MeV),
67 emaxT(fManagerHP->GetMaxEnergyDoppler()), 67 emaxT(fManagerHP->GetMaxEnergyDoppler()),
68 lab4mom(0.,0.,0.,0.) 68 lab4mom(0.,0.,0.,0.)
69 { 69 {
70 SetMaxEnergy( G4HadronicParameters::Instance 70 SetMaxEnergy( G4HadronicParameters::Instance()->GetMaxEnergy() );
71 theTableOfIons = G4ParticleTable::GetParticl 71 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
72 } 72 }
73 73
74 G4NeutronRadCaptureHP::~G4NeutronRadCaptureHP( 74 G4NeutronRadCaptureHP::~G4NeutronRadCaptureHP()
75 { 75 {
76 if (fLocalPE) { delete photonEvaporation; } 76 if (fLocalPE) { delete photonEvaporation; }
77 } 77 }
78 78
79 void G4NeutronRadCaptureHP::BuildPhysicsTable( 79 void G4NeutronRadCaptureHP::BuildPhysicsTable(const G4ParticleDefinition&)
80 { 80 {
81 if (photonEvaporation != nullptr) { return; 81 if (photonEvaporation != nullptr) { return; }
82 G4HadronicInteraction* p = 82 G4HadronicInteraction* p =
83 G4HadronicInteractionRegistry::Instance()- 83 G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
84 if (nullptr != p) { 84 if (nullptr != p) {
85 auto handler = 85 auto handler =
86 (static_cast<G4VPreCompoundModel*>(p))-> 86 (static_cast<G4VPreCompoundModel*>(p))->GetExcitationHandler();
87 if (nullptr != handler) 87 if (nullptr != handler)
88 photonEvaporation = handler->GetPhotonEv 88 photonEvaporation = handler->GetPhotonEvaporation();
89 } 89 }
90 G4DeexPrecoParameters* param = 90 G4DeexPrecoParameters* param =
91 G4NuclearLevelData::GetInstance()->GetPara 91 G4NuclearLevelData::GetInstance()->GetParameters();
92 minExcitation = param->GetMinExcitation(); 92 minExcitation = param->GetMinExcitation();
93 icID = G4PhysicsModelCatalog::GetModelID("mo 93 icID = G4PhysicsModelCatalog::GetModelID("model_e-InternalConversion");
94 secID = G4PhysicsModelCatalog::GetModelID("m 94 secID = G4PhysicsModelCatalog::GetModelID("model_" + GetModelName());
95 if (nullptr == photonEvaporation) { 95 if (nullptr == photonEvaporation) {
96 photonEvaporation = new G4PhotonEvaporatio 96 photonEvaporation = new G4PhotonEvaporation();
97 fLocalPE = true; 97 fLocalPE = true;
98 } 98 }
99 photonEvaporation->Initialise(); 99 photonEvaporation->Initialise();
100 photonEvaporation->SetICM(true); 100 photonEvaporation->SetICM(true);
101 } 101 }
102 102
103 G4HadFinalState* G4NeutronRadCaptureHP::ApplyY 103 G4HadFinalState* G4NeutronRadCaptureHP::ApplyYourself(
104 const G4HadProjectile& aTrack, G4Nucleus& 104 const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
105 { 105 {
106 theParticleChange.Clear(); 106 theParticleChange.Clear();
107 G4double ekin = aTrack.GetKineticEnergy(); 107 G4double ekin = aTrack.GetKineticEnergy();
108 if (ekin > emax) { 108 if (ekin > emax) {
109 return &theParticleChange; 109 return &theParticleChange;
110 } 110 }
111 111
112 theParticleChange.SetStatusChange(stopAndKil 112 theParticleChange.SetStatusChange(stopAndKill);
113 G4double T = aTrack.GetMaterial()->GetTemper 113 G4double T = aTrack.GetMaterial()->GetTemperature();
114 114
115 G4int A = theNucleus.GetA_asInt(); 115 G4int A = theNucleus.GetA_asInt();
116 G4int Z = theNucleus.GetZ_asInt(); 116 G4int Z = theNucleus.GetZ_asInt();
117 117
118 G4double time = aTrack.GetGlobalTime(); 118 G4double time = aTrack.GetGlobalTime();
119 119
120 // Create initial state 120 // Create initial state
121 G4double mass = G4NucleiProperties::GetNucle 121 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
122 122
123 // no Doppler broading 123 // no Doppler broading
124 G4double factT = T/CLHEP::STP_Temperature; 124 G4double factT = T/CLHEP::STP_Temperature;
125 125
126 if (ekin >= emaxT*factT || fManagerHP->GetNe 126 if (ekin >= emaxT*factT || fManagerHP->GetNeglectDoppler()) {
127 lab4mom.set(0.,0.,0.,mass); 127 lab4mom.set(0.,0.,0.,mass);
128 } else { 128 } else {
129 G4double lambda = 1.0/(CLHEP::k_Boltzmann* 129 G4double lambda = 1.0/(CLHEP::k_Boltzmann*T);
130 G4double erand = G4RandGamma::shoot(2.0, l 130 G4double erand = G4RandGamma::shoot(2.0, lambda);
131 auto mom = G4RandomDirection()*std::sqrt(2 131 auto mom = G4RandomDirection()*std::sqrt(2*mass*erand);
132 lab4mom.set(mom.x(), mom.y(), mom.z(), mas 132 lab4mom.set(mom.x(), mom.y(), mom.z(), mass + erand);
133 } 133 }
134 134
135 lab4mom += aTrack.Get4Momentum(); 135 lab4mom += aTrack.Get4Momentum();
136 136
137 G4double M = lab4mom.mag(); 137 G4double M = lab4mom.mag();
138 ++A; 138 ++A;
139 mass = G4NucleiProperties::GetNuclearMass(A, 139 mass = G4NucleiProperties::GetNuclearMass(A, Z);
140 //G4cout << "Capture start: Z= " << Z << " A 140 //G4cout << "Capture start: Z= " << Z << " A= " << A
141 // << " LabM= " << M << " Mcompound= " << 141 // << " LabM= " << M << " Mcompound= " << mass << G4endl;
142 142
143 // simplified method of 1 gamma emission 143 // simplified method of 1 gamma emission
144 if (A <= 4) { 144 if (A <= 4) {
145 145
146 if (verboseLevel > 1) { 146 if (verboseLevel > 1) {
147 G4cout << "G4NeutronRadCaptureHP::DoIt: 147 G4cout << "G4NeutronRadCaptureHP::DoIt: Eini(MeV)="
148 << ekin/MeV << " Eexc(MeV)= " 148 << ekin/MeV << " Eexc(MeV)= "
149 << (M - mass)/MeV 149 << (M - mass)/MeV
150 << " Z= " << Z << " A= " << A << G4en 150 << " Z= " << Z << " A= " << A << G4endl;
151 } 151 }
152 if (M - mass > lowestEnergyLimit) { 152 if (M - mass > lowestEnergyLimit) {
153 G4ThreeVector bst = lab4mom.boostVector( 153 G4ThreeVector bst = lab4mom.boostVector();
154 G4double e1 = (M - mass)*(M + mass)/(2*M 154 G4double e1 = (M - mass)*(M + mass)/(2*M);
155 G4LorentzVector lv2(e1*G4RandomDirection 155 G4LorentzVector lv2(e1*G4RandomDirection(),e1);
156 lv2.boost(bst); 156 lv2.boost(bst);
157 if (verboseLevel > 1) { 157 if (verboseLevel > 1) {
158 G4cout << "Gamma 4-mom: " << lv2 << " 158 G4cout << "Gamma 4-mom: " << lv2 << " Escm(MeV)=" << e1/CLHEP::MeV << G4endl;
159 } 159 }
160 lab4mom -= lv2; 160 lab4mom -= lv2;
161 G4HadSecondary* news = 161 G4HadSecondary* news =
162 new G4HadSecondary(new G4DynamicParticle(G4G 162 new G4HadSecondary(new G4DynamicParticle(G4Gamma::Gamma(), lv2));
163 news->SetTime(time); 163 news->SetTime(time);
164 news->SetCreatorModelID(secID); 164 news->SetCreatorModelID(secID);
165 theParticleChange.AddSecondary(*news); 165 theParticleChange.AddSecondary(*news);
166 delete news; 166 delete news;
167 } 167 }
168 168
169 const G4ParticleDefinition* theDef = nullp 169 const G4ParticleDefinition* theDef = nullptr;
170 170
171 if (Z == 1 && A == 2) {theDef = G4Deu 171 if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron(); }
172 else if (Z == 1 && A == 3) {theDef = G4Tri 172 else if (Z == 1 && A == 3) {theDef = G4Triton::Triton(); }
173 else if (Z == 2 && A == 3) {theDef = G4He3 173 else if (Z == 2 && A == 3) {theDef = G4He3::He3(); }
174 else if (Z == 2 && A == 4) {theDef = G4Alp 174 else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha(); }
175 else { theDef = theTableOfIons->GetIon(Z, 175 else { theDef = theTableOfIons->GetIon(Z, A, 0.0, noFloat, 0); }
176 176
177 if (nullptr != theDef) { 177 if (nullptr != theDef) {
178 G4HadSecondary* news = 178 G4HadSecondary* news =
179 new G4HadSecondary(new G4DynamicPartic 179 new G4HadSecondary(new G4DynamicParticle(theDef, lab4mom));
180 news->SetTime(time); 180 news->SetTime(time);
181 news->SetCreatorModelID(secID); 181 news->SetCreatorModelID(secID);
182 theParticleChange.AddSecondary(*news); 182 theParticleChange.AddSecondary(*news);
183 delete news; 183 delete news;
184 } 184 }
185 185
186 // Use photon evaporation 186 // Use photon evaporation
187 } else { 187 } else {
188 188
189 // protection against wrong kinematic 189 // protection against wrong kinematic
190 if (M < mass) { 190 if (M < mass) {
191 G4double etot = std::max(mass, lab4mom.e 191 G4double etot = std::max(mass, lab4mom.e());
192 G4double ptot = std::sqrt((etot - mass)* 192 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
193 G4ThreeVector v = lab4mom.vect().unit(); 193 G4ThreeVector v = lab4mom.vect().unit();
194 lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()* 194 lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*ptot,etot);
195 } 195 }
196 196
197 G4Fragment* aFragment = new G4Fragment(A, 197 G4Fragment* aFragment = new G4Fragment(A, Z, lab4mom);
198 198
199 if (verboseLevel > 1) { 199 if (verboseLevel > 1) {
200 G4cout << "G4NeutronRadCaptureHP::ApplyY 200 G4cout << "G4NeutronRadCaptureHP::ApplyYourself initial G4Fragmet:"
201 << G4endl; 201 << G4endl;
202 G4cout << aFragment << G4endl; 202 G4cout << aFragment << G4endl;
203 } 203 }
204 204
205 // 205 //
206 // Sample final state 206 // Sample final state
207 // 207 //
208 G4FragmentVector* fv = photonEvaporation-> 208 G4FragmentVector* fv = photonEvaporation->BreakUpFragment(aFragment);
209 if (nullptr == fv) { fv = new G4FragmentVe 209 if (nullptr == fv) { fv = new G4FragmentVector(); }
210 fv->push_back(aFragment); 210 fv->push_back(aFragment);
211 211
212 if (verboseLevel > 1) { 212 if (verboseLevel > 1) {
213 G4cout << "G4NeutronRadCaptureHP: " << f 213 G4cout << "G4NeutronRadCaptureHP: " << fv->size() << " final particle icID= "
214 << icID << G4endl; 214 << icID << G4endl;
215 } 215 }
216 for (auto const & f : *fv) { 216 for (auto const & f : *fv) {
217 G4double etot = f->GetMomentum().e(); 217 G4double etot = f->GetMomentum().e();
218 218
219 Z = f->GetZ_asInt(); 219 Z = f->GetZ_asInt();
220 A = f->GetA_asInt(); 220 A = f->GetA_asInt();
221 221
222 const G4ParticleDefinition* theDef; 222 const G4ParticleDefinition* theDef;
223 if (0 == Z && 0 == A) { theDef = f->Get 223 if (0 == Z && 0 == A) { theDef = f->GetParticleDefinition(); }
224 else if (Z == 1 && A == 2) { theDef = G4 224 else if (Z == 1 && A == 2) { theDef = G4Deuteron::Deuteron(); }
225 else if (Z == 1 && A == 3) { theDef = G4 225 else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
226 else if (Z == 2 && A == 3) { theDef = G4 226 else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
227 else if (Z == 2 && A == 4) { theDef = G4 227 else if (Z == 2 && A == 4) { theDef = G4Alpha::Alpha(); }
228 else { 228 else {
229 G4double eexc = f->GetExcitationEnergy 229 G4double eexc = f->GetExcitationEnergy();
230 if (eexc <= minExcitation) { eexc = 0.0; } 230 if (eexc <= minExcitation) { eexc = 0.0; }
231 theDef = theTableOfIons->GetIon(Z, A, eexc, 231 theDef = theTableOfIons->GetIon(Z, A, eexc, noFloat, 0);
232 /* 232 /*
233 G4cout << "### NC Find ion Z= " << Z << " A= 233 G4cout << "### NC Find ion Z= " << Z << " A= " << A
234 << " Eexc(MeV)= " << eexc/MeV << " " 234 << " Eexc(MeV)= " << eexc/MeV << " "
235 << theDef << G4endl; 235 << theDef << G4endl;
236 */ 236 */
237 } 237 }
238 ekin = std::max(0.0, etot - theDef->GetP 238 ekin = std::max(0.0, etot - theDef->GetPDGMass());
239 if (verboseLevel > 1) { 239 if (verboseLevel > 1) {
240 G4cout << theDef->GetParticleName() 240 G4cout << theDef->GetParticleName()
241 << " Ekin(MeV)= " << ekin/MeV 241 << " Ekin(MeV)= " << ekin/MeV
242 << " p: " << f->GetMomentum().vect() 242 << " p: " << f->GetMomentum().vect()
243 << G4endl; 243 << G4endl;
244 } 244 }
245 G4HadSecondary* news = 245 G4HadSecondary* news =
246 new G4HadSecondary(new G4DynamicPartic 246 new G4HadSecondary(new G4DynamicParticle(theDef,
247 247 f->GetMomentum().vect().unit(),
248 248 ekin));
249 G4double timeF = std::max(f->GetCreation 249 G4double timeF = std::max(f->GetCreationTime(), 0.0);
250 news->SetTime(time + timeF); 250 news->SetTime(time + timeF);
251 if (theDef == electron) { 251 if (theDef == electron) {
252 news->SetCreatorModelID(icID); 252 news->SetCreatorModelID(icID);
253 } else { 253 } else {
254 news->SetCreatorModelID(secID); 254 news->SetCreatorModelID(secID);
255 } 255 }
256 theParticleChange.AddSecondary(*news); 256 theParticleChange.AddSecondary(*news);
257 delete news; 257 delete news;
258 delete f; 258 delete f;
259 } 259 }
260 delete fv; 260 delete fv;
261 } 261 }
262 //G4cout << "Capture done" << G4endl; 262 //G4cout << "Capture done" << G4endl;
263 return &theParticleChange; 263 return &theParticleChange;
264 } 264 }
265 265
266 266