| 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 // 080801 Give a warning message for irregular 30 // 080801 Give a warning message for irregular mass value in data file by T. Koi
31 // Introduce theNDLDataA,Z which has A 31 // Introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
32 // 081024 G4NucleiPropertiesTable:: to G4Nucle 32 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
33 // 101111 Add Special treatment for Be9(n,2n)B 33 // 101111 Add Special treatment for Be9(n,2n)Be8(2a) case by T. Koi
34 // 34 //
35 // P. Arce, June-2014 Conversion neutron_hp to 35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 // 36 //
37 // June-2019 - E. Mendoza --> Added protection 37 // June-2019 - E. Mendoza --> Added protection against residual with Z<0 or A<Z + adjust_final_state
38 // is not applied when data is in MF=6 format 38 // is not applied when data is in MF=6 format (no correlated particle emission) + bug correction
39 // (add Q value info to G4ParticleHPNBodyPhase 39 // (add Q value info to G4ParticleHPNBodyPhaseSpace).
40 40
41 #include "G4ParticleHPInelasticBaseFS.hh" 41 #include "G4ParticleHPInelasticBaseFS.hh"
42 42
43 #include "G4Alpha.hh" 43 #include "G4Alpha.hh"
44 #include "G4Electron.hh" 44 #include "G4Electron.hh"
45 #include "G4He3.hh" 45 #include "G4He3.hh"
46 #include "G4IonTable.hh" 46 #include "G4IonTable.hh"
47 #include "G4NucleiProperties.hh" 47 #include "G4NucleiProperties.hh"
48 #include "G4Nucleus.hh" 48 #include "G4Nucleus.hh"
49 #include "G4ParticleHPDataUsed.hh" 49 #include "G4ParticleHPDataUsed.hh"
50 #include "G4ParticleHPManager.hh" 50 #include "G4ParticleHPManager.hh"
51 51
52 G4ParticleHPInelasticBaseFS::G4ParticleHPInela <<
53 { <<
54 hasXsec = true; <<
55 theXsection = new G4ParticleHPVector; <<
56 } <<
57 <<
58 G4ParticleHPInelasticBaseFS::~G4ParticleHPInel <<
59 { <<
60 delete theXsection; <<
61 delete theEnergyDistribution; <<
62 delete theFinalStatePhotons; <<
63 delete theEnergyAngData; <<
64 delete theAngularDistribution; <<
65 } <<
66 <<
67 void G4ParticleHPInelasticBaseFS::InitGammas(G 52 void G4ParticleHPInelasticBaseFS::InitGammas(G4double AR, G4double ZR)
68 { 53 {
69 G4int Z = G4lrint(ZR); <<
70 G4int A = G4lrint(AR); <<
71 std::ostringstream ost; 54 std::ostringstream ost;
72 ost << gammaPath << "z" << Z << ".a" << A; << 55 ost << gammaPath << "z" << ZR << ".a" << AR;
73 G4String aName = ost.str(); 56 G4String aName = ost.str();
74 std::ifstream from(aName, std::ios::in); 57 std::ifstream from(aName, std::ios::in);
75 58
76 if (!from) return; // no data found for thi 59 if (!from) return; // no data found for this isotope
77 std::ifstream theGammaData(aName, std::ios:: 60 std::ifstream theGammaData(aName, std::ios::in);
78 61
79 theNuclearMassDifference = G4NucleiPropertie << 62 G4double eps = 0.001;
80 - G4NucleiProperties::GetBindingEnergy(the << 63 theNuclearMassDifference =
>> 64 G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(AR + eps), static_cast<G4int>(ZR + eps))
>> 65 - G4NucleiProperties::GetBindingEnergy(static_cast<G4int>(theBaseA + eps),
>> 66 static_cast<G4int>(theBaseZ + eps));
81 theGammas.Init(theGammaData); 67 theGammas.Init(theGammaData);
82 } 68 }
83 69
84 void G4ParticleHPInelasticBaseFS::Init(G4doubl << 70 void G4ParticleHPInelasticBaseFS::Init(G4double A, G4double Z, G4int M, G4String& dirName,
85 const G << 71 G4String& bit, G4ParticleDefinition*)
86 const G <<
87 { 72 {
88 gammaPath = fManager->GetNeutronHPPath() + " << 73 gammaPath = "/Inelastic/Gammas/";
>> 74 if (!G4FindDataDir("G4NEUTRONHPDATA"))
>> 75 throw G4HadronicException(__FILE__, __LINE__,
>> 76 "Please setenv G4NEUTRONHPDATA to point to the neutron cross-section "
>> 77 "files where Inelastic/Gammas data is found.");
>> 78 G4String tBase = G4FindDataDir("G4NEUTRONHPDATA");
>> 79 gammaPath = tBase + gammaPath;
89 G4String tString = dirName; 80 G4String tString = dirName;
90 SetA_Z(A, Z, M); << 81 G4bool dbool;
91 G4bool dbool = true; <<
92 G4ParticleHPDataUsed aFile = 82 G4ParticleHPDataUsed aFile =
93 theNames.GetName(theBaseA, theBaseZ, M, tS << 83 theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, tString, bit, dbool);
94 SetAZMs(aFile); <<
95 G4String filename = aFile.GetName(); 84 G4String filename = aFile.GetName();
96 #ifdef G4VERBOSE << 85 #ifdef G4PHPDEBUG
97 if (fManager->GetDEBUG()) << 86 if (std::getenv("G4ParticleHPDebug"))
98 G4cout << " G4ParticleHPInelasticBaseFS::I 87 G4cout << " G4ParticleHPInelasticBaseFS::Init FILE " << filename << G4endl;
99 #endif 88 #endif
>> 89 SetAZMs(A, Z, M, aFile);
100 90
101 if (!dbool || (theBaseZ <= 2 && (theNDLDataZ << 91 if (!dbool
>> 92 || (Z < 2.5 && (std::abs(theNDLDataZ - Z) > 0.0001 || std::abs(theNDLDataA - A) > 0.0001)))
102 { 93 {
103 #ifdef G4PHPDEBUG 94 #ifdef G4PHPDEBUG
104 if (fManager->GetDEBUG()) << 95 if (std::getenv("G4ParticleHPDebug_NamesLogging"))
105 G4cout << "Skipped = " << filename << " 96 G4cout << "Skipped = " << filename << " " << A << " " << Z << G4endl;
106 #endif 97 #endif
107 hasAnyData = false; 98 hasAnyData = false;
108 hasFSData = false; 99 hasFSData = false;
109 hasXsec = false; 100 hasXsec = false;
110 return; 101 return;
111 } 102 }
112 103
113 std::istringstream theData(std::ios::in); 104 std::istringstream theData(std::ios::in);
114 fManager->GetDataStream(filename, theData); << 105 G4ParticleHPManager::GetInstance()->GetDataStream(filename, theData);
115 106
116 if (!theData) //"!" is a operator of ios 107 if (!theData) //"!" is a operator of ios
117 { 108 {
118 hasAnyData = false; 109 hasAnyData = false;
119 hasFSData = false; 110 hasFSData = false;
120 hasXsec = false; 111 hasXsec = false;
>> 112 // theData.close();
121 return; // no data for exactly this isoto 113 return; // no data for exactly this isotope and FS
122 } 114 }
123 // here we go 115 // here we go
124 G4int infoType, dataType, dummy = INT_MAX; 116 G4int infoType, dataType, dummy = INT_MAX;
125 hasFSData = false; 117 hasFSData = false;
126 while (theData >> infoType) // Loop checkin 118 while (theData >> infoType) // Loop checking, 11.05.2015, T. Koi
127 { 119 {
128 theData >> dataType; 120 theData >> dataType;
129 121
130 if (dummy == INT_MAX) theData >> Qvalue >> 122 if (dummy == INT_MAX) theData >> Qvalue >> dummy;
131 Qvalue *= CLHEP::eV; 123 Qvalue *= CLHEP::eV;
132 // In G4NDL4.5 this value is the MT number 124 // In G4NDL4.5 this value is the MT number (<1000),
133 // in others is que Q-value in eV 125 // in others is que Q-value in eV
134 126
135 if (dataType == 3) { 127 if (dataType == 3) {
136 G4int total; 128 G4int total;
137 theData >> total; 129 theData >> total;
138 theXsection->Init(theData, total, CLHEP: 130 theXsection->Init(theData, total, CLHEP::eV);
139 } 131 }
140 else if (dataType == 4) { 132 else if (dataType == 4) {
141 theAngularDistribution = new G4ParticleH 133 theAngularDistribution = new G4ParticleHPAngular;
142 theAngularDistribution->Init(theData); 134 theAngularDistribution->Init(theData);
143 hasFSData = true; 135 hasFSData = true;
144 } 136 }
145 else if (dataType == 5) { 137 else if (dataType == 5) {
146 theEnergyDistribution = new G4ParticleHP 138 theEnergyDistribution = new G4ParticleHPEnergyDistribution;
147 theEnergyDistribution->Init(theData); 139 theEnergyDistribution->Init(theData);
148 hasFSData = true; 140 hasFSData = true;
149 } 141 }
150 else if (dataType == 6) { 142 else if (dataType == 6) {
151 theEnergyAngData = new G4ParticleHPEnAng 143 theEnergyAngData = new G4ParticleHPEnAngCorrelation(theProjectile);
152 theEnergyAngData->Init(theData); 144 theEnergyAngData->Init(theData);
153 hasFSData = true; 145 hasFSData = true;
154 } 146 }
155 else if (dataType == 12) { 147 else if (dataType == 12) {
156 theFinalStatePhotons = new G4ParticleHPP 148 theFinalStatePhotons = new G4ParticleHPPhotonDist;
157 theFinalStatePhotons->InitMean(theData); 149 theFinalStatePhotons->InitMean(theData);
158 hasFSData = true; 150 hasFSData = true;
159 } 151 }
160 else if (dataType == 13) { 152 else if (dataType == 13) {
161 theFinalStatePhotons = new G4ParticleHPP 153 theFinalStatePhotons = new G4ParticleHPPhotonDist;
162 theFinalStatePhotons->InitPartials(theDa 154 theFinalStatePhotons->InitPartials(theData, theXsection);
163 hasFSData = true; 155 hasFSData = true;
164 } 156 }
165 else if (dataType == 14) { 157 else if (dataType == 14) {
166 theFinalStatePhotons->InitAngular(theDat 158 theFinalStatePhotons->InitAngular(theData);
167 hasFSData = true; 159 hasFSData = true;
168 } 160 }
169 else if (dataType == 15) { 161 else if (dataType == 15) {
170 theFinalStatePhotons->InitEnergies(theDa 162 theFinalStatePhotons->InitEnergies(theData);
171 hasFSData = true; 163 hasFSData = true;
172 } 164 }
173 else { 165 else {
174 G4ExceptionDescription ed; << 166 throw G4HadronicException(__FILE__, __LINE__,
175 ed << "Z=" << theBaseZ << " A=" << theBa << 167 "Data-type unknown to G4ParticleHPInelasticBaseFS");
176 << " projectile: " << theProjectile->GetPar <<
177 G4Exception("G4ParticleHPInelasticBaseFS <<
178 ed, "Data-type unknown"); <<
179 } 168 }
180 } 169 }
181 } 170 }
182 171
183 void G4ParticleHPInelasticBaseFS::BaseApply(co 172 void G4ParticleHPInelasticBaseFS::BaseApply(const G4HadProjectile& theTrack,
184 G4 << 173 G4ParticleDefinition** theDefs, G4int nDef)
185 G4 <<
186 { 174 {
187 // G4cout << "G4ParticleHPInelasticBaseFS::B <<
188 // prepare neutron 175 // prepare neutron
189 if (theResult.Get() == nullptr) theResult.Pu << 176 if (theResult.Get() == NULL) theResult.Put(new G4HadFinalState);
190 theResult.Get()->Clear(); 177 theResult.Get()->Clear();
191 G4double eKinetic = theTrack.GetKineticEnerg 178 G4double eKinetic = theTrack.GetKineticEnergy();
192 G4ReactionProduct incidReactionProduct(theTr << 179 const G4HadProjectile* hadProjectile = &theTrack;
193 incidReactionProduct.SetMomentum(theTrack.Ge << 180 G4ReactionProduct incidReactionProduct(
>> 181 const_cast<G4ParticleDefinition*>(hadProjectile->GetDefinition()));
>> 182 incidReactionProduct.SetMomentum(hadProjectile->Get4Momentum().vect());
194 incidReactionProduct.SetKineticEnergy(eKinet 183 incidReactionProduct.SetKineticEnergy(eKinetic);
195 184
196 // prepare target 185 // prepare target
197 G4double targetMass = G4NucleiProperties::Ge << 186 G4double targetMass;
198 /CLHEP::neutron_mass_c2; << 187 G4double eps = 0.0001;
199 /* << 188 targetMass = (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA + eps),
200 G4cout << " Target Z=" << theBaseZ << " A=" << 189 static_cast<G4int>(theBaseZ + eps)))
201 << " projectile: " << theTrack.GetDefinitio << 190 /
202 << " Ekin(MeV)=" << theTrack.GetKineticEner << 191 // theProjectile->GetPDGMass();
203 << " Mtar/Mn=" << targetMass << 192 G4Neutron::Neutron()->GetPDGMass();
204 << " hasFS=" << HasFSData() << G4endl; << 193
205 */ << 194 // give priority to ENDF vales for target mass
206 // give priority to ENDF vales for target ma << 195 if (theEnergyAngData != 0) {
207 if (theEnergyAngData != nullptr) { << 196 targetMass = theEnergyAngData->GetTargetMass();
208 G4double x = theEnergyAngData->GetTargetMa << 197 }
209 if(0.0 < x) { targetMass = x/CLHEP::neutro << 198 if (theAngularDistribution != 0) {
210 } << 199 targetMass = theAngularDistribution->GetTargetMass();
211 if (theAngularDistribution != nullptr) { <<
212 G4double x = theAngularDistribution->GetTa <<
213 if(0.0 < x) { targetMass = x/CLHEP::neutr <<
214 } 200 }
215 201
216 // 110512 TKDB ENDF-VII.0 21Sc45 has trouble 202 // 110512 TKDB ENDF-VII.0 21Sc45 has trouble in MF4MT22 (n,np) targetMass is not properly
217 // recorded. TKDB targetMass = 0; ENDF-VII.0 << 203 // recorded.
>> 204 if (targetMass == 0) {
>> 205 // G4cout << "TKDB targetMass = 0; ENDF-VII.0 21Sc45 has trouble in MF4MT22 (n,np) targetMass is
>> 206 // not properly recorded. This could be a similar situation." << G4endl; targetMass = (
>> 207 // G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps),
>> 208 // static_cast<G4int>(theBaseZ+eps))) / theProjectile->GetPDGMass();
>> 209 targetMass = (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA + eps),
>> 210 static_cast<G4int>(theBaseZ + eps)))
>> 211 / G4Neutron::Neutron()->GetPDGMass();
>> 212 }
218 213
219 G4Nucleus aNucleus; 214 G4Nucleus aNucleus;
220 G4ReactionProduct theTarget; 215 G4ReactionProduct theTarget;
221 216
222 G4ThreeVector neuVelo = 217 G4ThreeVector neuVelo =
223 incidReactionProduct.GetMomentum() / CLHEP << 218 (1. / G4Neutron::Neutron()->GetPDGMass()) * incidReactionProduct.GetMomentum();
224 <<
225 theTarget = 219 theTarget =
226 aNucleus.GetBiasedThermalNucleus(targetMas << 220 aNucleus.GetBiasedThermalNucleus(targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
227 theTrack.GetMaterial()->GetTemper << 221
228 theTarget.SetDefinition(ionTable->GetIon(the << 222 theTarget.SetDefinition(G4IonTable::GetIonTable()->GetIon(G4int(theBaseZ), G4int(theBaseA), 0.0));
229 223
230 // prepare energy in target rest frame 224 // prepare energy in target rest frame
231 G4ReactionProduct boosted; 225 G4ReactionProduct boosted;
232 boosted.Lorentz(incidReactionProduct, theTar 226 boosted.Lorentz(incidReactionProduct, theTarget);
233 eKinetic = boosted.GetKineticEnergy(); 227 eKinetic = boosted.GetKineticEnergy();
234 <<
235 G4double orgMomentum = boosted.GetMomentum() 228 G4double orgMomentum = boosted.GetMomentum().mag();
236 229
237 // Take N-body phase-space distribution, if 230 // Take N-body phase-space distribution, if no other data present.
238 if (!HasFSData()) // adding the residual is 231 if (!HasFSData()) // adding the residual is trivial here @@@
239 { 232 {
240 G4ParticleHPNBodyPhaseSpace thePhaseSpaceD 233 G4ParticleHPNBodyPhaseSpace thePhaseSpaceDistribution;
241 G4double aPhaseMass = 0; 234 G4double aPhaseMass = 0;
242 G4int ii; 235 G4int ii;
243 for (ii = 0; ii < nDef; ++ii) { 236 for (ii = 0; ii < nDef; ++ii) {
244 aPhaseMass += theDefs[ii]->GetPDGMass(); 237 aPhaseMass += theDefs[ii]->GetPDGMass();
245 } 238 }
246 239
247 //---------------------------------------- << 240 //----------------------------------------------------------------------------
248 if (std::abs(Qvalue) < CLHEP::keV) { << 241 if (Qvalue < 1. * CLHEP::keV && Qvalue > -1. * CLHEP::keV) {
249 // Not in the G4NDL lib or not calculate 242 // Not in the G4NDL lib or not calculated yet:
250 243
251 // Calculate residual: 244 // Calculate residual:
252 G4int ResidualA = theBaseA; 245 G4int ResidualA = theBaseA;
253 G4int ResidualZ = theBaseZ; 246 G4int ResidualZ = theBaseZ;
254 for (ii = 0; ii < nDef; ++ii) { 247 for (ii = 0; ii < nDef; ++ii) {
255 ResidualZ -= theDefs[ii]->GetAtomicNum 248 ResidualZ -= theDefs[ii]->GetAtomicNumber();
256 ResidualA -= theDefs[ii]->GetBaryonNum 249 ResidualA -= theDefs[ii]->GetBaryonNumber();
257 } 250 }
258 251
259 if (ResidualA > 0 && ResidualZ > 0) { 252 if (ResidualA > 0 && ResidualZ > 0) {
260 G4ParticleDefinition* resid = ionTable << 253 G4ParticleDefinition* resid = G4IonTable::GetIonTable()->GetIon(ResidualZ, ResidualA);
261 Qvalue = 254 Qvalue =
262 incidReactionProduct.GetMass() + the 255 incidReactionProduct.GetMass() + theTarget.GetMass() - aPhaseMass - resid->GetPDGMass();
263 } 256 }
264 257
265 if (std::abs(Qvalue) > 400 * CLHEP::MeV) << 258 if (Qvalue > 400 * CLHEP::MeV || Qvalue < -400 * CLHEP::MeV) {
266 // Then Q value is probably too large 259 // Then Q value is probably too large ...
267 Qvalue = 1.1 * CLHEP::keV; 260 Qvalue = 1.1 * CLHEP::keV;
268 } 261 }
269 } 262 }
270 //---------------------------------------- 263 //----------------------------------------------------------------------------
271 264
272 thePhaseSpaceDistribution.Init(aPhaseMass, 265 thePhaseSpaceDistribution.Init(aPhaseMass, nDef);
273 thePhaseSpaceDistribution.SetProjectileRP( 266 thePhaseSpaceDistribution.SetProjectileRP(&incidReactionProduct);
274 thePhaseSpaceDistribution.SetTarget(&theTa 267 thePhaseSpaceDistribution.SetTarget(&theTarget);
275 thePhaseSpaceDistribution.SetQValue(Qvalue 268 thePhaseSpaceDistribution.SetQValue(Qvalue);
276 269
277 for (ii = 0; ii < nDef; ++ii) { 270 for (ii = 0; ii < nDef; ++ii) {
278 G4double massCode = 1000. * std::abs(the 271 G4double massCode = 1000. * std::abs(theDefs[ii]->GetPDGCharge());
279 massCode += theDefs[ii]->GetBaryonNumber 272 massCode += theDefs[ii]->GetBaryonNumber();
280 G4double dummy = 0; 273 G4double dummy = 0;
281 G4ReactionProduct* aSec = thePhaseSpaceD 274 G4ReactionProduct* aSec = thePhaseSpaceDistribution.Sample(eKinetic, massCode, dummy);
282 aSec->Lorentz(*aSec, -1. * theTarget); 275 aSec->Lorentz(*aSec, -1. * theTarget);
283 auto aPart = new G4DynamicParticle(); << 276 G4DynamicParticle* aPart = new G4DynamicParticle();
284 aPart->SetDefinition(aSec->GetDefinition 277 aPart->SetDefinition(aSec->GetDefinition());
285 aPart->SetMomentum(aSec->GetMomentum()); 278 aPart->SetMomentum(aSec->GetMomentum());
286 delete aSec; 279 delete aSec;
287 theResult.Get()->AddSecondary(aPart, sec 280 theResult.Get()->AddSecondary(aPart, secID);
288 #ifdef G4VERBOSE << 281 #ifdef G4PHPDEBUG
289 if (fManager->GetDEBUG()) << 282 if (std::getenv("G4ParticleHPDebug"))
290 G4cout << this << " G4ParticleHPInelas 283 G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply NoFSData add secondary "
291 << aPart->GetParticleDefinition 284 << aPart->GetParticleDefinition()->GetParticleName()
292 << " E= " << aPart->GetKineticE 285 << " E= " << aPart->GetKineticEnergy() << " NSECO "
293 << theResult.Get()->GetNumberOf 286 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
294 #endif 287 #endif
295 } 288 }
296 <<
297 theResult.Get()->SetStatusChange(stopAndKi 289 theResult.Get()->SetStatusChange(stopAndKill);
298 // Final momentum check should be done bef 290 // Final momentum check should be done before return
299 const G4ParticleDefinition* targ_pd = ionT << 291 G4ParticleDefinition* targ_pd =
300 G4double mass = targ_pd->GetPDGMass(); << 292 G4IonTable::GetIonTable()->GetIon((G4int)theBaseZ, (G4int)theBaseA, 0.0);
301 G4LorentzVector targ_4p_lab(theTarget.GetM << 293 G4LorentzVector targ_4p_lab(
302 std::sqrt(mass*mass + theTarget.GetMo << 294 theTarget.GetMomentum(),
>> 295 std::sqrt(targ_pd->GetPDGMass() * targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2()));
303 G4LorentzVector proj_4p_lab = theTrack.Get 296 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
304 G4LorentzVector init_4p_lab = proj_4p_lab 297 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
305 adjust_final_state(init_4p_lab); 298 adjust_final_state(init_4p_lab);
>> 299
306 return; 300 return;
307 } 301 }
308 302
309 // set target and neutron in the relevant ex 303 // set target and neutron in the relevant exit channel
310 if (theAngularDistribution != nullptr) { << 304 if (theAngularDistribution != 0) {
311 theAngularDistribution->SetTarget(theTarge 305 theAngularDistribution->SetTarget(theTarget);
312 theAngularDistribution->SetProjectileRP(in 306 theAngularDistribution->SetProjectileRP(incidReactionProduct);
313 } 307 }
314 else if (theEnergyAngData != nullptr) { << 308 else if (theEnergyAngData != 0) {
315 theEnergyAngData->SetTarget(theTarget); 309 theEnergyAngData->SetTarget(theTarget);
316 theEnergyAngData->SetProjectileRP(incidRea 310 theEnergyAngData->SetProjectileRP(incidReactionProduct);
317 } 311 }
318 312
319 G4ReactionProductVector* tmpHadrons = nullpt << 313 G4ReactionProductVector* tmpHadrons = 0;
>> 314 #ifdef G4PHPDEBUG
>> 315 // To avoid compilation error around line 532.
>> 316 G4int ii(0);
>> 317 #endif
320 G4int dummy; 318 G4int dummy;
321 std::size_t i; 319 std::size_t i;
322 320
323 if (theEnergyAngData != nullptr) { << 321 if (theEnergyAngData != 0) {
324 tmpHadrons = theEnergyAngData->Sample(eKin 322 tmpHadrons = theEnergyAngData->Sample(eKinetic);
325 auto hproj = theTrack.GetDefinition(); <<
326 G4int projA = hproj->GetBaryonNumber(); <<
327 G4int projZ = G4lrint(hproj->GetPDGCharge( <<
328 323
329 if (!fManager->GetDoNotAdjustFinalState()) << 324 if (!G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) {
330 // Adjust A and Z in the case of miss mu 325 // Adjust A and Z in the case of miss much between selected data and target nucleus
331 if (tmpHadrons != nullptr) { 326 if (tmpHadrons != nullptr) {
332 G4int sumA = 0; 327 G4int sumA = 0;
333 G4int sumZ = 0; 328 G4int sumZ = 0;
334 G4int maxA = 0; 329 G4int maxA = 0;
335 G4int jAtMaxA = 0; 330 G4int jAtMaxA = 0;
336 for (G4int j = 0; j != (G4int)tmpHadro 331 for (G4int j = 0; j != (G4int)tmpHadrons->size(); ++j) {
337 auto had = ((*tmpHadrons)[j])->GetDefiniti << 332 // G4cout << __FILE__ << " " << __LINE__ << "th line:
338 if (had->GetBaryonNumber() > maxA) { << 333 // tmpHadrons->at(j)->GetDefinition()->GetParticleName() = " <<
339 maxA = had->GetBaryonNumber(); << 334 // tmpHadrons->at(j)->GetDefinition()->GetParticleName() << G4endl;
>> 335 if (tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber() > maxA) {
>> 336 maxA = tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber();
340 jAtMaxA = j; 337 jAtMaxA = j;
341 } 338 }
342 sumA += had->GetBaryonNumber(); << 339 sumA += tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber();
343 sumZ += G4lrint(had->GetPDGCharge()/ << 340 sumZ += G4int(tmpHadrons->at(j)->GetDefinition()->GetPDGCharge() + eps);
344 } 341 }
345 G4int dA = theBaseA + projA - sumA; << 342 G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
346 G4int dZ = theBaseZ + projZ - sumZ; << 343 G4int dZ =
>> 344 (G4int)theBaseZ + G4int(hadProjectile->GetDefinition()->GetPDGCharge() + eps) - sumZ;
347 if (dA < 0 || dZ < 0) { 345 if (dA < 0 || dZ < 0) {
348 auto p = ((*tmpHadrons)[jAtMaxA])->G << 346 G4int newA = tmpHadrons->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA;
349 G4int newA = p->GetBaryonNumber() + << 347 G4int newZ = G4int(tmpHadrons->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps) + dZ;
350 G4int newZ = G4lrint(p->GetPDGCharge <<
351 if (newA > newZ && newZ > 0) { 348 if (newA > newZ && newZ > 0) {
352 G4ParticleDefinition* pd = ionTabl << 349 G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon(newZ, newA);
353 ((*tmpHadrons)[jAtMaxA])->SetDefin << 350 tmpHadrons->at(jAtMaxA)->SetDefinition(pd);
354 } 351 }
355 } 352 }
356 } 353 }
357 } 354 }
358 } 355 }
359 else if (theAngularDistribution != nullptr) << 356 else if (theAngularDistribution != 0) {
360 << 357 G4bool* Done = new G4bool[nDef];
361 auto Done = new G4bool[nDef]; <<
362 G4int i0; 358 G4int i0;
363 for (i0 = 0; i0 < nDef; ++i0) 359 for (i0 = 0; i0 < nDef; ++i0)
364 Done[i0] = false; 360 Done[i0] = false;
365 361
366 tmpHadrons = new G4ReactionProductVector; 362 tmpHadrons = new G4ReactionProductVector;
367 G4ReactionProduct* aHadron; 363 G4ReactionProduct* aHadron;
368 G4double localMass = G4NucleiProperties::G << 364 G4double localMass = (G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA + eps),
>> 365 static_cast<G4int>(theBaseZ + eps)));
369 G4ThreeVector bufferedDirection(0, 0, 0); 366 G4ThreeVector bufferedDirection(0, 0, 0);
370 for (i0 = 0; i0 < nDef; ++i0) { 367 for (i0 = 0; i0 < nDef; ++i0) {
371 if (!Done[i0]) { 368 if (!Done[i0]) {
372 aHadron = new G4ReactionProduct; 369 aHadron = new G4ReactionProduct;
373 if (theEnergyDistribution != nullptr) << 370 if (theEnergyDistribution != 0) {
374 aHadron->SetDefinition(theDefs[i0]); 371 aHadron->SetDefinition(theDefs[i0]);
375 aHadron->SetKineticEnergy(theEnergyD 372 aHadron->SetKineticEnergy(theEnergyDistribution->Sample(eKinetic, dummy));
376 } 373 }
377 else if (nDef == 1) { 374 else if (nDef == 1) {
378 aHadron->SetDefinition(theDefs[i0]); 375 aHadron->SetDefinition(theDefs[i0]);
379 aHadron->SetKineticEnergy(eKinetic); 376 aHadron->SetKineticEnergy(eKinetic);
380 } 377 }
381 else if (nDef == 2) { 378 else if (nDef == 2) {
382 aHadron->SetDefinition(theDefs[i0]); 379 aHadron->SetDefinition(theDefs[i0]);
383 aHadron->SetKineticEnergy(50 * CLHEP 380 aHadron->SetKineticEnergy(50 * CLHEP::MeV);
384 } 381 }
385 else { 382 else {
386 throw G4HadronicException( 383 throw G4HadronicException(
387 __FILE__, __LINE__, 384 __FILE__, __LINE__,
388 "No energy distribution to sample 385 "No energy distribution to sample from in InelasticBaseFS::BaseApply");
389 } 386 }
390 theAngularDistribution->SampleAndUpdat 387 theAngularDistribution->SampleAndUpdate(*aHadron);
391 if (theEnergyDistribution == nullptr & << 388 if (theEnergyDistribution == 0 && nDef == 2) {
392 if (i0 == 0) { 389 if (i0 == 0) {
393 G4double mass1 = theDefs[0]->GetPD 390 G4double mass1 = theDefs[0]->GetPDGMass();
394 G4double mass2 = theDefs[1]->GetPD 391 G4double mass2 = theDefs[1]->GetPDGMass();
395 G4double massn = theProjectile->Ge 392 G4double massn = theProjectile->GetPDGMass();
396 G4int z1 = theBaseZ - << 393 G4int z1 = static_cast<G4int>(theBaseZ + eps - theDefs[0]->GetPDGCharge()
397 G4lrint((theDefs[0]->GetPDGCharge() + << 394 - theDefs[1]->GetPDGCharge());
398 G4int a1 = theBaseA - theDefs[0]-> << 395 G4int a1 = static_cast<G4int>(theBaseA + eps) - theDefs[0]->GetBaryonNumber()
>> 396 - theDefs[1]->GetBaryonNumber();
399 G4double concreteMass = G4NucleiPr 397 G4double concreteMass = G4NucleiProperties::GetNuclearMass(a1, z1);
400 G4double availableEnergy = eKineti << 398 G4double availableEnergy = eKinetic + massn + localMass - mass1 - mass2 - concreteMass;
401 - concreteMass; <<
402 // available kinetic energy in CMS 399 // available kinetic energy in CMS (non relativistic)
403 G4double emin = availableEnergy + << 400 G4double emin =
>> 401 availableEnergy + mass1 + mass2
404 - std::sqrt((mass1 + mass2) * (m 402 - std::sqrt((mass1 + mass2) * (mass1 + mass2) + orgMomentum * orgMomentum);
405 G4double p1 = std::sqrt(2. * mass2 403 G4double p1 = std::sqrt(2. * mass2 * emin);
406 bufferedDirection = p1 * aHadron-> 404 bufferedDirection = p1 * aHadron->GetMomentum().unit();
407 #ifdef G4PHPDEBUG 405 #ifdef G4PHPDEBUG
408 if (fManager->GetDEBUG()) // @@@@ << 406 if (std::getenv("G4ParticleHPDebug")) // @@@@@ verify the nucleon counting...
409 { 407 {
410 G4cout << "G4ParticleHPInelastic << 408 G4cout << "G4ParticleHPInelasticBaseFS " << z1 << " " << theBaseZ << " " << a1 << " "
411 << " " << a1 << " " << theBaseA << " << 409 << theBaseA << " " << availableEnergy << " " << emin << G4endl;
412 << " " << emin << G4endl; <<
413 } 410 }
414 #endif 411 #endif
415 } 412 }
416 else { 413 else {
417 bufferedDirection = -bufferedDirec 414 bufferedDirection = -bufferedDirection;
418 } 415 }
419 // boost from cms to lab 416 // boost from cms to lab
>> 417 #ifdef G4PHPDEBUG
>> 418 if (std::getenv("G4ParticleHPDebug")) {
>> 419 G4cout << " G4ParticleHPInelasticBaseFS " << bufferedDirection.mag2() << G4endl;
>> 420 }
>> 421 #endif
420 aHadron->SetTotalEnergy( 422 aHadron->SetTotalEnergy(
421 std::sqrt(aHadron->GetMass() * aHa 423 std::sqrt(aHadron->GetMass() * aHadron->GetMass() + bufferedDirection.mag2()));
422 aHadron->SetMomentum(bufferedDirecti 424 aHadron->SetMomentum(bufferedDirection);
423 aHadron->Lorentz(*aHadron, -1. * (th 425 aHadron->Lorentz(*aHadron, -1. * (theTarget + incidReactionProduct));
424 #ifdef G4PHPDEBUG 426 #ifdef G4PHPDEBUG
425 if (fManager->GetDEBUG()) { << 427 if (std::getenv("G4ParticleHPDebug")) {
426 G4cout << " G4ParticleHPInelasticB << 428 G4cout << " G4ParticleHPInelasticBaseFS " << aHadron->GetTotalEnergy() << " "
427 << " P=" << aHadron->GetMomentum() << 429 << aHadron->GetMomentum() << G4endl;
428 << " bufDir^2=" << bufferedDirection.ma <<
429 << G4endl; <<
430 } 430 }
431 #endif 431 #endif
432 } 432 }
433 tmpHadrons->push_back(aHadron); 433 tmpHadrons->push_back(aHadron);
434 #ifdef G4PHPDEBUG 434 #ifdef G4PHPDEBUG
435 if (fManager->GetDEBUG()) << 435 if (std::getenv("G4ParticleHPDebug"))
436 G4cout << " G4ParticleHPInelasticBas 436 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply FSData add secondary "
437 << aHadron->GetDefinition()-> 437 << aHadron->GetDefinition()->GetParticleName()
438 << " E= " << aHadron->GetKine 438 << " E= " << aHadron->GetKineticEnergy() << G4endl;
439 #endif 439 #endif
440 } 440 }
441 } 441 }
442 delete[] Done; 442 delete[] Done;
443 } 443 }
444 else { 444 else {
445 throw G4HadronicException(__FILE__, __LINE << 445 throw G4HadronicException(__FILE__, __LINE__, "No data to create the neutrons in NInelasticFS");
446 } 446 }
447 447
448 G4ReactionProductVector* thePhotons = nullpt 448 G4ReactionProductVector* thePhotons = nullptr;
449 if (theFinalStatePhotons != nullptr) { << 449 if (theFinalStatePhotons != 0) {
450 // the photon distributions are in the Nuc 450 // the photon distributions are in the Nucleus rest frame.
451 G4ReactionProduct boosted_tmp; 451 G4ReactionProduct boosted_tmp;
452 boosted_tmp.Lorentz(incidReactionProduct, 452 boosted_tmp.Lorentz(incidReactionProduct, theTarget);
453 G4double anEnergy = boosted_tmp.GetKinetic 453 G4double anEnergy = boosted_tmp.GetKineticEnergy();
454 thePhotons = theFinalStatePhotons->GetPhot 454 thePhotons = theFinalStatePhotons->GetPhotons(anEnergy);
455 if (thePhotons != nullptr) { << 455 if (thePhotons != 0) {
456 for (i = 0; i < thePhotons->size(); ++i) 456 for (i = 0; i < thePhotons->size(); ++i) {
457 // back to lab 457 // back to lab
458 thePhotons->operator[](i)->Lorentz(*(t 458 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1. * theTarget);
459 } 459 }
460 } 460 }
461 } 461 }
462 else if (theEnergyAngData != nullptr) { << 462 else if (theEnergyAngData != 0) {
463 // PA130927: do not create photons to adju 463 // PA130927: do not create photons to adjust binding energy
464 G4bool bAdjustPhotons{true}; << 464 G4bool bAdjustPhotons = true;
465 #ifdef PHP_AS_HP 465 #ifdef PHP_AS_HP
466 bAdjustPhotons = true; 466 bAdjustPhotons = true;
467 #else 467 #else
468 if (fManager->GetDoNotAdjustFinalState()) << 468 if (G4ParticleHPManager::GetInstance()->GetDoNotAdjustFinalState()) bAdjustPhotons = false;
469 #endif 469 #endif
470 470
471 if (bAdjustPhotons) { 471 if (bAdjustPhotons) {
472 G4double theGammaEnergy = theEnergyAngDa 472 G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
473 G4double anEnergy = boosted.GetKineticEn 473 G4double anEnergy = boosted.GetKineticEnergy();
474 theGammaEnergy = anEnergy - theGammaEner 474 theGammaEnergy = anEnergy - theGammaEnergy;
475 theGammaEnergy += theNuclearMassDifferen 475 theGammaEnergy += theNuclearMassDifference;
476 G4double eBindProducts = 0; 476 G4double eBindProducts = 0;
477 G4double eBindN = 0; 477 G4double eBindN = 0;
478 G4double eBindP = 0; 478 G4double eBindP = 0;
479 G4double eBindD = G4NucleiProperties::Ge 479 G4double eBindD = G4NucleiProperties::GetBindingEnergy(2, 1);
480 G4double eBindT = G4NucleiProperties::Ge 480 G4double eBindT = G4NucleiProperties::GetBindingEnergy(3, 1);
481 G4double eBindHe3 = G4NucleiProperties:: 481 G4double eBindHe3 = G4NucleiProperties::GetBindingEnergy(3, 2);
482 G4double eBindA = G4NucleiProperties::Ge 482 G4double eBindA = G4NucleiProperties::GetBindingEnergy(4, 2);
483 G4int ia = 0; 483 G4int ia = 0;
484 for (auto const & had : *tmpHadrons) { << 484 for (i = 0; i < tmpHadrons->size(); i++) {
485 if (had->GetDefinition() == G4Neutron: << 485 if (tmpHadrons->operator[](i)->GetDefinition() == G4Neutron::Neutron()) {
486 eBindProducts += eBindN; 486 eBindProducts += eBindN;
487 } 487 }
488 else if (had->GetDefinition() == G4Pro << 488 else if (tmpHadrons->operator[](i)->GetDefinition() == G4Proton::Proton()) {
489 eBindProducts += eBindP; 489 eBindProducts += eBindP;
490 } 490 }
491 else if (had->GetDefinition() == G4Deu << 491 else if (tmpHadrons->operator[](i)->GetDefinition() == G4Deuteron::Deuteron()) {
492 eBindProducts += eBindD; 492 eBindProducts += eBindD;
493 } 493 }
494 else if (had->GetDefinition() == G4Tri << 494 else if (tmpHadrons->operator[](i)->GetDefinition() == G4Triton::Triton()) {
495 eBindProducts += eBindT; 495 eBindProducts += eBindT;
496 } 496 }
497 else if (had->GetDefinition() == G4He3 << 497 else if (tmpHadrons->operator[](i)->GetDefinition() == G4He3::He3()) {
498 eBindProducts += eBindHe3; 498 eBindProducts += eBindHe3;
499 } 499 }
500 else if (had->GetDefinition() == G4Alp << 500 else if (tmpHadrons->operator[](i)->GetDefinition() == G4Alpha::Alpha()) {
501 eBindProducts += eBindA; 501 eBindProducts += eBindA;
502 ia++; 502 ia++;
503 } 503 }
504 } 504 }
>> 505
505 theGammaEnergy += eBindProducts; 506 theGammaEnergy += eBindProducts;
506 507
507 #ifdef G4PHPDEBUG 508 #ifdef G4PHPDEBUG
508 if (fManager->GetDEBUG()) << 509 if (std::getenv("G4ParticleHPDebug"))
509 G4cout << " G4ParticleHPInelasticBaseF 510 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply gamma Energy " << theGammaEnergy
510 << " eBindProducts " << eBindPr 511 << " eBindProducts " << eBindProducts << G4endl;
511 #endif 512 #endif
512 513
513 // Special treatment for Be9 + n -> 2n + 514 // Special treatment for Be9 + n -> 2n + Be8 -> 2n + a + a
514 if (theBaseZ == 4 && theBaseA == 9) { << 515 if ((G4int)(theBaseZ + eps) == 4 && (G4int)(theBaseA + eps) == 9) {
515 // This only valid for G4NDL3.13,,, 516 // This only valid for G4NDL3.13,,,
516 if (ia == 2 && std::abs(G4NucleiProper << 517 if (std::abs(theNuclearMassDifference
517 G4NucleiProperties::GetBindingEnergy(9 << 518 - (G4NucleiProperties::GetBindingEnergy(8, 4)
518 theNuclearMassDifference) < CLHEP::keV << 519 - G4NucleiProperties::GetBindingEnergy(9, 4)))
>> 520 < 1 * CLHEP::keV
>> 521 && ia == 2)
519 { 522 {
520 theGammaEnergy -= (2 * eBindA); 523 theGammaEnergy -= (2 * eBindA);
521 } 524 }
522 } 525 }
523 526
524 if (theGammaEnergy > 0.0) { 527 if (theGammaEnergy > 0.0) {
525 for (G4int iLevel = theGammas.GetNumberOfLev 528 for (G4int iLevel = theGammas.GetNumberOfLevels() - 1; iLevel > 0; --iLevel) {
526 G4double e = theGammas.GetLevelEnerg 529 G4double e = theGammas.GetLevelEnergy(iLevel);
527 if (e < theGammaEnergy) { 530 if (e < theGammaEnergy) {
528 thePhotons = theGammas.GetDecayGam 531 thePhotons = theGammas.GetDecayGammas(iLevel);
529 theGammaEnergy -= e; 532 theGammaEnergy -= e;
530 break; 533 break;
531 } 534 }
532 } 535 }
533 } 536 }
534 } 537 }
535 } 538 }
536 // fill the result 539 // fill the result
537 std::size_t nSecondaries = tmpHadrons->size( 540 std::size_t nSecondaries = tmpHadrons->size();
538 std::size_t nPhotons = 0; 541 std::size_t nPhotons = 0;
539 if (thePhotons != nullptr) { << 542 if (thePhotons != 0) {
540 nPhotons = thePhotons->size(); 543 nPhotons = thePhotons->size();
541 } 544 }
542 nSecondaries += nPhotons; 545 nSecondaries += nPhotons;
543 <<
544 G4DynamicParticle* theSec; 546 G4DynamicParticle* theSec;
545 #ifdef G4PHPDEBUG 547 #ifdef G4PHPDEBUG
546 if (fManager->GetDEBUG()) << 548 if (std::getenv("G4ParticleHPDebug"))
547 G4cout << " G4ParticleHPInelasticBaseFS::B 549 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N hadrons " << nSecondaries - nPhotons
548 << G4endl; 550 << G4endl;
549 #endif 551 #endif
550 552
551 for (i = 0; i < nSecondaries - nPhotons; ++i 553 for (i = 0; i < nSecondaries - nPhotons; ++i) {
552 theSec = new G4DynamicParticle; 554 theSec = new G4DynamicParticle;
553 theSec->SetDefinition(tmpHadrons->operator 555 theSec->SetDefinition(tmpHadrons->operator[](i)->GetDefinition());
554 theSec->SetMomentum(tmpHadrons->operator[] 556 theSec->SetMomentum(tmpHadrons->operator[](i)->GetMomentum());
555 theResult.Get()->AddSecondary(theSec, secI 557 theResult.Get()->AddSecondary(theSec, secID);
556 #ifdef G4PHPDEBUG 558 #ifdef G4PHPDEBUG
557 if (fManager->GetDEBUG()) << 559 if (std::getenv("G4ParticleHPDebug"))
558 G4cout << " G4ParticleHPInelasticBaseFS: << 560 G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply add secondary2 "
559 << theSec->GetParticleDefinition( 561 << theSec->GetParticleDefinition()->GetParticleName()
560 << " E=" << theSec->GetKineticEne << 562 << " E= " << theSec->GetKineticEnergy() << " NSECO "
561 << theResult.Get()->GetNumberOfSe 563 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
562 #endif 564 #endif
563 delete (*tmpHadrons)[i]; << 565 delete tmpHadrons->operator[](i);
564 } 566 }
565 #ifdef G4PHPDEBUG 567 #ifdef G4PHPDEBUG
566 if (fManager->GetDEBUG()) << 568 if (std::getenv("G4ParticleHPDebug"))
567 G4cout << " G4ParticleHPInelasticBaseFS::B 569 G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N photons " << nPhotons << G4endl;
568 #endif 570 #endif
569 if (thePhotons != nullptr) { << 571 if (thePhotons != 0) {
570 for (i = 0; i < nPhotons; ++i) { 572 for (i = 0; i < nPhotons; ++i) {
571 theSec = new G4DynamicParticle; 573 theSec = new G4DynamicParticle;
572 theSec->SetDefinition((*thePhotons)[i]-> << 574 theSec->SetDefinition(thePhotons->operator[](i)->GetDefinition());
573 theSec->SetMomentum((*thePhotons)[i]->Ge << 575 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
574 theResult.Get()->AddSecondary(theSec, se 576 theResult.Get()->AddSecondary(theSec, secID);
575 #ifdef G4PHPDEBUG 577 #ifdef G4PHPDEBUG
576 if (fManager->GetDEBUG()) << 578 if (std::getenv("G4ParticleHPDebug"))
577 G4cout << " G4ParticleHPInelasticBaseF << 579 G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply add secondary3 "
578 << theSec->GetParticleDefinitio 580 << theSec->GetParticleDefinition()->GetParticleName()
579 << " E=" << theSec->GetKineticE << 581 << " E= " << theSec->GetKineticEnergy() << " NSECO "
580 << theResult.Get()->GetNumberOf 582 << theResult.Get()->GetNumberOfSecondaries() << G4endl;
581 #endif 583 #endif
582 delete thePhotons->operator[](i); 584 delete thePhotons->operator[](i);
583 } 585 }
584 } 586 }
585 587
586 // some garbage collection 588 // some garbage collection
587 delete thePhotons; 589 delete thePhotons;
588 delete tmpHadrons; 590 delete tmpHadrons;
589 591
590 G4ParticleDefinition* targ_pd = ionTable->Ge << 592 G4ParticleDefinition* targ_pd =
>> 593 G4IonTable::GetIonTable()->GetIon((G4int)theBaseZ, (G4int)theBaseA, 0.0);
591 G4LorentzVector targ_4p_lab( 594 G4LorentzVector targ_4p_lab(
592 theTarget.GetMomentum(), 595 theTarget.GetMomentum(),
593 std::sqrt(targ_pd->GetPDGMass() * targ_pd- 596 std::sqrt(targ_pd->GetPDGMass() * targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2()));
594 G4LorentzVector proj_4p_lab = theTrack.Get4M 597 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
595 G4LorentzVector init_4p_lab = proj_4p_lab + 598 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
596 599
597 // if data in MF=6 format (no correlated par 600 // if data in MF=6 format (no correlated particle emission), then adjust_final_state can give
598 // severe errors: 601 // severe errors:
599 if (theEnergyAngData == nullptr) { << 602 if (theEnergyAngData == 0) {
600 adjust_final_state(init_4p_lab); 603 adjust_final_state(init_4p_lab);
601 } 604 }
602 605
603 // clean up the primary neutron 606 // clean up the primary neutron
604 theResult.Get()->SetStatusChange(stopAndKill 607 theResult.Get()->SetStatusChange(stopAndKill);
605 } 608 }
606 609