| 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 // Physics model class G4NeutronRadCapture
29 // Created: 31 August 2009
30 // Author V.Ivanchenko
31 //
32 // Modified:
33 // 09.09.2010 V.Ivanchenko added usage of G4PhotonEvaporation
34 //
35
36 #include "G4NeutronRadCapture.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4ParticleDefinition.hh"
39 #include "G4Fragment.hh"
40 #include "G4FragmentVector.hh"
41 #include "G4NucleiProperties.hh"
42 #include "G4VEvaporationChannel.hh"
43 #include "G4PhotonEvaporation.hh"
44 #include "G4DynamicParticle.hh"
45 #include "G4ParticleTable.hh"
46 #include "G4IonTable.hh"
47 #include "G4Electron.hh"
48 #include "G4Deuteron.hh"
49 #include "G4Triton.hh"
50 #include "G4He3.hh"
51 #include "G4Alpha.hh"
52 #include "G4RandomDirection.hh"
53 #include "G4HadronicParameters.hh"
54 #include "G4PhysicsModelCatalog.hh"
55
56 G4NeutronRadCapture::G4NeutronRadCapture()
57 : G4HadronicInteraction("nRadCapture"),
58 photonEvaporation(nullptr),lab4mom(0.,0.,0.,0.)
59 {
60 lowestEnergyLimit = 10*CLHEP::eV;
61 minExcitation = 0.1*CLHEP::keV;
62 secID = -1;
63 theTableOfIons = G4ParticleTable::GetParticleTable()->GetIonTable();
64 }
65
66 G4NeutronRadCapture::~G4NeutronRadCapture()
67 {
68 delete photonEvaporation;
69 }
70
71 void G4NeutronRadCapture::InitialiseModel()
72 {
73 if(photonEvaporation != nullptr) { return; }
74 G4DeexPrecoParameters* param =
75 G4NuclearLevelData::GetInstance()->GetParameters();
76 minExcitation = param->GetMinExcitation();
77 secID = G4PhysicsModelCatalog::GetModelID("model_" + GetModelName());
78 photonEvaporation = new G4PhotonEvaporation();
79 photonEvaporation->Initialise();
80 photonEvaporation->SetICM(true);
81 }
82
83 G4HadFinalState* G4NeutronRadCapture::ApplyYourself(
84 const G4HadProjectile& aTrack, G4Nucleus& theNucleus)
85 {
86 theParticleChange.Clear();
87 theParticleChange.SetStatusChange(stopAndKill);
88
89 G4int A = theNucleus.GetA_asInt();
90 G4int Z = theNucleus.GetZ_asInt();
91
92 G4double time = aTrack.GetGlobalTime();
93
94 // Create initial state
95 lab4mom.set(0.,0.,0.,G4NucleiProperties::GetNuclearMass(A, Z));
96 lab4mom += aTrack.Get4Momentum();
97
98 G4double M = lab4mom.mag();
99 ++A;
100 G4double mass = G4NucleiProperties::GetNuclearMass(A, Z);
101 //G4cout << "Capture start: Z= " << Z << " A= " << A
102 // << " LabM= " << M << " Mcompound= " << mass << G4endl;
103
104 // simplified method of 1 gamma emission
105 if(A <= 4) {
106
107 G4ThreeVector bst = lab4mom.boostVector();
108
109 if(M - mass <= lowestEnergyLimit) {
110 return &theParticleChange;
111 }
112
113 if (verboseLevel > 1) {
114 G4cout << "G4NeutronRadCapture::DoIt: Eini(MeV)="
115 << aTrack.GetKineticEnergy()/MeV << " Eexc(MeV)= "
116 << (M - mass)/MeV
117 << " Z= " << Z << " A= " << A << G4endl;
118 }
119 G4double e1 = (M - mass)*(M + mass)/(2*M);
120 G4LorentzVector lv2(e1*G4RandomDirection(),e1);
121 lv2.boost(bst);
122 G4HadSecondary* news =
123 new G4HadSecondary(new G4DynamicParticle(G4Gamma::Gamma(), lv2));
124 news->SetTime(time);
125 news->SetCreatorModelID(secID);
126 theParticleChange.AddSecondary(*news);
127 delete news;
128
129 const G4ParticleDefinition* theDef = 0;
130
131 lab4mom -= lv2;
132 if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
133 else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
134 else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
135 else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
136 else { theDef = theTableOfIons->GetIon(Z,A,0.0,noFloat,0); }
137
138 if (verboseLevel > 1) {
139 G4cout << "Gamma 4-mom: " << lv2 << " "
140 << theDef->GetParticleName() << " " << lab4mom << G4endl;
141 }
142 if(theDef) {
143 news = new G4HadSecondary(new G4DynamicParticle(theDef, lab4mom));
144 news->SetTime(time);
145 news->SetCreatorModelID(secID);
146 theParticleChange.AddSecondary(*news);
147 delete news;
148 }
149
150 // Use photon evaporation
151 } else {
152
153 // protection against wrong kinematic
154 if(M < mass) {
155 G4double etot = std::max(mass, lab4mom.e());
156 G4double ptot = std::sqrt((etot - mass)*(etot + mass));
157 G4ThreeVector v = lab4mom.vect().unit();
158 lab4mom.set(v.x()*ptot,v.y()*ptot,v.z()*ptot,etot);
159 }
160
161 G4Fragment* aFragment = new G4Fragment(A, Z, lab4mom);
162
163 if (verboseLevel > 1) {
164 G4cout << "G4NeutronRadCapture::ApplyYourself initial G4Fragmet:"
165 << G4endl;
166 G4cout << aFragment << G4endl;
167 }
168
169 //
170 // Sample final state
171 //
172 G4FragmentVector* fv = photonEvaporation->BreakUpFragment(aFragment);
173 if (nullptr == fv) { fv = new G4FragmentVector(); }
174 fv->push_back(aFragment);
175 std::size_t n = fv->size();
176
177 if (verboseLevel > 1) {
178 G4cout << "G4NeutronRadCapture: " << n << " final particles" << G4endl;
179 }
180 for(std::size_t i=0; i<n; ++i) {
181
182 G4Fragment* f = (*fv)[i];
183 G4double etot = f->GetMomentum().e();
184
185 Z = f->GetZ_asInt();
186 A = f->GetA_asInt();
187
188 const G4ParticleDefinition* theDef;
189 if(0 == Z && 0 == A) {theDef = f->GetParticleDefinition();}
190 else if (Z == 1 && A == 2) {theDef = G4Deuteron::Deuteron();}
191 else if (Z == 1 && A == 3) {theDef = G4Triton::Triton();}
192 else if (Z == 2 && A == 3) {theDef = G4He3::He3();}
193 else if (Z == 2 && A == 4) {theDef = G4Alpha::Alpha();}
194 else {
195 G4double eexc = f->GetExcitationEnergy();
196 if(eexc <= minExcitation) { eexc = 0.0; }
197 theDef = theTableOfIons->GetIon(Z, A, eexc, noFloat, 0);
198 /*
199 G4cout << "### NC Find ion Z= " << Z << " A= " << A
200 << " Eexc(MeV)= " << eexc/MeV << " "
201 << theDef << G4endl;
202 */
203 }
204 G4double ekin = std::max(0.0,etot - theDef->GetPDGMass());
205 if (verboseLevel > 1) {
206 G4cout << i << ". " << theDef->GetParticleName()
207 << " Ekin(MeV)= " << etot/MeV
208 << " p: " << f->GetMomentum().vect()
209 << G4endl;
210 }
211 G4HadSecondary* news = new G4HadSecondary(
212 new G4DynamicParticle(theDef,
213 f->GetMomentum().vect().unit(),
214 ekin));
215 G4double timeF = f->GetCreationTime();
216 if(timeF < 0.0) { timeF = 0.0; }
217 news->SetTime(time + timeF);
218 news->SetCreatorModelID(secID);
219 theParticleChange.AddSecondary(*news);
220 delete news;
221 delete f;
222 }
223 delete fv;
224 }
225 //G4cout << "Capture done" << G4endl;
226 return &theParticleChange;
227 }
228
229