| 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 // -------------------------------------------------------------------
28 //
29 // GEANT4 source file
30 //
31 // File name: G4NuDEXNeutronCaptureModel
32 //
33 // Author: E.Mendoza & A.Ribon
34 //
35 // Creation date: 29 May 2024
36 //
37 // Description: This class (a proxy of the class G4NuDEX) uses
38 // the NuDEX model to produce gammas and internal
39 // conversion electrons from neutron capture.
40 // Whenever NuDEX is not applicable, G4PhotonEvaporation
41 // is used.
42 // The implementation of this class follows the code
43 // of the class G4NeutronRadCapture.
44 //
45 // Modifications:
46 //
47 // -------------------------------------------------------------------
48 //
49
50 #include "G4NuDEXNeutronCaptureModel.hh"
51 #include "G4NuDEXStatisticalNucleus.hh"
52 #include "Randomize.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4LorentzVector.hh"
56 #include "G4Gamma.hh"
57 #include "G4Electron.hh"
58 #include "G4Positron.hh"
59 #include "G4Deuteron.hh"
60 #include "G4Triton.hh"
61 #include "G4He3.hh"
62 #include "G4Alpha.hh"
63 #include "G4NucleiProperties.hh"
64 #include "G4IonTable.hh"
65 #include "G4ParticleTable.hh"
66 #include "G4HadronicParameters.hh"
67 #include "G4DeexPrecoParameters.hh"
68 #include "G4NuclearLevelData.hh"
69 #include "G4PhotonEvaporation.hh"
70 #include "G4PhysicsModelCatalog.hh"
71
72
73 G4NuDEXNeutronCaptureModel::G4NuDEXNeutronCaptureModel() : G4HadronicInteraction( "nuDEX_neutronCapture" ) {
74 for ( G4int i = 0; i < G4NUDEX_MAXZA; i++ ) {
75 theStatisticalNucleus[i] = nullptr;
76 HasData[i] = 0;
77 }
78 BrOption = -1;
79 BandWidth = 0;
80 NuDEXLibDirectory = "";
81 photonEvaporation = nullptr;
82 auto ch = G4FindDataDir( "G4NUDEXLIBDATA" );
83 if ( ch == nullptr ) {
84 G4Exception( "G4NuDEXNeutronCaptureModel()", "had0707", FatalException, "Environment variable G4NUDEXLIBDATA is not defined" );
85 } else {
86 NuDEXLibDirectory = G4String(ch);
87 }
88 }
89
90
91 void G4NuDEXNeutronCaptureModel::InitialiseModel() {
92 if ( photonEvaporation != nullptr ) return;
93 G4DeexPrecoParameters* param = G4NuclearLevelData::GetInstance()->GetParameters();
94 minExcitation = param->GetMinExcitation();
95 photonEvaporation = new G4PhotonEvaporation;
96 photonEvaporation->Initialise();
97 photonEvaporation->SetICM( true );
98 secID = G4PhysicsModelCatalog::GetModelID( "model_" + GetModelName() );
99 lowestEnergyLimit = 10.0*CLHEP::eV;
100 minExcitation = 0.1*CLHEP::keV;
101 }
102
103
104 G4NuDEXNeutronCaptureModel::~G4NuDEXNeutronCaptureModel(){
105 for ( G4int i = 0; i < G4NUDEX_MAXZA; i++ ) {
106 if ( theStatisticalNucleus[i] ) delete theStatisticalNucleus[i];
107 }
108 }
109
110
111 G4HadFinalState* G4NuDEXNeutronCaptureModel::ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &theNucleus ) {
112 theParticleChange.Clear();
113 theParticleChange.SetStatusChange( stopAndKill );
114 G4int A = theNucleus.GetA_asInt();
115 G4int Z = theNucleus.GetZ_asInt();
116 G4double time = aTrack.GetGlobalTime(); // Time in the lab frame
117 // Create initial state
118 G4LorentzVector lab4mom( 0.0, 0.0, 0.0, G4NucleiProperties::GetNuclearMass(A, Z) );
119 lab4mom += aTrack.Get4Momentum();
120 G4double systemMass = lab4mom.mag();
121 ++A; // Compound nucleus: target nucleus + neutron
122 G4double compoundMass = G4NucleiProperties::GetNuclearMass(A, Z);
123 // If the energy available is to small to do anything interesting, gives up
124 if ( systemMass - compoundMass <= lowestEnergyLimit ) return &theParticleChange;
125 G4ThreeVector boostFromCMtoLAB = lab4mom.boostVector();
126 G4double neutronEnergy = aTrack.GetKineticEnergy();
127
128 // Try to apply NuDEX
129 //G4int lspin = 0; // l-spin = 0, 1, 2 --> s-wave, p-wave, d-wave ...
130 //G4int jspinx2 = -1; // A negative value of jspinx2 means that is sampled according to the 2J+1 rule.
131 //G4int initialLevel = SelectInitialLevel( Z, A, neutronEnergy, lspin, jspinx2 );
132 G4int initialLevel = -1; // thermal neutron capture
133 std::vector< char > pType;
134 std::vector< G4double > pEnergy, pTime;
135 G4int npar = GenerateNeutronCaptureCascade( Z, A, neutronEnergy, initialLevel, pType, pEnergy, pTime );
136 if ( npar > 0 ) { // NuDEX can be applied
137
138 G4LorentzVector remainingLab4mom = lab4mom;
139 G4double latestEmission = time;
140 // Loop over the EM particles produced by 'GenerateNeutronCaptureCascade' and add them to the
141 // theParticleChange as secondaries. These particles are produced by NuDEX in the nucleus' rest-frame.
142 for ( G4int i = 0; i < npar; i++ ) {
143 G4ParticleDefinition* particleDef = nullptr;
144 if ( pType.at(i) == 'g' ) {
145 particleDef = G4Gamma::Definition();
146 } else if (pType.at(i) == 'e' ) {
147 particleDef = G4Electron::Definition();
148 } else if ( pType.at(i) == 'p' ) {
149 particleDef = G4Positron::Definition();
150 } else {
151 G4Exception( "G4NUDEXNeutronCaptureModel::ApplyYourself()", "had0707", FatalException, "Unknown particle type" );
152 }
153 G4double phi = G4UniformRand()*twopi;
154 G4double costheta = 2.0*G4UniformRand() - 1.0;
155 G4double sintheta = std::sqrt( 1.0 - costheta*costheta );
156 G4ThreeVector direction( sintheta*std::cos(phi), sintheta*std::sin(phi), costheta );
157 G4double mass = particleDef->GetPDGMass();
158 G4double particle3momMod = std::sqrt( pEnergy.at(i) * ( pEnergy.at(i) + 2.0*mass ) );
159 G4LorentzVector particle4mom( particle3momMod*direction, mass + pEnergy.at(i) ); // In the center-of-mass frame
160 particle4mom.boost( boostFromCMtoLAB ); // Now in the Lab frame
161 G4HadSecondary* secondary = new G4HadSecondary( new G4DynamicParticle( particleDef, particle4mom ) );
162 remainingLab4mom -= particle4mom;
163 // For simplicity, we neglect below the different frames of time (Lab) and pTime (center-of-mass)
164 secondary->SetTime( time + pTime.at(i) );
165 if ( latestEmission < time + pTime.at(i) ) latestEmission = time + pTime.at(i);
166 secondary->SetCreatorModelID( secID );
167 theParticleChange.AddSecondary( *secondary );
168 delete secondary;
169 }
170 // Treat now the residual nucleus (which is neglected by NuDEX)
171 const G4ParticleDefinition* resNuclDef = nullptr;
172 if ( Z == 1 && A == 2 ) resNuclDef = G4Deuteron::Definition();
173 else if ( Z == 1 && A == 3 ) resNuclDef = G4Triton::Definition();
174 else if ( Z == 2 && A == 3 ) resNuclDef = G4He3::Definition();
175 else if ( Z == 2 && A == 4 ) resNuclDef = G4Alpha::Alpha();
176 else resNuclDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon( Z, A, 0.0, noFloat, 0 );
177 if ( resNuclDef ) {
178 // To conserve energy-momentum, remainingLab4mom should be the Lorentz 4-momentum of the residual nucleus.
179 // Imposing the mass 'compoundMass' to the residual nucleus, and trying to conserve the total energy
180 // while keeping as low as possible the violation of the 3-momentum; in the case that the total energy
181 // cannot be conserved, the residual nucleus is produced at rest (in the Lab frame).
182 G4double resNuclLabEkin = std::max( remainingLab4mom.e() - compoundMass, 0.0 );
183 G4double resNuclLab3momMod = 0.0;
184 G4ThreeVector resNuclLabDir( 0.0, 0.0, 0.0 );
185 if ( resNuclLabEkin > 0.0 ) {
186 resNuclLab3momMod = std::sqrt( resNuclLabEkin * ( resNuclLabEkin + 2.0*compoundMass ) );
187 resNuclLabDir = remainingLab4mom.vect().unit();
188 }
189 G4LorentzVector resNuclLab4mom( resNuclLab3momMod*resNuclLabDir, resNuclLabEkin + compoundMass );
190 G4HadSecondary* secondary = new G4HadSecondary( new G4DynamicParticle( resNuclDef, resNuclLab4mom ) );
191 secondary->SetTime( latestEmission );
192 secondary->SetCreatorModelID( secID );
193 theParticleChange.AddSecondary( *secondary );
194 delete secondary;
195 }
196
197 } else { // NuDEX cannot be applied: use G4PhotonEvaporation
198
199 // Code taken from G4NeutronRadCapture
200
201 G4Fragment* aFragment = new G4Fragment( A, Z, lab4mom );
202 G4FragmentVector* fv = photonEvaporation->BreakUpFragment( aFragment );
203 if ( fv == nullptr ) fv = new G4FragmentVector;
204 fv->push_back( aFragment );
205 size_t n = fv->size();
206 for ( size_t i = 0; i < n; ++i ) {
207 G4Fragment* f = (*fv)[i];
208 G4double etot = f->GetMomentum().e();
209 Z = f->GetZ_asInt();
210 A = f->GetA_asInt();
211 const G4ParticleDefinition* theDef = nullptr;
212 if ( Z == 0 && A == 0 ) { theDef = f->GetParticleDefinition(); }
213 else if ( Z == 1 && A == 2 ) { theDef = G4Deuteron::Definition(); }
214 else if ( Z == 1 && A == 3 ) { theDef = G4Triton::Definition(); }
215 else if ( Z == 2 && A == 3 ) { theDef = G4He3::Definition(); }
216 else if ( Z == 2 && A == 4 ) { theDef = G4Alpha::Definition(); }
217 else {
218 G4double eexc = f->GetExcitationEnergy();
219 if ( eexc <= minExcitation ) eexc = 0.0;
220 theDef = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon( Z, A, eexc, noFloat, 0 );
221 }
222 G4double ekin = std::max( 0.0, etot - theDef->GetPDGMass() );
223 G4HadSecondary* news = new G4HadSecondary( new G4DynamicParticle( theDef, f->GetMomentum().vect().unit(), ekin ) );
224 G4double timeF = f->GetCreationTime();
225 if ( timeF < 0.0 ) timeF = 0.0;
226 news->SetTime( time + timeF );
227 news->SetCreatorModelID( secID );
228 theParticleChange.AddSecondary( *news );
229 delete news;
230 delete f;
231 }
232 delete fv;
233 }
234
235 return &theParticleChange;
236 }
237
238
239 G4int G4NuDEXNeutronCaptureModel::GenerateNeutronCaptureCascade( G4int theZ, G4int theA, G4double NeutronEnergy, G4int InitialLevel,
240 std::vector< char >& pType, std::vector< G4double >& pEnergy,
241 std::vector< G4double >& pTime ) {
242 // Returns the number of emitted particles. Returns -1 if the nucleus is not in the database.
243 G4int theZA = 1000*theZ + theA;
244 G4int check = Init( theZA );
245 if ( check < 0 ) return -1;
246 G4double Sn, I0;
247 theStatisticalNucleus[theZA]->GetSnAndI0( Sn, I0 ); Sn *= MeV; // I0 is the spin of the A-1 nucleus in the g.s.
248 G4double ExcitationEnergy = Sn + (theA-1.0)/(G4double)theA*NeutronEnergy;
249 G4int nPar = theStatisticalNucleus[theZA]->GenerateCascade( InitialLevel, ExcitationEnergy/MeV, pType, pEnergy, pTime );
250 for ( G4int i = 0; i < nPar; i++ ) {
251 pEnergy.at(i) *= MeV;
252 pTime.at(i) *= s;
253 }
254 return nPar;
255 }
256
257
258 G4int G4NuDEXNeutronCaptureModel::Init( G4int theZA, unsigned int seed1, unsigned int seed2, unsigned int seed3 ) {
259 if ( HasData[theZA] == -1 ) return -1;
260 if ( HasData[theZA] == 1 ) return 0;
261 if ( theStatisticalNucleus[theZA] == 0 ) {
262 G4int theZ = theZA/1000;
263 G4int theA = theZA-1000*theZ;
264 theStatisticalNucleus[theZA] = new G4NuDEXStatisticalNucleus( theZ, theA );
265 if ( BandWidth != 0 ) theStatisticalNucleus[theZA]->SetBandWidth( BandWidth );
266 theStatisticalNucleus[theZA]->SetBrOption( BrOption );
267 if ( seed1 > 0 ) theStatisticalNucleus[theZA]->SetRandom1Seed( seed1 );
268 if ( seed2 > 0 ) theStatisticalNucleus[theZA]->SetRandom1Seed( seed2 );
269 if ( seed3 > 0 ) theStatisticalNucleus[theZA]->SetRandom1Seed( seed3 );
270 G4int check = theStatisticalNucleus[theZA]->Init( NuDEXLibDirectory.c_str() );
271 if ( check < 0 ) {
272 HasData[theZA] = -1;
273 return -1;
274 } else {
275 HasData[theZA] = 1;
276 }
277 }
278 return 0;
279 }
280
281
282 G4int G4NuDEXNeutronCaptureModel::SelectInitialLevel( G4int theCompoundZ, G4int theCompoundA, G4double NeutronEnergy, G4int lspin, G4int jspinx2 ) {
283 // Initial level for neutron capture. If jspinx2 < 0 it is sampled according to the 2J+1 rule.
284 // l-spin = 0, 1, 2 --> s-wave, p-wave, d-wave ...
285 G4int theZ = theCompoundZ;
286 G4int theA = theCompoundA;
287 G4int theZA = 1000*theZ + theA;
288 G4int check = Init( theZA );
289 if ( check < 0 ) return -1;
290 G4double Sn, I0;
291 theStatisticalNucleus[theZA]->GetSnAndI0( Sn, I0 ); Sn *= MeV; // I0 is the spin of the A-1 nucleus in the g.s.
292 G4double ExcitationEnergy = Sn + (theA-1.0)/(G4double)theA*NeutronEnergy;
293 if ( lspin < 0 ) lspin = 0;
294 if ( jspinx2 < 0 ) jspinx2 = SampleJ( theZ, theA, lspin );
295 G4bool parity = false;
296 if ( ( I0 >= 0 && (lspin%2) == 0 ) || ( I0 < 0 && (lspin%2) == 1 ) ) parity = true;
297 G4int InitialLevel = theStatisticalNucleus[theZA]->GetClosestLevel( ExcitationEnergy/MeV, jspinx2, parity );
298 return InitialLevel;
299 }
300
301
302 G4int G4NuDEXNeutronCaptureModel ::SampleJ( G4int theCompoundZ, G4int theCompoundA, G4int lspin ) {
303 // Samples J for this l-spin (l-spin = 0, 1, 2 --> s-wave, p-wave, d-wave ...)
304 // The probability will be proportional to 2J+1
305 // Returns Jx2
306 G4int AllowedJx2[100];
307 G4int NAllowedJvals = GetAllowedJx2values( theCompoundZ, theCompoundA, lspin, AllowedJx2 );
308 G4double AllowedJx2CumulProb[100], TotalCumul = 0.0;
309 for ( G4int i = 0; i < NAllowedJvals; i++ ) {
310 AllowedJx2CumulProb[i] = AllowedJx2[i] + 1.0;
311 TotalCumul += AllowedJx2CumulProb[i];
312 }
313 for ( G4int i = 0; i < NAllowedJvals; i++ ) {
314 AllowedJx2CumulProb[i] /= TotalCumul;
315 if ( i > 0 ) AllowedJx2CumulProb[i] += AllowedJx2CumulProb[i-1];
316 }
317 G4double rand = G4UniformRand();
318 G4int i_result = -1;
319 for ( G4int i = 0; i < NAllowedJvals; i++ ) {
320 if ( rand < AllowedJx2CumulProb[i] ) {
321 i_result = i; break;
322 }
323 }
324 if ( i_result < 0 ) {
325 G4cerr << " ############ Error in " << __FILE__ << ", line " << __LINE__ << " ############"<< G4endl;
326 exit(1);
327 }
328 G4int jspinx2 = AllowedJx2[i_result];
329 return jspinx2;
330 }
331
332
333 G4int G4NuDEXNeutronCaptureModel::GetAllowedJx2values( G4int theCompoundZ, G4int theCompoundA, G4int lspin, G4int* jx2vals ) {
334 // Provides the allowed jx2 values in neutron capture for a certain l-spin (l-spin = 0, 1, 2 --> s-wave, p-wave, d-wave ...)
335 G4int theZA = 1000*theCompoundZ + theCompoundA;
336 G4int check = Init( theZA );
337 if ( check < 0 ) return -1;
338 G4double Sn, I0;
339 theStatisticalNucleus[theZA]->GetSnAndI0( Sn, I0 ); Sn *= MeV; // I0 is the spin of the A-1 nucleus in the g.s.
340 G4int Ix2 = (G4int)( ( std::fabs(I0) + 0.1 )*2.0 );
341 G4int Jx2min = std::min( std::abs( Ix2-1-2*lspin ), std::abs( Ix2+1-2*lspin ) );
342 G4int Jx2max = Ix2 + 1 + 2*lspin;
343 G4int NAllowedJvals = 0;
344 for ( G4int Jx2 = Jx2min; Jx2 <= Jx2max; Jx2 += 2 ) {
345 if ( Jx2 >= 0 ) {
346 jx2vals[NAllowedJvals] = Jx2;
347 NAllowedJvals++;
348 }
349 }
350 return NAllowedJvals;
351 }
352