| 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 // particle_hp -- source file 26 // particle_hp -- source file
27 // J.P. Wellisch, Nov-1996 27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans 28 // A prototype of the low energy neutron transport model.
29 // 29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. How 30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 070606 bug fix and migrate to enable to Par 31 // 070606 bug fix and migrate to enable to Partial cases by T. Koi
32 // 080603 bug fix for Hadron Hyper News #932 b 32 // 080603 bug fix for Hadron Hyper News #932 by T. Koi
33 // 080612 bug fix contribution from Benoit Pir 33 // 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
34 // 080717 bug fix of calculation of residual m 34 // 080717 bug fix of calculation of residual momentum by T. Koi
35 // 080801 protect negative available energy by 35 // 080801 protect negative available energy by T. Koi
36 // introduce theNDLDataA,Z which has A 36 // introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
37 // 081024 G4NucleiPropertiesTable:: to G4Nucle 37 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
38 // 090514 Fix bug in IC electron emission case 38 // 090514 Fix bug in IC electron emission case
39 // Contribution from Chao Zhang (Chao.Z 39 // Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
40 // 100406 "nothingWasKnownOnHadron=1" then sam 40 // 100406 "nothingWasKnownOnHadron=1" then sample mu isotropic in CM
41 // add two_body_reaction 41 // add two_body_reaction
42 // 100909 add safty 42 // 100909 add safty
43 // 101111 add safty for _nat_ data case in Bin 43 // 101111 add safty for _nat_ data case in Binary reaction, but break conservation
44 // 110430 add Reaction Q value and break up fl 44 // 110430 add Reaction Q value and break up flag (MF3::QI and LR)
45 // 45 //
46 // P. Arce, June-2014 Conversion neutron_hp to 46 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
47 // June-2019 - E. Mendoza - re-build "two_body 47 // June-2019 - E. Mendoza - re-build "two_body_reaction", to be used by incident charged particles
48 // (now isotropic emission in the CMS). Also r 48 // (now isotropic emission in the CMS). Also restrict nresp use below 20 MeV (for future
49 // developments). Add photon emission when no 49 // developments). Add photon emission when no data available.
50 // 50 //
51 // nresp71_m03.hh and nresp71_m02.hh are alike 51 // nresp71_m03.hh and nresp71_m02.hh are alike. The only difference between m02 and m03
52 // is in the total carbon cross section that i 52 // is in the total carbon cross section that is properly included in the latter.
53 // These data are not used in nresp71_m0*.hh. 53 // These data are not used in nresp71_m0*.hh.
54 // 54 //
55 55
56 #include "G4ParticleHPInelasticCompFS.hh" 56 #include "G4ParticleHPInelasticCompFS.hh"
57 57
58 #include "G4Alpha.hh" 58 #include "G4Alpha.hh"
59 #include "G4Electron.hh" 59 #include "G4Electron.hh"
60 #include "G4He3.hh" 60 #include "G4He3.hh"
61 #include "G4IonTable.hh" 61 #include "G4IonTable.hh"
62 #include "G4NRESP71M03.hh" 62 #include "G4NRESP71M03.hh"
63 #include "G4NucleiProperties.hh" 63 #include "G4NucleiProperties.hh"
64 #include "G4Nucleus.hh" 64 #include "G4Nucleus.hh"
65 #include "G4ParticleHPDataUsed.hh" 65 #include "G4ParticleHPDataUsed.hh"
66 #include "G4ParticleHPManager.hh" 66 #include "G4ParticleHPManager.hh"
67 #include "G4RandomDirection.hh" 67 #include "G4RandomDirection.hh"
68 #include "G4SystemOfUnits.hh" 68 #include "G4SystemOfUnits.hh"
69 69
70 #include <fstream> 70 #include <fstream>
71 71
72 G4ParticleHPInelasticCompFS::G4ParticleHPInela 72 G4ParticleHPInelasticCompFS::G4ParticleHPInelasticCompFS()
73 : G4ParticleHPFinalState() 73 : G4ParticleHPFinalState()
74 { 74 {
75 QI.resize(51); 75 QI.resize(51);
76 LR.resize(51); 76 LR.resize(51);
77 for (G4int i = 0; i < 51; ++i) { 77 for (G4int i = 0; i < 51; ++i) {
78 hasXsec = true; 78 hasXsec = true;
79 theXsection[i] = nullptr; 79 theXsection[i] = nullptr;
80 theEnergyDistribution[i] = nullptr; 80 theEnergyDistribution[i] = nullptr;
81 theAngularDistribution[i] = nullptr; 81 theAngularDistribution[i] = nullptr;
82 theEnergyAngData[i] = nullptr; 82 theEnergyAngData[i] = nullptr;
83 theFinalStatePhotons[i] = nullptr; 83 theFinalStatePhotons[i] = nullptr;
84 QI[i] = 0.0; 84 QI[i] = 0.0;
85 LR[i] = 0; 85 LR[i] = 0;
86 } 86 }
87 } 87 }
88 88
89 G4ParticleHPInelasticCompFS::~G4ParticleHPInel 89 G4ParticleHPInelasticCompFS::~G4ParticleHPInelasticCompFS()
90 { 90 {
91 for (G4int i = 0; i < 51; ++i) { 91 for (G4int i = 0; i < 51; ++i) {
92 if (theXsection[i] != nullptr) delete theX 92 if (theXsection[i] != nullptr) delete theXsection[i];
93 if (theEnergyDistribution[i] != nullptr) d 93 if (theEnergyDistribution[i] != nullptr) delete theEnergyDistribution[i];
94 if (theAngularDistribution[i] != nullptr) 94 if (theAngularDistribution[i] != nullptr) delete theAngularDistribution[i];
95 if (theEnergyAngData[i] != nullptr) delete 95 if (theEnergyAngData[i] != nullptr) delete theEnergyAngData[i];
96 if (theFinalStatePhotons[i] != nullptr) de 96 if (theFinalStatePhotons[i] != nullptr) delete theFinalStatePhotons[i];
97 } 97 }
98 } 98 }
99 99
100 void G4ParticleHPInelasticCompFS::InitDistribu 100 void G4ParticleHPInelasticCompFS::InitDistributionInitialState(G4ReactionProduct& inPart,
101 G4ReactionPr 101 G4ReactionProduct& aTarget, G4int it)
102 { 102 {
103 if (theAngularDistribution[it] != nullptr) { 103 if (theAngularDistribution[it] != nullptr) {
104 theAngularDistribution[it]->SetTarget(aTar 104 theAngularDistribution[it]->SetTarget(aTarget);
105 theAngularDistribution[it]->SetProjectileR 105 theAngularDistribution[it]->SetProjectileRP(inPart);
106 } 106 }
107 107
108 if (theEnergyAngData[it] != nullptr) { 108 if (theEnergyAngData[it] != nullptr) {
109 theEnergyAngData[it]->SetTarget(aTarget); 109 theEnergyAngData[it]->SetTarget(aTarget);
110 theEnergyAngData[it]->SetProjectileRP(inPa 110 theEnergyAngData[it]->SetProjectileRP(inPart);
111 } 111 }
112 } 112 }
113 113
114 void G4ParticleHPInelasticCompFS::InitGammas(G 114 void G4ParticleHPInelasticCompFS::InitGammas(G4double AR, G4double ZR)
115 { 115 {
116 G4int Z = G4lrint(ZR); 116 G4int Z = G4lrint(ZR);
117 G4int A = G4lrint(AR); 117 G4int A = G4lrint(AR);
118 std::ostringstream ost; 118 std::ostringstream ost;
119 ost << gammaPath << "z" << Z << ".a" << A; 119 ost << gammaPath << "z" << Z << ".a" << A;
120 G4String aName = ost.str(); 120 G4String aName = ost.str();
121 std::ifstream from(aName, std::ios::in); 121 std::ifstream from(aName, std::ios::in);
122 122
123 if (!from) return; // no data found for thi 123 if (!from) return; // no data found for this isotope
124 std::ifstream theGammaData(aName, std::ios:: 124 std::ifstream theGammaData(aName, std::ios::in);
125 125
126 theGammas.Init(theGammaData); 126 theGammas.Init(theGammaData);
127 } 127 }
128 128
129 void G4ParticleHPInelasticCompFS::Init(G4doubl << 129 void G4ParticleHPInelasticCompFS::Init(G4double A, G4double Z, G4int M, G4String& dirName,
130 const G << 130 G4String& aFSType, G4ParticleDefinition*)
131 { 131 {
132 gammaPath = fManager->GetNeutronHPPath() + " 132 gammaPath = fManager->GetNeutronHPPath() + "/Inelastic/Gammas/";
133 const G4String& tString = dirName; << 133 G4String tString = dirName;
134 SetA_Z(A, Z, M); 134 SetA_Z(A, Z, M);
135 G4bool dbool; 135 G4bool dbool;
136 const G4ParticleHPDataUsed& aFile = << 136 G4ParticleHPDataUsed aFile =
137 theNames.GetName(theBaseA, theBaseZ, M, tS 137 theNames.GetName(theBaseA, theBaseZ, M, tString, aFSType, dbool);
138 SetAZMs(aFile); 138 SetAZMs(aFile);
139 const G4String& filename = aFile.GetName(); << 139 G4String filename = aFile.GetName();
140 #ifdef G4VERBOSE 140 #ifdef G4VERBOSE
141 if (fManager->GetDEBUG()) 141 if (fManager->GetDEBUG())
142 G4cout << " G4ParticleHPInelasticCompFS::I 142 G4cout << " G4ParticleHPInelasticCompFS::Init FILE " << filename << G4endl;
143 #endif 143 #endif
144 144
145 SetAZMs(A, Z, M, aFile); 145 SetAZMs(A, Z, M, aFile);
146 146
147 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ 147 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ != theBaseZ || theNDLDataA != theBaseA)))
148 { 148 {
149 #ifdef G4VERBOSE 149 #ifdef G4VERBOSE
150 if (fManager->GetDEBUG()) 150 if (fManager->GetDEBUG())
151 G4cout << "Skipped = " << filename << " 151 G4cout << "Skipped = " << filename << " " << A << " " << Z << G4endl;
152 #endif 152 #endif
153 hasAnyData = false; 153 hasAnyData = false;
154 hasFSData = false; 154 hasFSData = false;
155 hasXsec = false; 155 hasXsec = false;
156 return; 156 return;
157 } 157 }
158 std::istringstream theData(std::ios::in); 158 std::istringstream theData(std::ios::in);
159 fManager->GetDataStream(filename, theData); 159 fManager->GetDataStream(filename, theData);
160 if (!theData) //"!" is a operator of ios 160 if (!theData) //"!" is a operator of ios
161 { 161 {
162 hasAnyData = false; 162 hasAnyData = false;
163 hasFSData = false; 163 hasFSData = false;
164 hasXsec = false; 164 hasXsec = false;
165 return; 165 return;
166 } 166 }
167 // here we go 167 // here we go
168 G4int infoType, dataType, dummy; 168 G4int infoType, dataType, dummy;
169 G4int sfType, it; 169 G4int sfType, it;
170 hasFSData = false; 170 hasFSData = false;
171 while (theData >> infoType) // Loop checkin 171 while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
172 { 172 {
173 hasFSData = true; 173 hasFSData = true;
174 theData >> dataType; 174 theData >> dataType;
175 theData >> sfType >> dummy; 175 theData >> sfType >> dummy;
176 it = 50; 176 it = 50;
177 if (sfType >= 600 || (sfType < 100 && sfTy 177 if (sfType >= 600 || (sfType < 100 && sfType >= 50))
178 it = sfType % 50; 178 it = sfType % 50;
179 if (dataType == 3) 179 if (dataType == 3)
180 { 180 {
181 G4double dqi; 181 G4double dqi;
182 G4int ilr; 182 G4int ilr;
183 theData >> dqi >> ilr; 183 theData >> dqi >> ilr;
184 184
185 QI[it] = dqi * CLHEP::eV; 185 QI[it] = dqi * CLHEP::eV;
186 LR[it] = ilr; 186 LR[it] = ilr;
187 theXsection[it] = new G4ParticleHPVector 187 theXsection[it] = new G4ParticleHPVector;
188 G4int total; 188 G4int total;
189 theData >> total; 189 theData >> total;
190 theXsection[it]->Init(theData, total, CL 190 theXsection[it]->Init(theData, total, CLHEP::eV);
191 } 191 }
192 else if (dataType == 4) { 192 else if (dataType == 4) {
193 theAngularDistribution[it] = new G4Parti 193 theAngularDistribution[it] = new G4ParticleHPAngular;
194 theAngularDistribution[it]->Init(theData 194 theAngularDistribution[it]->Init(theData);
195 } 195 }
196 else if (dataType == 5) { 196 else if (dataType == 5) {
197 theEnergyDistribution[it] = new G4Partic 197 theEnergyDistribution[it] = new G4ParticleHPEnergyDistribution;
198 theEnergyDistribution[it]->Init(theData) 198 theEnergyDistribution[it]->Init(theData);
199 } 199 }
200 else if (dataType == 6) { 200 else if (dataType == 6) {
201 theEnergyAngData[it] = new G4ParticleHPE 201 theEnergyAngData[it] = new G4ParticleHPEnAngCorrelation(theProjectile);
202 // G4cout << this << " CompFS theEn 202 // G4cout << this << " CompFS theEnergyAngData " << it << theEnergyAngData[it] << G4endl;
203 theEnergyAngData[it]->Init(theData); 203 theEnergyAngData[it]->Init(theData);
204 } 204 }
205 else if (dataType == 12) { 205 else if (dataType == 12) {
206 theFinalStatePhotons[it] = new G4Particl 206 theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
207 theFinalStatePhotons[it]->InitMean(theDa 207 theFinalStatePhotons[it]->InitMean(theData);
208 } 208 }
209 else if (dataType == 13) { 209 else if (dataType == 13) {
210 theFinalStatePhotons[it] = new G4Particl 210 theFinalStatePhotons[it] = new G4ParticleHPPhotonDist;
211 theFinalStatePhotons[it]->InitPartials(t 211 theFinalStatePhotons[it]->InitPartials(theData, theXsection[50]);
212 } 212 }
213 else if (dataType == 14) { 213 else if (dataType == 14) {
214 theFinalStatePhotons[it]->InitAngular(th 214 theFinalStatePhotons[it]->InitAngular(theData);
215 } 215 }
216 else if (dataType == 15) { 216 else if (dataType == 15) {
217 theFinalStatePhotons[it]->InitEnergies(t 217 theFinalStatePhotons[it]->InitEnergies(theData);
218 } 218 }
219 else { 219 else {
220 G4ExceptionDescription ed; 220 G4ExceptionDescription ed;
221 ed << "Z=" << theBaseZ << " A=" << theBa 221 ed << "Z=" << theBaseZ << " A=" << theBaseA << " dataType=" << dataType
222 << " projectile: " << theProjectile->GetPar 222 << " projectile: " << theProjectile->GetParticleName();
223 G4Exception("G4ParticleHPInelasticCompFS 223 G4Exception("G4ParticleHPInelasticCompFS::Init", "hadr01", JustWarning,
224 ed, "Data-type unknown"); 224 ed, "Data-type unknown");
225 } 225 }
226 } 226 }
227 } 227 }
228 228
229 G4int G4ParticleHPInelasticCompFS::SelectExitC 229 G4int G4ParticleHPInelasticCompFS::SelectExitChannel(G4double eKinetic)
230 { 230 {
231 G4double running[50]; 231 G4double running[50];
232 running[0] = 0; 232 running[0] = 0;
233 G4int i; 233 G4int i;
234 for (i = 0; i < 50; ++i) { 234 for (i = 0; i < 50; ++i) {
235 if (i != 0) running[i] = running[i - 1]; 235 if (i != 0) running[i] = running[i - 1];
236 if (theXsection[i] != nullptr) { 236 if (theXsection[i] != nullptr) {
237 running[i] += std::max(0., theXsection[i 237 running[i] += std::max(0., theXsection[i]->GetXsec(eKinetic));
238 } 238 }
239 } 239 }
240 G4double random = G4UniformRand(); 240 G4double random = G4UniformRand();
241 G4double sum = running[49]; 241 G4double sum = running[49];
242 G4int it = 50; 242 G4int it = 50;
243 if (0 != sum) { 243 if (0 != sum) {
244 G4int i0; 244 G4int i0;
245 for (i0 = 0; i0 < 50; ++i0) { 245 for (i0 = 0; i0 < 50; ++i0) {
246 it = i0; 246 it = i0;
247 if (random < running[i0] / sum) break; 247 if (random < running[i0] / sum) break;
248 } 248 }
249 } 249 }
250 return it; 250 return it;
251 } 251 }
252 252
253 // n,p,d,t,he3,a 253 // n,p,d,t,he3,a
254 void G4ParticleHPInelasticCompFS::CompositeApp 254 void G4ParticleHPInelasticCompFS::CompositeApply(const G4HadProjectile& theTrack,
255 255 G4ParticleDefinition* aDefinition)
256 { 256 {
257 // prepare neutron 257 // prepare neutron
258 if (theResult.Get() == nullptr) theResult.Pu 258 if (theResult.Get() == nullptr) theResult.Put(new G4HadFinalState);
259 theResult.Get()->Clear(); 259 theResult.Get()->Clear();
260 G4double eKinetic = theTrack.GetKineticEnerg 260 G4double eKinetic = theTrack.GetKineticEnergy();
261 const G4HadProjectile* hadProjectile = &theT 261 const G4HadProjectile* hadProjectile = &theTrack;
262 G4ReactionProduct incidReactionProduct(hadPr 262 G4ReactionProduct incidReactionProduct(hadProjectile->GetDefinition());
263 incidReactionProduct.SetMomentum(hadProjecti 263 incidReactionProduct.SetMomentum(hadProjectile->Get4Momentum().vect());
264 incidReactionProduct.SetKineticEnergy(eKinet 264 incidReactionProduct.SetKineticEnergy(eKinetic);
265 265
266 // prepare target 266 // prepare target
267 for (G4int i = 0; i < 50; ++i) { 267 for (G4int i = 0; i < 50; ++i) {
268 if (theXsection[i] != nullptr) { 268 if (theXsection[i] != nullptr) {
269 break; 269 break;
270 } 270 }
271 } 271 }
272 272
273 G4double targetMass = G4NucleiProperties::Ge 273 G4double targetMass = G4NucleiProperties::GetNuclearMass(theBaseA, theBaseZ);
274 #ifdef G4VERBOSE 274 #ifdef G4VERBOSE
275 if (fManager->GetDEBUG()) 275 if (fManager->GetDEBUG())
276 G4cout << "G4ParticleHPInelasticCompFS::Co 276 G4cout << "G4ParticleHPInelasticCompFS::CompositeApply A=" << theBaseA << " Z="
277 << theBaseZ << " incident " << hadProject 277 << theBaseZ << " incident " << hadProjectile->GetDefinition()->GetParticleName()
278 << G4endl; 278 << G4endl;
279 #endif 279 #endif
280 G4ReactionProduct theTarget; 280 G4ReactionProduct theTarget;
281 G4Nucleus aNucleus; 281 G4Nucleus aNucleus;
282 // G4ThreeVector neuVelo = 282 // G4ThreeVector neuVelo =
283 // (1./hadProjectile->GetDefinition()->GetPD 283 // (1./hadProjectile->GetDefinition()->GetPDGMass())*incidReactionProduct.GetMomentum(); theTarget
284 // = aNucleus.GetBiasedThermalNucleus( targe 284 // = aNucleus.GetBiasedThermalNucleus( targetMass/hadProjectile->GetDefinition()->GetPDGMass() ,
285 // neuVelo, theTrack.GetMaterial()->GetTempe 285 // neuVelo, theTrack.GetMaterial()->GetTemperature()); G4Nucleus::GetBiasedThermalNucleus requests
286 // normalization of mass and velocity in neu 286 // normalization of mass and velocity in neutron mass
287 G4ThreeVector neuVelo = incidReactionProduct 287 G4ThreeVector neuVelo = incidReactionProduct.GetMomentum() / CLHEP::neutron_mass_c2;
288 theTarget = aNucleus.GetBiasedThermalNucleus 288 theTarget = aNucleus.GetBiasedThermalNucleus(targetMass / CLHEP::neutron_mass_c2,
289 289 neuVelo, theTrack.GetMaterial()->GetTemperature());
290 290
291 theTarget.SetDefinition(G4IonTable::GetIonTa 291 theTarget.SetDefinition(G4IonTable::GetIonTable()->GetIon(theBaseZ, theBaseA, 0.0));
292 292
293 // prepare the residual mass 293 // prepare the residual mass
294 G4double residualMass = 0; 294 G4double residualMass = 0;
295 G4int residualZ = theBaseZ + 295 G4int residualZ = theBaseZ +
296 G4lrint((theProjectile->GetPDGCharge() - a 296 G4lrint((theProjectile->GetPDGCharge() - aDefinition->GetPDGCharge())/CLHEP::eplus);
297 G4int residualA = theBaseA + theProjectile-> 297 G4int residualA = theBaseA + theProjectile->GetBaryonNumber() - aDefinition->GetBaryonNumber();
298 residualMass = G4NucleiProperties::GetNuclea 298 residualMass = G4NucleiProperties::GetNuclearMass(residualA, residualZ);
299 299
300 // prepare energy in target rest frame 300 // prepare energy in target rest frame
301 G4ReactionProduct boosted; 301 G4ReactionProduct boosted;
302 boosted.Lorentz(incidReactionProduct, theTar 302 boosted.Lorentz(incidReactionProduct, theTarget);
303 eKinetic = boosted.GetKineticEnergy(); 303 eKinetic = boosted.GetKineticEnergy();
304 304
305 // select exit channel for composite FS clas 305 // select exit channel for composite FS class.
306 G4int it = SelectExitChannel(eKinetic); 306 G4int it = SelectExitChannel(eKinetic);
307 307
308 // E. Mendoza (2018) -- to use JENDL/AN-2005 308 // E. Mendoza (2018) -- to use JENDL/AN-2005
309 if (theEnergyDistribution[it] == nullptr && 309 if (theEnergyDistribution[it] == nullptr && theAngularDistribution[it] == nullptr
310 && theEnergyAngData[it] == nullptr) 310 && theEnergyAngData[it] == nullptr)
311 { 311 {
312 if (theEnergyDistribution[50] != nullptr | 312 if (theEnergyDistribution[50] != nullptr || theAngularDistribution[50] != nullptr
313 || theEnergyAngData[50] != nullptr) 313 || theEnergyAngData[50] != nullptr)
314 { 314 {
315 it = 50; 315 it = 50;
316 } 316 }
317 } 317 }
318 318
319 // set target and neutron in the relevant ex 319 // set target and neutron in the relevant exit channel
320 InitDistributionInitialState(incidReactionPr 320 InitDistributionInitialState(incidReactionProduct, theTarget, it);
321 321
322 //------------------------------------------ 322 //---------------------------------------------------------------------//
323 // Hook for NRESP71MODEL 323 // Hook for NRESP71MODEL
324 if (fManager->GetUseNRESP71Model() && eKinet 324 if (fManager->GetUseNRESP71Model() && eKinetic < 20 * CLHEP::MeV) {
325 if (theBaseZ == 6) // If the reaction is 325 if (theBaseZ == 6) // If the reaction is with Carbon...
326 { 326 {
327 if (theProjectile == G4Neutron::Definiti 327 if (theProjectile == G4Neutron::Definition()) {
328 if (use_nresp71_model(aDefinition, it, 328 if (use_nresp71_model(aDefinition, it, theTarget, boosted)) return;
329 } 329 }
330 } 330 }
331 } 331 }
332 //------------------------------------------ 332 //---------------------------------------------------------------------//
333 333
334 G4ReactionProductVector* thePhotons = nullpt 334 G4ReactionProductVector* thePhotons = nullptr;
335 G4ReactionProductVector* theParticles = null 335 G4ReactionProductVector* theParticles = nullptr;
336 G4ReactionProduct aHadron; 336 G4ReactionProduct aHadron;
337 aHadron.SetDefinition(aDefinition); // what 337 aHadron.SetDefinition(aDefinition); // what if only cross-sections exist ==> Na 23 11 @@@@
338 G4double availableEnergy = incidReactionProd 338 G4double availableEnergy = incidReactionProduct.GetKineticEnergy()
339 + incidReactionPr 339 + incidReactionProduct.GetMass() - aHadron.GetMass()
340 + (targetMass - r 340 + (targetMass - residualMass);
341 341
342 if (availableEnergy < 0) { 342 if (availableEnergy < 0) {
343 availableEnergy = 0; 343 availableEnergy = 0;
344 } 344 }
345 G4int nothingWasKnownOnHadron = 0; 345 G4int nothingWasKnownOnHadron = 0;
346 G4double eGamm = 0; 346 G4double eGamm = 0;
347 G4int iLevel = -1; 347 G4int iLevel = -1;
348 // max gamma energy and index 348 // max gamma energy and index
349 G4int imaxEx = theGammas.GetNumberOfLevels() 349 G4int imaxEx = theGammas.GetNumberOfLevels() - 1;
350 350
351 // without photon has it = 0 351 // without photon has it = 0
352 if (50 == it) { 352 if (50 == it) {
353 // Excitation level is not determined 353 // Excitation level is not determined
354 aHadron.SetKineticEnergy(availableEnergy * 354 aHadron.SetKineticEnergy(availableEnergy * residualMass / (aHadron.GetMass() + residualMass));
355 355
356 // TK add safty 100909 356 // TK add safty 100909
357 G4double p2 = 357 G4double p2 =
358 (aHadron.GetTotalEnergy() * aHadron.GetT 358 (aHadron.GetTotalEnergy() * aHadron.GetTotalEnergy() - aHadron.GetMass() * aHadron.GetMass());
359 G4double p = (p2 > 0.0) ? std::sqrt(p2) : 359 G4double p = (p2 > 0.0) ? std::sqrt(p2) : 0.0;
360 aHadron.SetMomentum(p * incidReactionProdu 360 aHadron.SetMomentum(p * incidReactionProduct.GetMomentum() /
361 incidReactionProduct.GetTotalMomentum()) 361 incidReactionProduct.GetTotalMomentum());
362 } 362 }
363 else { 363 else {
364 iLevel = imaxEx; 364 iLevel = imaxEx;
365 } 365 }
366 366
367 if (theAngularDistribution[it] != nullptr) 367 if (theAngularDistribution[it] != nullptr) // MF4
368 { 368 {
369 if (theEnergyDistribution[it] != nullptr) 369 if (theEnergyDistribution[it] != nullptr) // MF5
370 { 370 {
371 //************************************** 371 //************************************************************
372 /* 372 /*
373 aHadron.SetKineticEnergy(theEnergy 373 aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
374 G4double eSecN = aHadron.GetKineti 374 G4double eSecN = aHadron.GetKineticEnergy();
375 */ 375 */
376 //************************************** 376 //************************************************************
377 // EMendoza --> maximum allowable energy 377 // EMendoza --> maximum allowable energy should be taken into account.
378 G4double dqi = 0.0; 378 G4double dqi = 0.0;
379 if (QI[it] < 0 || 849 < QI[it]) 379 if (QI[it] < 0 || 849 < QI[it])
380 dqi = QI[it]; // For backword compati 380 dqi = QI[it]; // For backword compatibility QI introduced since G4NDL3.15
381 G4double MaxEne = eKinetic + dqi; 381 G4double MaxEne = eKinetic + dqi;
382 G4double eSecN = 0.; 382 G4double eSecN = 0.;
383 383
384 G4int icounter = 0; 384 G4int icounter = 0;
385 G4int icounter_max = 1024; 385 G4int icounter_max = 1024;
386 G4int dummy = 0; 386 G4int dummy = 0;
387 do { 387 do {
388 ++icounter; 388 ++icounter;
389 if (icounter > icounter_max) { 389 if (icounter > icounter_max) {
390 G4cout << "Loop-counter exceeded the 390 G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
391 << __FILE__ << "." << G4endl; 391 << __FILE__ << "." << G4endl;
392 break; 392 break;
393 } 393 }
394 eSecN = theEnergyDistribution[it]->Sam 394 eSecN = theEnergyDistribution[it]->Sample(eKinetic, dummy);
395 } while (eSecN > MaxEne); // Loop check 395 } while (eSecN > MaxEne); // Loop checking, 11.05.2015, T. Koi
396 aHadron.SetKineticEnergy(eSecN); 396 aHadron.SetKineticEnergy(eSecN);
397 //************************************** 397 //************************************************************
398 eGamm = eKinetic - eSecN; 398 eGamm = eKinetic - eSecN;
399 for (iLevel = imaxEx; iLevel >= 0; --iLe 399 for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
400 if (theGammas.GetLevelEnergy(iLevel) < 400 if (theGammas.GetLevelEnergy(iLevel) < eGamm) break;
401 } 401 }
402 if (iLevel < imaxEx && iLevel >= 0) { 402 if (iLevel < imaxEx && iLevel >= 0) {
403 if (G4UniformRand() > 0.5) { 403 if (G4UniformRand() > 0.5) {
404 ++iLevel; 404 ++iLevel;
405 } 405 }
406 } 406 }
407 } 407 }
408 else { 408 else {
409 G4double eExcitation = 0; 409 G4double eExcitation = 0;
410 for (iLevel = imaxEx; iLevel >= 0; --iLe 410 for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
411 if (theGammas.GetLevelEnergy(iLevel) < 411 if (theGammas.GetLevelEnergy(iLevel) < eKinetic) break;
412 } 412 }
413 413
414 // Use QI value for calculating excitati 414 // Use QI value for calculating excitation energy of residual.
415 G4bool useQI = false; 415 G4bool useQI = false;
416 G4double dqi = QI[it]; 416 G4double dqi = QI[it];
417 if (dqi < 0 || 849 < dqi) useQI = true; 417 if (dqi < 0 || 849 < dqi) useQI = true; // Former libraries do not have values in this range
418 418
419 if (useQI) { 419 if (useQI) {
420 eExcitation = std::max(0., QI[0] - QI[ 420 eExcitation = std::max(0., QI[0] - QI[it]); // Bug fix 2333
421 421
422 // Re-evaluate iLevel based on this eE 422 // Re-evaluate iLevel based on this eExcitation
423 iLevel = 0; 423 iLevel = 0;
424 G4bool find = false; 424 G4bool find = false;
425 const G4double level_tolerance = 1.0 * 425 const G4double level_tolerance = 1.0 * CLHEP::keV;
426 426
427 // VI: the first level is ground 427 // VI: the first level is ground
428 if (0 < imaxEx) { 428 if (0 < imaxEx) {
429 for (iLevel = 1; iLevel <= imaxEx; + 429 for (iLevel = 1; iLevel <= imaxEx; ++iLevel) {
430 G4double elevel = theGammas.GetLev 430 G4double elevel = theGammas.GetLevelEnergy(iLevel);
431 if (std::abs(eExcitation - elevel) 431 if (std::abs(eExcitation - elevel) < level_tolerance) {
432 find = true; 432 find = true;
433 break; 433 break;
434 } 434 }
435 if (eExcitation < elevel) { 435 if (eExcitation < elevel) {
436 find = true; 436 find = true;
437 iLevel = std::max(iLevel - 1, 0) 437 iLevel = std::max(iLevel - 1, 0);
438 break; 438 break;
439 } 439 }
440 } 440 }
441 441
442 // If proper level cannot be found, 442 // If proper level cannot be found, use the maximum level
443 if (!find) iLevel = imaxEx; 443 if (!find) iLevel = imaxEx;
444 } 444 }
445 } 445 }
446 446
447 if (fManager->GetDEBUG() && eKinetic - e 447 if (fManager->GetDEBUG() && eKinetic - eExcitation < 0) {
448 throw G4HadronicException( 448 throw G4HadronicException(
449 __FILE__, __LINE__, 449 __FILE__, __LINE__,
450 "SEVERE: InelasticCompFS: Consistenc 450 "SEVERE: InelasticCompFS: Consistency of data not good enough, please file report");
451 } 451 }
452 if (eKinetic - eExcitation < 0) eExcitat 452 if (eKinetic - eExcitation < 0) eExcitation = 0;
453 if (iLevel != -1) aHadron.SetKineticEner 453 if (iLevel != -1) aHadron.SetKineticEnergy(eKinetic - eExcitation);
454 } 454 }
455 theAngularDistribution[it]->SampleAndUpdat 455 theAngularDistribution[it]->SampleAndUpdate(aHadron);
456 456
457 if (theFinalStatePhotons[it] == nullptr) { 457 if (theFinalStatePhotons[it] == nullptr) {
458 thePhotons = theGammas.GetDecayGammas(iL 458 thePhotons = theGammas.GetDecayGammas(iLevel);
459 eGamm -= theGammas.GetLevelEnergy(iLevel 459 eGamm -= theGammas.GetLevelEnergy(iLevel);
460 } 460 }
461 } 461 }
462 else if (theEnergyAngData[it] != nullptr) / 462 else if (theEnergyAngData[it] != nullptr) // MF6
463 { 463 {
464 theParticles = theEnergyAngData[it]->Sampl 464 theParticles = theEnergyAngData[it]->Sample(eKinetic);
465 465
466 // Adjust A and Z in the case of miss much 466 // Adjust A and Z in the case of miss much between selected data and target nucleus
467 if (theParticles != nullptr) { 467 if (theParticles != nullptr) {
468 G4int sumA = 0; 468 G4int sumA = 0;
469 G4int sumZ = 0; 469 G4int sumZ = 0;
470 G4int maxA = 0; 470 G4int maxA = 0;
471 G4int jAtMaxA = 0; 471 G4int jAtMaxA = 0;
472 for (G4int j = 0; j != (G4int)theParticl 472 for (G4int j = 0; j != (G4int)theParticles->size(); ++j) {
473 auto ptr = theParticles->at(j); 473 auto ptr = theParticles->at(j);
474 G4int barnum = ptr->GetDefinition()->GetBary 474 G4int barnum = ptr->GetDefinition()->GetBaryonNumber();
475 if (barnum > maxA) { 475 if (barnum > maxA) {
476 maxA = barnum; 476 maxA = barnum;
477 jAtMaxA = j; 477 jAtMaxA = j;
478 } 478 }
479 sumA += barnum; 479 sumA += barnum;
480 sumZ += G4lrint(ptr->GetDefinition()-> 480 sumZ += G4lrint(ptr->GetDefinition()->GetPDGCharge()/CLHEP::eplus);
481 } 481 }
482 G4int dA = theBaseA + hadProjectile->Get 482 G4int dA = theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
483 G4int dZ = theBaseZ + 483 G4int dZ = theBaseZ +
484 G4lrint(hadProjectile->GetDefinition()->GetP 484 G4lrint(hadProjectile->GetDefinition()->GetPDGCharge()/CLHEP::eplus) - sumZ;
485 if (dA < 0 || dZ < 0) { 485 if (dA < 0 || dZ < 0) {
486 G4int newA = theParticles->at(jAtMaxA) 486 G4int newA = theParticles->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA;
487 G4int newZ = 487 G4int newZ =
488 G4lrint(theParticles->at(jAtMaxA)->GetDefi 488 G4lrint(theParticles->at(jAtMaxA)->GetDefinition()->GetPDGCharge()/CLHEP::eplus) + dZ;
489 G4ParticleDefinition* pd = ionTable->G 489 G4ParticleDefinition* pd = ionTable->GetIon(newZ, newA);
490 theParticles->at(jAtMaxA)->SetDefiniti 490 theParticles->at(jAtMaxA)->SetDefinition(pd);
491 } 491 }
492 } 492 }
493 } 493 }
494 else { 494 else {
495 // @@@ what to do, if we have photon data, 495 // @@@ what to do, if we have photon data, but no info on the hadron itself
496 nothingWasKnownOnHadron = 1; 496 nothingWasKnownOnHadron = 1;
497 } 497 }
498 498
499 if (theFinalStatePhotons[it] != nullptr) { 499 if (theFinalStatePhotons[it] != nullptr) {
500 // the photon distributions are in the Nuc 500 // the photon distributions are in the Nucleus rest frame.
501 // TK residual rest frame 501 // TK residual rest frame
502 G4ReactionProduct boosted_tmp; 502 G4ReactionProduct boosted_tmp;
503 boosted_tmp.Lorentz(incidReactionProduct, 503 boosted_tmp.Lorentz(incidReactionProduct, theTarget);
504 G4double anEnergy = boosted_tmp.GetKinetic 504 G4double anEnergy = boosted_tmp.GetKineticEnergy();
505 thePhotons = theFinalStatePhotons[it]->Get 505 thePhotons = theFinalStatePhotons[it]->GetPhotons(anEnergy);
506 G4double aBaseEnergy = theFinalStatePhoton 506 G4double aBaseEnergy = theFinalStatePhotons[it]->GetLevelEnergy();
507 G4double testEnergy = 0; 507 G4double testEnergy = 0;
508 if (thePhotons != nullptr && !thePhotons-> 508 if (thePhotons != nullptr && !thePhotons->empty()) {
509 aBaseEnergy -= (*thePhotons)[0]->GetTota 509 aBaseEnergy -= (*thePhotons)[0]->GetTotalEnergy();
510 } 510 }
511 if (theFinalStatePhotons[it]->NeedsCascade 511 if (theFinalStatePhotons[it]->NeedsCascade()) {
512 while (aBaseEnergy > 0.01 * CLHEP::keV) 512 while (aBaseEnergy > 0.01 * CLHEP::keV) // Loop checking, 11.05.2015, T. Koi
513 { 513 {
514 // cascade down the levels 514 // cascade down the levels
515 G4bool foundMatchingLevel = false; 515 G4bool foundMatchingLevel = false;
516 G4int closest = 2; 516 G4int closest = 2;
517 G4double deltaEold = -1; 517 G4double deltaEold = -1;
518 for (G4int j = 1; j < it; ++j) { 518 for (G4int j = 1; j < it; ++j) {
519 if (theFinalStatePhotons[j] != nullp 519 if (theFinalStatePhotons[j] != nullptr) {
520 testEnergy = theFinalStatePhotons[ 520 testEnergy = theFinalStatePhotons[j]->GetLevelEnergy();
521 } 521 }
522 else { 522 else {
523 testEnergy = 0; 523 testEnergy = 0;
524 } 524 }
525 G4double deltaE = std::abs(testEnerg 525 G4double deltaE = std::abs(testEnergy - aBaseEnergy);
526 if (deltaE < 0.1 * CLHEP::keV) { 526 if (deltaE < 0.1 * CLHEP::keV) {
527 G4ReactionProductVector* theNext = 527 G4ReactionProductVector* theNext = theFinalStatePhotons[j]->GetPhotons(anEnergy);
528 if (thePhotons != nullptr) thePhot 528 if (thePhotons != nullptr) thePhotons->push_back(theNext->operator[](0));
529 aBaseEnergy = testEnergy - theNext 529 aBaseEnergy = testEnergy - theNext->operator[](0)->GetTotalEnergy();
530 delete theNext; 530 delete theNext;
531 foundMatchingLevel = true; 531 foundMatchingLevel = true;
532 break; // ===> 532 break; // ===>
533 } 533 }
534 if (theFinalStatePhotons[j] != nullp 534 if (theFinalStatePhotons[j] != nullptr && (deltaE < deltaEold || deltaEold < 0.)) {
535 closest = j; 535 closest = j;
536 deltaEold = deltaE; 536 deltaEold = deltaE;
537 } 537 }
538 } // <=== the break goes here. 538 } // <=== the break goes here.
539 if (!foundMatchingLevel) { 539 if (!foundMatchingLevel) {
540 G4ReactionProductVector* theNext = t 540 G4ReactionProductVector* theNext = theFinalStatePhotons[closest]->GetPhotons(anEnergy);
541 if (thePhotons != nullptr) thePhoton 541 if (thePhotons != nullptr) thePhotons->push_back(theNext->operator[](0));
542 aBaseEnergy = aBaseEnergy - theNext- 542 aBaseEnergy = aBaseEnergy - theNext->operator[](0)->GetTotalEnergy();
543 delete theNext; 543 delete theNext;
544 } 544 }
545 } 545 }
546 } 546 }
547 } 547 }
548 548
549 if (thePhotons != nullptr) { 549 if (thePhotons != nullptr) {
550 for (auto const & p : *thePhotons) { 550 for (auto const & p : *thePhotons) {
551 // back to lab 551 // back to lab
552 p->Lorentz(*p, -1. * theTarget); 552 p->Lorentz(*p, -1. * theTarget);
553 } 553 }
554 } 554 }
555 if (nothingWasKnownOnHadron != 0) { 555 if (nothingWasKnownOnHadron != 0) {
556 // In this case, hadron should be isotropi 556 // In this case, hadron should be isotropic in CM
557 // Next 12 lines are Emilio's replacement 557 // Next 12 lines are Emilio's replacement
558 // G4double QM=(incidReactionProduct.GetMa 558 // G4double QM=(incidReactionProduct.GetMass()+targetMass)-(aHadron.GetMass()+residualMass);
559 // G4double eExcitation = QM-QI[it]; 559 // G4double eExcitation = QM-QI[it];
560 // G4double eExcitation = QI[0] - QI[it]; 560 // G4double eExcitation = QI[0] - QI[it]; // Fix of bug #1838
561 // if(eExcitation<20*CLHEP::keV){eExcitati 561 // if(eExcitation<20*CLHEP::keV){eExcitation=0;}
562 562
563 G4double eExcitation = std::max(0., QI[0] 563 G4double eExcitation = std::max(0., QI[0] - QI[it]); // Fix of bug #2333
564 564
565 two_body_reaction(&incidReactionProduct, & 565 two_body_reaction(&incidReactionProduct, &theTarget, &aHadron, eExcitation);
566 if (thePhotons == nullptr && eExcitation > 566 if (thePhotons == nullptr && eExcitation > 0) {
567 for (iLevel = imaxEx; iLevel >= 0; --iLe 567 for (iLevel = imaxEx; iLevel >= 0; --iLevel) {
568 if (theGammas.GetLevelEnergy(iLevel) < 568 if (theGammas.GetLevelEnergy(iLevel) < eExcitation + 5 * keV) break; // 5 keV tolerance
569 } 569 }
570 thePhotons = theGammas.GetDecayGammas(iL 570 thePhotons = theGammas.GetDecayGammas(iLevel);
571 } 571 }
572 } 572 }
573 573
574 // fill the result 574 // fill the result
575 // Beware - the recoil is not necessarily in 575 // Beware - the recoil is not necessarily in the particles...
576 // Can be calculated from momentum conservat 576 // Can be calculated from momentum conservation?
577 // The idea is that the particles ar emitted 577 // The idea is that the particles ar emitted forst, and the gammas only once the
578 // recoil is on the residual; assumption is 578 // recoil is on the residual; assumption is that gammas do not contribute to
579 // the recoil. 579 // the recoil.
580 // This needs more design @@@ 580 // This needs more design @@@
581 581
582 G4bool needsSeparateRecoil = false; 582 G4bool needsSeparateRecoil = false;
583 G4int totalBaryonNumber = 0; 583 G4int totalBaryonNumber = 0;
584 G4int totalCharge = 0; 584 G4int totalCharge = 0;
585 G4ThreeVector totalMomentum(0); 585 G4ThreeVector totalMomentum(0);
586 if (theParticles != nullptr) { 586 if (theParticles != nullptr) {
587 const G4ParticleDefinition* aDef; 587 const G4ParticleDefinition* aDef;
588 for (std::size_t ii0 = 0; ii0 < theParticl 588 for (std::size_t ii0 = 0; ii0 < theParticles->size(); ++ii0) {
589 aDef = (*theParticles)[ii0]->GetDefiniti 589 aDef = (*theParticles)[ii0]->GetDefinition();
590 totalBaryonNumber += aDef->GetBaryonNumb 590 totalBaryonNumber += aDef->GetBaryonNumber();
591 totalCharge += G4lrint(aDef->GetPDGCharg 591 totalCharge += G4lrint(aDef->GetPDGCharge()/CLHEP::eplus);
592 totalMomentum += (*theParticles)[ii0]->G 592 totalMomentum += (*theParticles)[ii0]->GetMomentum();
593 } 593 }
594 if (totalBaryonNumber 594 if (totalBaryonNumber
595 != theBaseA + hadProjectile->GetDefin 595 != theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber())
596 { 596 {
597 needsSeparateRecoil = true; 597 needsSeparateRecoil = true;
598 residualA = theBaseA + hadProjectile->Ge 598 residualA = theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber()
599 - totalBaryonNumber; 599 - totalBaryonNumber;
600 residualZ = theBaseZ + 600 residualZ = theBaseZ +
601 G4lrint((hadProjectile->GetDefinition()->Get 601 G4lrint((hadProjectile->GetDefinition()->GetPDGCharge() - totalCharge)/CLHEP::eplus);
602 } 602 }
603 } 603 }
604 604
605 std::size_t nPhotons = 0; 605 std::size_t nPhotons = 0;
606 if (thePhotons != nullptr) { 606 if (thePhotons != nullptr) {
607 nPhotons = thePhotons->size(); 607 nPhotons = thePhotons->size();
608 } 608 }
609 609
610 G4DynamicParticle* theSec; 610 G4DynamicParticle* theSec;
611 611
612 if (theParticles == nullptr) { 612 if (theParticles == nullptr) {
613 theSec = new G4DynamicParticle; 613 theSec = new G4DynamicParticle;
614 theSec->SetDefinition(aHadron.GetDefinitio 614 theSec->SetDefinition(aHadron.GetDefinition());
615 theSec->SetMomentum(aHadron.GetMomentum()) 615 theSec->SetMomentum(aHadron.GetMomentum());
616 theResult.Get()->AddSecondary(theSec, secI 616 theResult.Get()->AddSecondary(theSec, secID);
617 #ifdef G4VERBOSE 617 #ifdef G4VERBOSE
618 if (fManager->GetDEBUG()) 618 if (fManager->GetDEBUG())
619 G4cout << " G4ParticleHPInelasticCompFS: 619 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary1 "
620 << theSec->GetParticleDefinition( 620 << theSec->GetParticleDefinition()->GetParticleName()
621 << " E= " << theSec->GetKineticEn 621 << " E= " << theSec->GetKineticEnergy() << " NSECO "
622 << theResult.Get()->GetNumberOfSe 622 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
623 #endif 623 #endif
624 624
625 aHadron.Lorentz(aHadron, theTarget); 625 aHadron.Lorentz(aHadron, theTarget);
626 G4ReactionProduct theResidual; 626 G4ReactionProduct theResidual;
627 theResidual.SetDefinition(ionTable->GetIon 627 theResidual.SetDefinition(ionTable->GetIon(residualZ, residualA, 0));
628 theResidual.SetKineticEnergy(aHadron.GetKi 628 theResidual.SetKineticEnergy(aHadron.GetKineticEnergy() * aHadron.GetMass()
629 / theResidual 629 / theResidual.GetMass());
630 630
631 // 080612TK contribution from Benoit Pirar 631 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
632 // theResidual.SetMomentum(-1.*aHadron.Get 632 // theResidual.SetMomentum(-1.*aHadron.GetMomentum());
633 G4ThreeVector incidentNeutronMomentum = in 633 G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
634 theResidual.SetMomentum(incidentNeutronMom 634 theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
635 635
636 theResidual.Lorentz(theResidual, -1. * the 636 theResidual.Lorentz(theResidual, -1. * theTarget);
637 G4ThreeVector totalPhotonMomentum(0, 0, 0) 637 G4ThreeVector totalPhotonMomentum(0, 0, 0);
638 if (thePhotons != nullptr) { 638 if (thePhotons != nullptr) {
639 for (std::size_t i = 0; i < nPhotons; ++ 639 for (std::size_t i = 0; i < nPhotons; ++i) {
640 totalPhotonMomentum += (*thePhotons)[i 640 totalPhotonMomentum += (*thePhotons)[i]->GetMomentum();
641 } 641 }
642 } 642 }
643 theSec = new G4DynamicParticle; 643 theSec = new G4DynamicParticle;
644 theSec->SetDefinition(theResidual.GetDefin 644 theSec->SetDefinition(theResidual.GetDefinition());
645 theSec->SetMomentum(theResidual.GetMomentu 645 theSec->SetMomentum(theResidual.GetMomentum() - totalPhotonMomentum);
646 theResult.Get()->AddSecondary(theSec, secI 646 theResult.Get()->AddSecondary(theSec, secID);
647 #ifdef G4VERBOSE 647 #ifdef G4VERBOSE
648 if (fManager->GetDEBUG()) 648 if (fManager->GetDEBUG())
649 G4cout << this << " G4ParticleHPInelasti 649 G4cout << this << " G4ParticleHPInelasticCompFS::BaseApply add secondary2 "
650 << theSec->GetParticleDefinition( 650 << theSec->GetParticleDefinition()->GetParticleName()
651 << " E= " << theSec->GetKineticEn 651 << " E= " << theSec->GetKineticEnergy() << " NSECO "
652 << theResult.Get()->GetNumberOfSe 652 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
653 #endif 653 #endif
654 } 654 }
655 else { 655 else {
656 for (std::size_t i0 = 0; i0 < theParticles 656 for (std::size_t i0 = 0; i0 < theParticles->size(); ++i0) {
657 theSec = new G4DynamicParticle; 657 theSec = new G4DynamicParticle;
658 theSec->SetDefinition((*theParticles)[i0 658 theSec->SetDefinition((*theParticles)[i0]->GetDefinition());
659 theSec->SetMomentum((*theParticles)[i0]- 659 theSec->SetMomentum((*theParticles)[i0]->GetMomentum());
660 theResult.Get()->AddSecondary(theSec, se 660 theResult.Get()->AddSecondary(theSec, secID);
661 #ifdef G4VERBOSE 661 #ifdef G4VERBOSE
662 if (fManager->GetDEBUG()) 662 if (fManager->GetDEBUG())
663 G4cout << " G4ParticleHPInelasticCompF 663 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary3 "
664 << theSec->GetParticleDefinitio 664 << theSec->GetParticleDefinition()->GetParticleName()
665 << " E= " << theSec->GetKinetic 665 << " E= " << theSec->GetKineticEnergy() << " NSECO "
666 << theResult.Get()->GetNumberOf 666 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
667 #endif 667 #endif
668 delete (*theParticles)[i0]; 668 delete (*theParticles)[i0];
669 } 669 }
670 delete theParticles; 670 delete theParticles;
671 if (needsSeparateRecoil && residualZ != 0) 671 if (needsSeparateRecoil && residualZ != 0) {
672 G4ReactionProduct theResidual; 672 G4ReactionProduct theResidual;
673 theResidual.SetDefinition(ionTable->GetI 673 theResidual.SetDefinition(ionTable->GetIon(residualZ, residualA, 0));
674 G4double resiualKineticEnergy = theResid 674 G4double resiualKineticEnergy = theResidual.GetMass() * theResidual.GetMass();
675 resiualKineticEnergy += totalMomentum * 675 resiualKineticEnergy += totalMomentum * totalMomentum;
676 resiualKineticEnergy = std::sqrt(resiual 676 resiualKineticEnergy = std::sqrt(resiualKineticEnergy) - theResidual.GetMass();
677 theResidual.SetKineticEnergy(resiualKine 677 theResidual.SetKineticEnergy(resiualKineticEnergy);
678 678
679 // 080612TK contribution from Benoit Pir 679 // 080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
680 // theResidual.SetMomentum(-1.*totalMome 680 // theResidual.SetMomentum(-1.*totalMomentum);
681 // G4ThreeVector incidentNeutronMomentum 681 // G4ThreeVector incidentNeutronMomentum = incidReactionProduct.GetMomentum();
682 // theResidual.SetMomentum(incidentNeutr 682 // theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
683 // 080717 TK Comment still do NOT includ 683 // 080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
684 theResidual.SetMomentum(incidReactionPro 684 theResidual.SetMomentum(incidReactionProduct.GetMomentum() + theTarget.GetMomentum()
685 - totalMomentum) 685 - totalMomentum);
686 686
687 theSec = new G4DynamicParticle; 687 theSec = new G4DynamicParticle;
688 theSec->SetDefinition(theResidual.GetDef 688 theSec->SetDefinition(theResidual.GetDefinition());
689 theSec->SetMomentum(theResidual.GetMomen 689 theSec->SetMomentum(theResidual.GetMomentum());
690 theResult.Get()->AddSecondary(theSec, se 690 theResult.Get()->AddSecondary(theSec, secID);
691 #ifdef G4VERBOSE 691 #ifdef G4VERBOSE
692 if (fManager->GetDEBUG()) 692 if (fManager->GetDEBUG())
693 G4cout << " G4ParticleHPInelasticCompF 693 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary4 "
694 << theSec->GetParticleDefinitio 694 << theSec->GetParticleDefinition()->GetParticleName()
695 << " E= " << theSec->GetKinetic 695 << " E= " << theSec->GetKineticEnergy() << " NSECO "
696 << theResult.Get()->GetNumberOf 696 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
697 #endif 697 #endif
698 } 698 }
699 } 699 }
700 if (thePhotons != nullptr) { 700 if (thePhotons != nullptr) {
701 for (std::size_t i = 0; i < nPhotons; ++i) 701 for (std::size_t i = 0; i < nPhotons; ++i) {
702 theSec = new G4DynamicParticle; 702 theSec = new G4DynamicParticle;
703 // Bug reported Chao Zhang (Chao.Zhang@u 703 // Bug reported Chao Zhang (Chao.Zhang@usd.edu), Dongming Mei(Dongming.Mei@usd.edu) Feb. 25,
704 // 2009 theSec->SetDefinition(G4Gamma::G 704 // 2009 theSec->SetDefinition(G4Gamma::Gamma());
705 theSec->SetDefinition((*thePhotons)[i]-> 705 theSec->SetDefinition((*thePhotons)[i]->GetDefinition());
706 // But never cause real effect at least 706 // But never cause real effect at least with G4NDL3.13 TK
707 theSec->SetMomentum((*thePhotons)[i]->Ge 707 theSec->SetMomentum((*thePhotons)[i]->GetMomentum());
708 theResult.Get()->AddSecondary(theSec, se 708 theResult.Get()->AddSecondary(theSec, secID);
709 #ifdef G4VERBOSE 709 #ifdef G4VERBOSE
710 if (fManager->GetDEBUG()) 710 if (fManager->GetDEBUG())
711 G4cout << " G4ParticleHPInelasticCompF 711 G4cout << " G4ParticleHPInelasticCompFS::BaseApply add secondary5 "
712 << theSec->GetParticleDefinitio 712 << theSec->GetParticleDefinition()->GetParticleName()
713 << " E= " << theSec->GetKinetic 713 << " E= " << theSec->GetKineticEnergy() << " NSECO "
714 << theResult.Get()->GetNumberOf 714 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
715 #endif 715 #endif
716 716
717 delete thePhotons->operator[](i); 717 delete thePhotons->operator[](i);
718 } 718 }
719 // some garbage collection 719 // some garbage collection
720 delete thePhotons; 720 delete thePhotons;
721 } 721 }
722 722
723 G4ParticleDefinition* targ_pd = ionTable->Ge 723 G4ParticleDefinition* targ_pd = ionTable->GetIon(theBaseZ, theBaseA, 0.0);
724 G4LorentzVector targ_4p_lab( 724 G4LorentzVector targ_4p_lab(
725 theTarget.GetMomentum(), 725 theTarget.GetMomentum(),
726 std::sqrt(targ_pd->GetPDGMass() * targ_pd- 726 std::sqrt(targ_pd->GetPDGMass() * targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2()));
727 G4LorentzVector proj_4p_lab = theTrack.Get4M 727 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
728 G4LorentzVector init_4p_lab = proj_4p_lab + 728 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
729 adjust_final_state(init_4p_lab); 729 adjust_final_state(init_4p_lab);
730 730
731 // clean up the primary neutron 731 // clean up the primary neutron
732 theResult.Get()->SetStatusChange(stopAndKill 732 theResult.Get()->SetStatusChange(stopAndKill);
733 } 733 }
734 734
735 // Re-implemented by E. Mendoza (2019). Isotro 735 // Re-implemented by E. Mendoza (2019). Isotropic emission in the CMS:
736 // proj: projectile in target-rest-frame (inp 736 // proj: projectile in target-rest-frame (input)
737 // targ: target in target-rest-frame (input) 737 // targ: target in target-rest-frame (input)
738 // product: secondary particle in target-rest 738 // product: secondary particle in target-rest-frame (output)
739 // resExcitationEnergy: excitation energy of 739 // resExcitationEnergy: excitation energy of the residual nucleus
740 // 740 //
741 void G4ParticleHPInelasticCompFS::two_body_rea 741 void G4ParticleHPInelasticCompFS::two_body_reaction(G4ReactionProduct* proj,
742 742 G4ReactionProduct* targ,
743 743 G4ReactionProduct* product,
744 744 G4double resExcitationEnergy)
745 { 745 {
746 // CMS system: 746 // CMS system:
747 G4ReactionProduct theCMS = *proj + *targ; 747 G4ReactionProduct theCMS = *proj + *targ;
748 748
749 // Residual definition: 749 // Residual definition:
750 G4int resZ = G4lrint((proj->GetDefinition()- 750 G4int resZ = G4lrint((proj->GetDefinition()->GetPDGCharge() + targ->GetDefinition()->GetPDGCharge()
751 - product->GetDefinition()->GetPDGCharge 751 - product->GetDefinition()->GetPDGCharge())/CLHEP::eplus);
752 G4int resA = proj->GetDefinition()->GetBaryo 752 G4int resA = proj->GetDefinition()->GetBaryonNumber() + targ->GetDefinition()->GetBaryonNumber()
753 - product->GetDefinition()->Get 753 - product->GetDefinition()->GetBaryonNumber();
754 G4ReactionProduct theResidual; 754 G4ReactionProduct theResidual;
755 theResidual.SetDefinition(ionTable->GetIon(r 755 theResidual.SetDefinition(ionTable->GetIon(resZ, resA, 0.0));
756 756
757 // CMS system: 757 // CMS system:
758 G4ReactionProduct theCMSproj; 758 G4ReactionProduct theCMSproj;
759 G4ReactionProduct theCMStarg; 759 G4ReactionProduct theCMStarg;
760 theCMSproj.Lorentz(*proj, theCMS); 760 theCMSproj.Lorentz(*proj, theCMS);
761 theCMStarg.Lorentz(*targ, theCMS); 761 theCMStarg.Lorentz(*targ, theCMS);
762 // final Momentum in the CMS: 762 // final Momentum in the CMS:
763 G4double totE = std::sqrt(theCMSproj.GetMass 763 G4double totE = std::sqrt(theCMSproj.GetMass() * theCMSproj.GetMass()
764 + theCMSproj.GetTo 764 + theCMSproj.GetTotalMomentum() * theCMSproj.GetTotalMomentum())
765 + std::sqrt(theCMStarg.GetMa 765 + std::sqrt(theCMStarg.GetMass() * theCMStarg.GetMass()
766 + theCMStarg.Get 766 + theCMStarg.GetTotalMomentum() * theCMStarg.GetTotalMomentum());
767 G4double prodmass = product->GetMass(); 767 G4double prodmass = product->GetMass();
768 G4double resmass = theResidual.GetMass() + r 768 G4double resmass = theResidual.GetMass() + resExcitationEnergy;
769 G4double fmomsquared = (totE * totE - (prodm 769 G4double fmomsquared = (totE * totE - (prodmass - resmass) * (prodmass - resmass)) *
770 (totE * totE - (prodmass + resmass) * (pro 770 (totE * totE - (prodmass + resmass) * (prodmass + resmass)) / (4.*totE*totE);
771 G4double fmom = (fmomsquared > 0) ? std::sqr 771 G4double fmom = (fmomsquared > 0) ? std::sqrt(fmomsquared) : 0.0;
772 772
773 // random (isotropic direction): 773 // random (isotropic direction):
774 product->SetMomentum(fmom * G4RandomDirectio 774 product->SetMomentum(fmom * G4RandomDirection());
775 product->SetTotalEnergy(std::sqrt(prodmass * 775 product->SetTotalEnergy(std::sqrt(prodmass * prodmass + fmom * fmom)); // CMS
776 // Back to the LAB system: 776 // Back to the LAB system:
777 product->Lorentz(*product, -1. * theCMS); 777 product->Lorentz(*product, -1. * theCMS);
778 } 778 }
779 779
780 G4bool G4ParticleHPInelasticCompFS::use_nresp7 780 G4bool G4ParticleHPInelasticCompFS::use_nresp71_model(const G4ParticleDefinition* aDefinition,
781 781 const G4int itt,
782 782 const G4ReactionProduct& theTarget,
783 783 G4ReactionProduct& boosted)
784 { 784 {
785 if (aDefinition == G4Neutron::Definition()) 785 if (aDefinition == G4Neutron::Definition()) // If the outgoing particle is a neutron...
786 { 786 {
787 // LR: flag LR in ENDF. It indicates wheth 787 // LR: flag LR in ENDF. It indicates whether there is breakup of the residual nucleus or not.
788 // it: exit channel (index of the carbon e 788 // it: exit channel (index of the carbon excited state)
789 789
790 // Added by A. R. Garcia (CIEMAT) to inclu 790 // Added by A. R. Garcia (CIEMAT) to include the physics of C(N,N'3A) reactions from NRESP71.
791 791
792 if (LR[itt] > 0) // If there is breakup of 792 if (LR[itt] > 0) // If there is breakup of the residual nucleus LR(flag LR in ENDF)>0 (i.e. Z=6
793 // MT=52-91 (it=MT-50)). 793 // MT=52-91 (it=MT-50)).
794 { 794 {
795 // Defining carbon as the target in the 795 // Defining carbon as the target in the reference frame at rest.
796 G4ReactionProduct theCarbon(theTarget); 796 G4ReactionProduct theCarbon(theTarget);
797 797
798 theCarbon.SetMomentum(G4ThreeVector()); 798 theCarbon.SetMomentum(G4ThreeVector());
799 theCarbon.SetKineticEnergy(0.); 799 theCarbon.SetKineticEnergy(0.);
800 800
801 // Creating four reaction products. 801 // Creating four reaction products.
802 G4ReactionProduct theProds[4]; 802 G4ReactionProduct theProds[4];
803 803
804 // Applying C(N,N'3A) reaction mechanism 804 // Applying C(N,N'3A) reaction mechanisms in the target rest frame.
805 if (itt == 41) { 805 if (itt == 41) {
806 // QI=QM=-7.275 MeV for C-0(N,N')C-C(3 806 // QI=QM=-7.275 MeV for C-0(N,N')C-C(3A) in ENDF/B-VII.1.
807 // This is not the value of the QI of 807 // This is not the value of the QI of the first step according
808 // to the model. So we don't take it. 808 // to the model. So we don't take it. Instead, we set the one
809 // we have calculated: QI=(mn+m12C)-(m 809 // we have calculated: QI=(mn+m12C)-(ma+m9Be+Ex9Be)=-8.130 MeV.
810 nresp71_model.ApplyMechanismI_NBeA2A(b 810 nresp71_model.ApplyMechanismI_NBeA2A(boosted, theCarbon, theProds, -8.130 /*QI[it]*/);
811 // N+C --> A[0]+9BE* | 9BE* --> N[1]+8 811 // N+C --> A[0]+9BE* | 9BE* --> N[1]+8BE | 8BE --> 2*A[2,3].
812 } 812 }
813 else { 813 else {
814 nresp71_model.ApplyMechanismII_ACN2A(b 814 nresp71_model.ApplyMechanismII_ACN2A(boosted, theCarbon, theProds, QI[itt]);
815 // N+C --> N'[0]+C* | C* --> A[1]+8BE 815 // N+C --> N'[0]+C* | C* --> A[1]+8BE | 8BE --> 2*A[2,3].
816 } 816 }
817 817
818 // Returning to the reference frame wher 818 // Returning to the reference frame where the target was in motion.
819 for (auto& theProd : theProds) { 819 for (auto& theProd : theProds) {
820 theProd.Lorentz(theProd, -1. * theTarg 820 theProd.Lorentz(theProd, -1. * theTarget);
821 theResult.Get()->AddSecondary( 821 theResult.Get()->AddSecondary(
822 new G4DynamicParticle(theProd.GetDef 822 new G4DynamicParticle(theProd.GetDefinition(), theProd.GetMomentum()), secID);
823 } 823 }
824 824
825 // Killing the primary neutron. 825 // Killing the primary neutron.
826 theResult.Get()->SetStatusChange(stopAnd 826 theResult.Get()->SetStatusChange(stopAndKill);
827 827
828 return true; 828 return true;
829 } 829 }
830 } 830 }
831 else if (aDefinition == G4Alpha::Definition( 831 else if (aDefinition == G4Alpha::Definition()) // If the outgoing particle is an alpha, ...
832 { 832 {
833 // Added by A. R. Garcia (CIEMAT) to inclu 833 // Added by A. R. Garcia (CIEMAT) to include the physics of C(N,A)9BE reactions from NRESP71.
834 834
835 if (LR[itt] == 0) // If Z=6, an alpha part 835 if (LR[itt] == 0) // If Z=6, an alpha particle is emitted and there is no breakup of the
836 // residual nucleus LR(f 836 // residual nucleus LR(flag LR in ENDF)==0.
837 { 837 {
838 // Defining carbon as the target in the 838 // Defining carbon as the target in the reference frame at rest.
839 G4ReactionProduct theCarbon(theTarget); 839 G4ReactionProduct theCarbon(theTarget);
840 theCarbon.SetMomentum(G4ThreeVector()); 840 theCarbon.SetMomentum(G4ThreeVector());
841 theCarbon.SetKineticEnergy(0.); 841 theCarbon.SetKineticEnergy(0.);
842 842
843 // Creating four reaction products. 843 // Creating four reaction products.
844 G4ReactionProduct theProds[2]; 844 G4ReactionProduct theProds[2];
845 845
846 // Applying C(N,A)9BE reaction mechanism 846 // Applying C(N,A)9BE reaction mechanism.
847 nresp71_model.ApplyMechanismABE(boosted, 847 nresp71_model.ApplyMechanismABE(boosted, theCarbon, theProds);
848 // N+C --> A[0]+9BE[1]. 848 // N+C --> A[0]+9BE[1].
849 849
850 for (auto& theProd : theProds) { 850 for (auto& theProd : theProds) {
851 // Returning to the system of referenc 851 // Returning to the system of reference where the target was in motion.
852 theProd.Lorentz(theProd, -1. * theTarg 852 theProd.Lorentz(theProd, -1. * theTarget);
853 theResult.Get()->AddSecondary( 853 theResult.Get()->AddSecondary(
854 new G4DynamicParticle(theProd.GetDef 854 new G4DynamicParticle(theProd.GetDefinition(), theProd.GetMomentum()), secID);
855 } 855 }
856 856
857 // Killing the primary neutron. 857 // Killing the primary neutron.
858 theResult.Get()->SetStatusChange(stopAnd 858 theResult.Get()->SetStatusChange(stopAndKill);
859 return true; 859 return true;
860 } 860 }
861 G4Exception("G4ParticleHPInelasticCompFS:: 861 G4Exception("G4ParticleHPInelasticCompFS::CompositeApply()", "G4ParticleInelasticCompFS.cc",
862 FatalException, "Alpha product 862 FatalException, "Alpha production with LR!=0.");
863 } 863 }
864 return false; 864 return false;
865 } 865 }
866 866