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