| 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 //
28 // File: G4LENDInelastic.cc
29 // Date: 24 March 2020
30 // Author: Dennis Wright
31 //
32 // Description: model for inelastic scatterin
33 // gammas at energies of order 2
34 // This model uses GIDI particle
35 // in spectrum mode. In this mo
36 // for each possible particle ty
37 // given interaction. Unlike Ge
38 // consideration of event-by-eve
39 // Indeed, forcing such conserva
40 // introduces correlations and d
41 // spectra which are not present
42 //
43 // In order to use GIDI data wit
44 // minimal event-by-event baryon
45 // enforced which allows deviati
46 // giving warnings. Neither cha
47 // conservation is enforced. Un
48 // fragment (n, p, d, t, alpha)
49 // after a large number of event
50 // conservation also approach th
51 // this limit. The mass, charge
52 // large fragments, however, are
53 // data very well. This is a re
54 // baryon number conservation an
55 // fragment spectra are correct.
56 //
57 //////////////////////////////////////////////
58
59 #include "G4LENDInelastic.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "G4Nucleus.hh"
62 #include "G4IonTable.hh"
63 #include <algorithm>
64 #include <random>
65
66 G4HadFinalState* G4LENDInelastic::ApplyYoursel
67
68 {
69 G4ThreeVector projMom = aTrack.Get4Momentum(
70 G4double temp = aTrack.GetMaterial()->GetTem
71
72 G4int iZ = aTarg.GetZ_asInt();
73 G4int iA = aTarg.GetA_asInt();
74 G4int iM = 0;
75 if (aTarg.GetIsotope() != nullptr) iM = aTar
76
77 G4double ke = aTrack.GetKineticEnergy();
78
79 G4HadFinalState* theResult = &theParticleCha
80 theResult->Clear();
81
82 G4GIDI_target* aGIDITarget =
83 get_target_from_map(lend_manager->GetNucle
84 if (aGIDITarget == nullptr) {
85 // G4cout << " No target found " << G4endl;
86 theParticleChange.Clear();
87 theParticleChange.SetStatusChange(isAlive)
88 theParticleChange.SetEnergyChange(aTrack.G
89 theParticleChange.SetMomentumChange(aTrack
90 return &theParticleChange;
91 }
92 // return returnUnchanged(aTrack, theResult);
93
94 // Get GIDI final state products for givern
95 G4int loop(0);
96 G4int loopMax = 1000;
97 std::vector<G4GIDI_Product>* products;
98 do {
99 products = aGIDITarget->getOthersFinalStat
100 loop++;
101 } while (products == nullptr && loop < loopM
102
103 // G4LENDInelastic accepts all light fragmen
104 // and removes any heavy fragments which cau
105 // Charge and energy non-conservation still
106 // of events, this improves on average.
107
108 if (loop > loopMax - 1) {
109 // G4cout << " too many loops, return initi
110
111 theParticleChange.Clear();
112 theParticleChange.SetStatusChange(isAlive)
113 theParticleChange.SetEnergyChange(aTrack.G
114 theParticleChange.SetMomentumChange(aTrack
115 return &theParticleChange;
116
117 // if (aTrack.GetDefinition() == G4Proton::
118 // aTrack.GetDefinition() == G4Neutron:
119 // theResult = preco->ApplyYourself(aTra
120 // } else {
121 // theResult = returnUnchanged(aTrack, t
122 // }
123
124 } else {
125 G4int iTotZ = iZ + aTrack.GetDefinition()-
126 G4int iTotA = iA + aTrack.GetDefinition()-
127
128 // Loop over GIDI products and separate li
129 G4int GZtot(0);
130 G4int GAtot(0);
131 G4int productA(0);
132 G4int productZ(0);
133 std::vector<G4int> lightProductIndex;
134 std::vector<G4int> heavyProductIndex;
135 for (G4int i = 0; i < int( products->size(
136 productA = (*products)[i].A;
137 if (productA < 5) {
138 lightProductIndex.push_back(i);
139 GZtot += (*products)[i].Z;
140 GAtot += productA;
141 } else {
142 heavyProductIndex.push_back(i);
143 }
144 }
145
146 // Randomize order of heavies to correct s
147 // std::random_shuffle(heavyProductIndex.b
148 // std::cout << " Heavy product index befo
149 // for (G4int i = 0; i < int(heavyProductI
150 // std::cout << std::endl;
151
152 auto rng = std::default_random_engine {};
153 std::shuffle(heavyProductIndex.begin(), he
154
155 // std::cout << " Heavy product index afte
156 // for (G4int i = 0; i < int(heavyProductI
157 // std::cout << std::endl;
158
159 std::vector<G4int> savedHeavyIndex;
160 G4int itest(0);
161 for (G4int i = 0; i < int(heavyProductInde
162 itest = heavyProductIndex[i];
163 productA = (*products)[itest].A;
164 productZ = (*products)[itest].Z;
165 if ((GAtot + productA <= iTotA) && (GZto
166 savedHeavyIndex.push_back(itest);
167 GZtot += productZ;
168 GAtot += productA;
169 }
170 }
171
172 /*
173 G4cout << " saved light products = ";
174 for (G4int k = 0; k < int(lightProductInde
175 itest = lightProductIndex[k];
176 G4cout << "(" << (*products)[itest].Z <<
177 }
178 G4cout << G4endl;
179
180 G4cout << " saved heavy products = ";
181 for (G4int k = 0; k < int(savedHeavyIndex.
182 itest = savedHeavyIndex[k];
183 G4cout << "(" << (*products)[itest].Z <<
184 }
185 G4cout << G4endl;
186 */
187 // Now convert saved products to Geant4 pa
188 // Note that, at least for heavy fragments
189 // have slightly different values.
190
191 G4DynamicParticle* theSec = nullptr;
192 G4ThreeVector Psum;
193 for (G4int i = 0; i < int(lightProductInde
194 itest = lightProductIndex[i];
195 productZ = (*products)[itest].Z;
196 productA = (*products)[itest].A;
197 theSec = new G4DynamicParticle();
198 if (productA == 1 && productZ == 0) {
199 theSec->SetDefinition(G4Neutron::Neutr
200 } else if (productA == 1 && productZ ==
201 theSec->SetDefinition(G4Proton::Proton
202 } else if (productA == 2 && productZ ==
203 theSec->SetDefinition(G4Deuteron::Deut
204 } else if (productA == 3 && productZ ==
205 theSec->SetDefinition(G4Triton::Triton
206 } else if (productA == 4 && productZ ==
207 theSec->SetDefinition(G4Alpha::Alpha()
208 } else {
209 theSec->SetDefinition(G4Gamma::Gamma()
210 }
211
212 G4ThreeVector momentum((*products)[itest
213 (*products)[itest
214 (*products)[itest
215 Psum += momentum;
216 theSec->SetMomentum(momentum);
217 // theResult->AddSecondary(theSec);
218 theParticleChange.AddSecondary(theSec, s
219 }
220
221 G4int productM(0);
222 for (G4int i = 0; i < int(savedHeavyIndex.
223 itest = savedHeavyIndex[i];
224 productZ = (*products)[itest].Z;
225 productA = (*products)[itest].A;
226 productM = (*products)[itest].m;
227 theSec = new G4DynamicParticle();
228 theSec->SetDefinition(G4IonTable::GetIon
229
230
231 G4ThreeVector momentum((*products)[itest
232 (*products)[itest
233 (*products)[itest
234 Psum += momentum;
235 theSec->SetMomentum(momentum);
236 // theResult->AddSecondary(theSec);
237 theParticleChange.AddSecondary(theSec, s
238 }
239
240 // Create heavy fragment if necessary to t
241 // Note: this step is only required to pre
242 // where "catastrophic" non-conserva
243 // The residual generated will not n
244 // occur in the actual reaction.
245 if (iTotA - GAtot > 1) {
246 theSec = new G4DynamicParticle();
247 if (iTotZ == GZtot) {
248 // Special case when a nucleus of only
249 // Violate charge conservation and set
250 // G4cout << " Z = 1, A = "<< iTotA -
251 theSec->SetDefinition(G4IonTable::GetI
252 } else {
253 theSec->SetDefinition(G4IonTable::GetI
254 }
255 theSec->SetMomentum(projMom - Psum);
256 // theResult->AddSecondary(theSec);
257 theParticleChange.AddSecondary(theSec, s
258 }
259 } // loop OK
260
261 delete products;
262 // theResult->SetStatusChange( stopAndKill );
263 theParticleChange.SetStatusChange( stopAndKi
264 // return theResult;
265 return &theParticleChange;
266 }
267