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