| 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 if ( tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber() > maxA ) {
336 for (G4int j = 0; j != (G4int)tmpHadro << 310 maxA = tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber();
337 auto had = ((*tmpHadrons)[j])->GetDefiniti << 311 jAtMaxA = j;
338 if (had->GetBaryonNumber() > maxA) { << 312 }
339 maxA = had->GetBaryonNumber(); << 313 sumA += tmpHadrons->at(j)->GetDefinition()->GetBaryonNumber();
340 jAtMaxA = j; << 314 sumZ += G4int( tmpHadrons->at(j)->GetDefinition()->GetPDGCharge() + eps );
341 } << 315 }
342 sumA += had->GetBaryonNumber(); << 316 G4int dA = (G4int)theBaseA + hadProjectile->GetDefinition()->GetBaryonNumber() - sumA;
343 sumZ += G4lrint(had->GetPDGCharge()/ << 317 G4int dZ = (G4int)theBaseZ + G4int( hadProjectile->GetDefinition()->GetPDGCharge() + eps ) - sumZ;
344 } << 318 if ( dA < 0 || dZ < 0 ) {
345 G4int dA = theBaseA + projA - sumA; << 319 G4int newA = tmpHadrons->at(jAtMaxA)->GetDefinition()->GetBaryonNumber() + dA ;
346 G4int dZ = theBaseZ + projZ - sumZ; << 320 G4int newZ = G4int( tmpHadrons->at(jAtMaxA)->GetDefinition()->GetPDGCharge() + eps ) + dZ;
347 if (dA < 0 || dZ < 0) { << 321 G4ParticleDefinition* pd = G4IonTable::GetIonTable()->GetIon ( newZ , newA );
348 auto p = ((*tmpHadrons)[jAtMaxA])->G << 322 tmpHadrons->at( jAtMaxA )->SetDefinition( pd );
349 G4int newA = p->GetBaryonNumber() + << 323 }
350 G4int newZ = G4lrint(p->GetPDGCharge <<
351 if (newA > newZ && newZ > 0) { <<
352 G4ParticleDefinition* pd = ionTabl <<
353 ((*tmpHadrons)[jAtMaxA])->SetDefin <<
354 } <<
355 } <<
356 } 324 }
>> 325 //141017 Fix END
357 } 326 }
358 } 327 }
359 else if (theAngularDistribution != nullptr) << 328 else if(theAngularDistribution!= 0)
360 << 329 {
361 auto Done = new G4bool[nDef]; << 330 G4bool * Done = new G4bool[nDef];
362 G4int i0; 331 G4int i0;
363 for (i0 = 0; i0 < nDef; ++i0) << 332 for(i0=0; i0<nDef; i0++) Done[i0] = false;
364 Done[i0] = false; <<
365 333
366 tmpHadrons = new G4ReactionProductVector; << 334 //Following lines are commented out to fix coverity defeat 58622
367 G4ReactionProduct* aHadron; << 335 // if(tmpHadrons == 0)
368 G4double localMass = G4NucleiProperties::G << 336 // {
369 G4ThreeVector bufferedDirection(0, 0, 0); << 337 tmpHadrons = new G4ReactionProductVector;
370 for (i0 = 0; i0 < nDef; ++i0) { << 338 // }
371 if (!Done[i0]) { << 339 // else
>> 340 // {
>> 341 // for(i=0; i<tmpHadrons->size(); i++)
>> 342 // {
>> 343 // for(ii=0; ii<nDef; ii++)
>> 344 // if(!Done[ii] && tmpHadrons->operator[](i)->GetDefinition() == theDefs[ii])
>> 345 // Done[ii] = true;
>> 346 // }
>> 347 #ifdef G4PHPDEBUG
>> 348 // if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply secondary previously added " << tmpHadrons->operator[](i)->GetDefinition()->GetParticleName() << " E= " << tmpHadrons->operator[](i)->GetKineticEnergy() << G4endl;
>> 349 #endif
>> 350 // }
>> 351 G4ReactionProduct * aHadron;
>> 352 G4double localMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps)));
>> 353 G4ThreeVector bufferedDirection(0,0,0);
>> 354 for(i0=0; i0<nDef; i0++)
>> 355 {
>> 356 if(!Done[i0])
>> 357 {
372 aHadron = new G4ReactionProduct; 358 aHadron = new G4ReactionProduct;
373 if (theEnergyDistribution != nullptr) << 359 if(theEnergyDistribution!=0)
374 aHadron->SetDefinition(theDefs[i0]); << 360 {
375 aHadron->SetKineticEnergy(theEnergyD << 361 aHadron->SetDefinition(theDefs[i0]);
376 } << 362 aHadron->SetKineticEnergy(theEnergyDistribution->Sample(eKinetic, dummy));
377 else if (nDef == 1) { << 363 }
378 aHadron->SetDefinition(theDefs[i0]); << 364 else if(nDef == 1)
379 aHadron->SetKineticEnergy(eKinetic); << 365 {
380 } << 366 aHadron->SetDefinition(theDefs[i0]);
381 else if (nDef == 2) { << 367 aHadron->SetKineticEnergy(eKinetic);
382 aHadron->SetDefinition(theDefs[i0]); << 368 }
383 aHadron->SetKineticEnergy(50 * CLHEP << 369 else if(nDef == 2)
384 } << 370 {
385 else { << 371 aHadron->SetDefinition(theDefs[i0]);
386 throw G4HadronicException( << 372 aHadron->SetKineticEnergy(50*CLHEP::MeV);
387 __FILE__, __LINE__, << 373 }
388 "No energy distribution to sample << 374 else
389 } << 375 {
390 theAngularDistribution->SampleAndUpdat << 376 throw G4HadronicException(__FILE__, __LINE__, "No energy distribution to sample from in InelasticBaseFS::BaseApply");
391 if (theEnergyDistribution == nullptr & << 377 }
392 if (i0 == 0) { << 378 theAngularDistribution->SampleAndUpdate(*aHadron);
393 G4double mass1 = theDefs[0]->GetPD << 379 if(theEnergyDistribution==0 && nDef == 2)
394 G4double mass2 = theDefs[1]->GetPD << 380 {
395 G4double massn = theProjectile->Ge << 381 if(i0==0)
396 G4int z1 = theBaseZ - << 382 {
397 G4lrint((theDefs[0]->GetPDGCharge() + << 383 G4double mass1 = theDefs[0]->GetPDGMass();
398 G4int a1 = theBaseA - theDefs[0]-> << 384 G4double mass2 = theDefs[1]->GetPDGMass();
399 G4double concreteMass = G4NucleiPr << 385 G4double massn = theProjectile->GetPDGMass();
400 G4double availableEnergy = eKineti << 386 G4int z1 = static_cast<G4int>(theBaseZ+eps-theDefs[0]->GetPDGCharge()-theDefs[1]->GetPDGCharge());
401 - concreteMass; << 387 G4int a1 = static_cast<G4int>(theBaseA+eps)-theDefs[0]->GetBaryonNumber()-theDefs[1]->GetBaryonNumber();
402 // available kinetic energy in CMS << 388 G4double concreteMass = G4NucleiProperties::GetNuclearMass(a1, z1);
403 G4double emin = availableEnergy + << 389 G4double availableEnergy = eKinetic+massn+localMass-mass1-mass2-concreteMass;
404 - std::sqrt((mass1 + mass2) * (m << 390 // available kinetic energy in CMS (non relativistic)
405 G4double p1 = std::sqrt(2. * mass2 << 391 G4double emin = availableEnergy+mass1+mass2 - std::sqrt((mass1+mass2)*(mass1+mass2)+orgMomentum*orgMomentum);
406 bufferedDirection = p1 * aHadron-> << 392 G4double p1=std::sqrt(2.*mass2*emin);
407 #ifdef G4PHPDEBUG << 393 bufferedDirection = p1*aHadron->GetMomentum().unit();
408 if (fManager->GetDEBUG()) // @@@@ << 394 #ifdef G4PHPDEBUG
409 { << 395 if(getenv("G4ParticleHPDebug")) // @@@@@ verify the nucleon counting...
410 G4cout << "G4ParticleHPInelastic << 396 {
411 << " " << a1 << " " << theBaseA << " << 397 G4cout << "G4ParticleHPInelasticBaseFS "<<z1<<" "<<theBaseZ<<" "<<a1<<" "<<theBaseA<<" "<<availableEnergy<<" "
412 << " " << emin << G4endl; << 398 << emin<<G4endl;
413 } 399 }
414 #endif 400 #endif
415 } << 401 }
416 else { << 402 else
417 bufferedDirection = -bufferedDirec << 403 {
418 } << 404 bufferedDirection = -bufferedDirection;
419 // boost from cms to lab << 405 }
420 aHadron->SetTotalEnergy( << 406 // boost from cms to lab
421 std::sqrt(aHadron->GetMass() * aHa << 407 #ifdef G4PHPDEBUG
422 aHadron->SetMomentum(bufferedDirecti << 408 if(getenv("G4ParticleHPDebug"))
423 aHadron->Lorentz(*aHadron, -1. * (th << 409 {
424 #ifdef G4PHPDEBUG << 410 G4cout << " G4ParticleHPInelasticBaseFS "<<bufferedDirection.mag2()<<G4endl;
425 if (fManager->GetDEBUG()) { << 411 }
426 G4cout << " G4ParticleHPInelasticB << 412 #endif
427 << " P=" << aHadron->GetMomentum() << 413 aHadron->SetTotalEnergy( std::sqrt(aHadron->GetMass()*aHadron->GetMass()
428 << " bufDir^2=" << bufferedDirection.ma << 414 +bufferedDirection.mag2()) );
429 << G4endl; << 415 aHadron->SetMomentum(bufferedDirection);
430 } << 416 aHadron->Lorentz(*aHadron, -1.*(theTarget+incidReactionProduct));
431 #endif << 417 #ifdef G4PHPDEBUG
432 } << 418 if(getenv("G4ParticleHPDebug"))
433 tmpHadrons->push_back(aHadron); << 419 {
434 #ifdef G4PHPDEBUG << 420 G4cout << " G4ParticleHPInelasticBaseFS "<<aHadron->GetTotalEnergy()<<" "<<aHadron->GetMomentum()<<G4endl;
435 if (fManager->GetDEBUG()) << 421 }
436 G4cout << " G4ParticleHPInelasticBas << 422 #endif
437 << aHadron->GetDefinition()-> << 423 }
438 << " E= " << aHadron->GetKine << 424 tmpHadrons->push_back(aHadron);
>> 425 #ifdef G4PHPDEBUG
>> 426 if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply FSData add secondary " << aHadron->GetDefinition()->GetParticleName() << " E= " << aHadron->GetKineticEnergy() << G4endl;
439 #endif 427 #endif
440 } 428 }
441 } 429 }
442 delete[] Done; << 430 delete [] Done;
443 } 431 }
444 else { << 432 else
445 throw G4HadronicException(__FILE__, __LINE << 433 {
>> 434 throw G4HadronicException(__FILE__, __LINE__, "No data to create the neutrons in NInelasticFS");
446 } 435 }
447 436
448 G4ReactionProductVector* thePhotons = nullpt << 437 G4ReactionProductVector * thePhotons = 0;
449 if (theFinalStatePhotons != nullptr) { << 438 if(theFinalStatePhotons!=0)
>> 439 {
450 // the photon distributions are in the Nuc 440 // the photon distributions are in the Nucleus rest frame.
451 G4ReactionProduct boosted_tmp; 441 G4ReactionProduct boosted_tmp;
452 boosted_tmp.Lorentz(incidReactionProduct, 442 boosted_tmp.Lorentz(incidReactionProduct, theTarget);
453 G4double anEnergy = boosted_tmp.GetKinetic 443 G4double anEnergy = boosted_tmp.GetKineticEnergy();
454 thePhotons = theFinalStatePhotons->GetPhot 444 thePhotons = theFinalStatePhotons->GetPhotons(anEnergy);
455 if (thePhotons != nullptr) { << 445 if(thePhotons!=0)
456 for (i = 0; i < thePhotons->size(); ++i) << 446 {
>> 447 for(i=0; i<thePhotons->size(); i++)
>> 448 {
457 // back to lab 449 // back to lab
458 thePhotons->operator[](i)->Lorentz(*(t << 450 thePhotons->operator[](i)->Lorentz(*(thePhotons->operator[](i)), -1.*theTarget);
459 } 451 }
460 } 452 }
461 } 453 }
462 else if (theEnergyAngData != nullptr) { << 454 else if(theEnergyAngData!=0)
463 // PA130927: do not create photons to adju << 455 {
464 G4bool bAdjustPhotons{true}; << 456
465 #ifdef PHP_AS_HP << 457 // PA130927: do not create photons to adjust binding energy
466 bAdjustPhotons = true; << 458 G4bool bAdjustPhotons = true;
>> 459 #ifdef PHP_AS_HP
>> 460 bAdjustPhotons = true;
467 #else 461 #else
468 if (fManager->GetDoNotAdjustFinalState()) << 462 if ( getenv( "G4PHP_DO_NOT_ADJUST_FINAL_STATE" ) ) bAdjustPhotons = false;
469 #endif 463 #endif
470 << 464
471 if (bAdjustPhotons) { << 465 if( bAdjustPhotons ) {
472 G4double theGammaEnergy = theEnergyAngDa 466 G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
473 G4double anEnergy = boosted.GetKineticEn 467 G4double anEnergy = boosted.GetKineticEnergy();
474 theGammaEnergy = anEnergy - theGammaEner << 468 theGammaEnergy = anEnergy-theGammaEnergy;
475 theGammaEnergy += theNuclearMassDifferen 469 theGammaEnergy += theNuclearMassDifference;
476 G4double eBindProducts = 0; 470 G4double eBindProducts = 0;
477 G4double eBindN = 0; 471 G4double eBindN = 0;
478 G4double eBindP = 0; 472 G4double eBindP = 0;
479 G4double eBindD = G4NucleiProperties::Ge << 473 G4double eBindD = G4NucleiProperties::GetBindingEnergy(2,1);
480 G4double eBindT = G4NucleiProperties::Ge << 474 G4double eBindT = G4NucleiProperties::GetBindingEnergy(3,1);
481 G4double eBindHe3 = G4NucleiProperties:: << 475 G4double eBindHe3 = G4NucleiProperties::GetBindingEnergy(3,2);
482 G4double eBindA = G4NucleiProperties::Ge << 476 G4double eBindA = G4NucleiProperties::GetBindingEnergy(4,2);
483 G4int ia = 0; << 477 G4int ia=0;
484 for (auto const & had : *tmpHadrons) { << 478 for(i=0; i<tmpHadrons->size(); i++)
485 if (had->GetDefinition() == G4Neutron: << 479 {
486 eBindProducts += eBindN; << 480 if(tmpHadrons->operator[](i)->GetDefinition() == G4Neutron::Neutron())
487 } << 481 {
488 else if (had->GetDefinition() == G4Pro << 482 eBindProducts+=eBindN;
489 eBindProducts += eBindP; << 483 }
490 } << 484 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Proton::Proton())
491 else if (had->GetDefinition() == G4Deu << 485 {
492 eBindProducts += eBindD; << 486 eBindProducts+=eBindP;
493 } << 487 }
494 else if (had->GetDefinition() == G4Tri << 488 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Deuteron::Deuteron())
495 eBindProducts += eBindT; << 489 {
496 } << 490 eBindProducts+=eBindD;
497 else if (had->GetDefinition() == G4He3 << 491 }
498 eBindProducts += eBindHe3; << 492 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Triton::Triton())
499 } << 493 {
500 else if (had->GetDefinition() == G4Alp << 494 eBindProducts+=eBindT;
501 eBindProducts += eBindA; << 495 }
502 ia++; << 496 else if(tmpHadrons->operator[](i)->GetDefinition() == G4He3::He3())
503 } << 497 {
504 } << 498 eBindProducts+=eBindHe3;
>> 499 }
>> 500 else if(tmpHadrons->operator[](i)->GetDefinition() == G4Alpha::Alpha())
>> 501 {
>> 502 eBindProducts+=eBindA;
>> 503 ia++;
>> 504 }
>> 505 }
>> 506
505 theGammaEnergy += eBindProducts; 507 theGammaEnergy += eBindProducts;
506 << 508
507 #ifdef G4PHPDEBUG 509 #ifdef G4PHPDEBUG
508 if (fManager->GetDEBUG()) << 510 if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply gamma Energy " << theGammaEnergy << " eBindProducts " << eBindProducts << G4endl;
509 G4cout << " G4ParticleHPInelasticBaseF <<
510 << " eBindProducts " << eBindPr <<
511 #endif 511 #endif
512 << 512
513 // Special treatment for Be9 + n -> 2n + << 513 //101111
514 if (theBaseZ == 4 && theBaseA == 9) { << 514 //Special treatment for Be9 + n -> 2n + Be8 -> 2n + a + a
515 // This only valid for G4NDL3.13,,, << 515 if ( (G4int)(theBaseZ+eps) == 4 && (G4int)(theBaseA+eps) == 9 )
516 if (ia == 2 && std::abs(G4NucleiProper << 516 {
517 G4NucleiProperties::GetBindingEnergy(9 << 517 // This only valid for G4NDL3.13,,,
518 theNuclearMassDifference) < CLHEP::keV << 518 if ( std::abs( theNuclearMassDifference -
519 { << 519 ( G4NucleiProperties::GetBindingEnergy( 8 , 4 ) -
520 theGammaEnergy -= (2 * eBindA); << 520 G4NucleiProperties::GetBindingEnergy( 9 , 4 ) ) ) < 1*CLHEP::keV
521 } << 521 && ia == 2 )
522 } << 522 {
523 << 523 theGammaEnergy -= (2*eBindA);
524 if (theGammaEnergy > 0.0) { << 524 }
525 for (G4int iLevel = theGammas.GetNumberOfLev << 525 }
526 G4double e = theGammas.GetLevelEnerg << 526
527 if (e < theGammaEnergy) { << 527 G4ReactionProductVector * theOtherPhotons = 0;
528 thePhotons = theGammas.GetDecayGam << 528 G4int iLevel;
529 theGammaEnergy -= e; << 529 while(theGammaEnergy>=theGammas.GetLevelEnergy(0)) // Loop checking, 11.05.2015, T. Koi
530 break; << 530 {
531 } << 531 for(iLevel=theGammas.GetNumberOfLevels()-1; iLevel>=0; iLevel--)
532 } << 532 {
533 } << 533 if(theGammas.GetLevelEnergy(iLevel)<theGammaEnergy) break;
>> 534 }
>> 535 if(iLevel==0||iLevel==theGammas.GetNumberOfLevels()-1)
>> 536 {
>> 537 theOtherPhotons = theGammas.GetDecayGammas(iLevel);
>> 538 #ifdef G4PHPDEBUG
>> 539 if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply adding gamma from level " << iLevel << theOtherPhotons->operator[](ii)->GetKineticEnergy() << G4endl;
>> 540 #endif
>> 541 }
>> 542 else
>> 543 {
>> 544 G4double random = G4UniformRand();
>> 545 G4double eLow = theGammas.GetLevelEnergy(iLevel);
>> 546 G4double eHigh = theGammas.GetLevelEnergy(iLevel+1);
>> 547 if(random > (eHigh-eLow)/(theGammaEnergy-eLow)) iLevel++;
>> 548 theOtherPhotons = theGammas.GetDecayGammas(iLevel);
>> 549 }
>> 550 if(thePhotons==0) thePhotons = new G4ReactionProductVector;
>> 551 if(theOtherPhotons != 0)
>> 552 {
>> 553 for(unsigned int iii=0; iii<theOtherPhotons->size(); iii++)
>> 554 {
>> 555 thePhotons->push_back(theOtherPhotons->operator[](iii));
>> 556 #ifdef G4PHPDEBUG
>> 557 if( getenv("G4ParticleHPDebug"))
>> 558 G4cout << iii << " G4ParticleHPInelasticBaseFS::BaseApply adding gamma " << theOtherPhotons->operator[](iii)->GetKineticEnergy() << G4endl;
>> 559 #endif
>> 560 }
>> 561 delete theOtherPhotons;
>> 562 }
>> 563 theGammaEnergy -= theGammas.GetLevelEnergy(iLevel);
>> 564 if(iLevel == -1) break;
>> 565 }
534 } 566 }
535 } 567 }
>> 568
536 // fill the result 569 // fill the result
537 std::size_t nSecondaries = tmpHadrons->size( << 570 unsigned int nSecondaries = tmpHadrons->size();
538 std::size_t nPhotons = 0; << 571 unsigned int nPhotons = 0;
539 if (thePhotons != nullptr) { << 572 if(thePhotons!=0) { nPhotons = thePhotons->size(); }
540 nPhotons = thePhotons->size(); <<
541 } <<
542 nSecondaries += nPhotons; 573 nSecondaries += nPhotons;
543 << 574 G4DynamicParticle * theSec;
544 G4DynamicParticle* theSec; <<
545 #ifdef G4PHPDEBUG 575 #ifdef G4PHPDEBUG
546 if (fManager->GetDEBUG()) << 576 if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N hadrons " << nSecondaries-nPhotons << G4endl;
547 G4cout << " G4ParticleHPInelasticBaseFS::B << 577 #endif
548 << G4endl; << 578
549 #endif << 579 for(i=0; i<nSecondaries-nPhotons; i++)
550 << 580 {
551 for (i = 0; i < nSecondaries - nPhotons; ++i << 581 theSec = new G4DynamicParticle;
552 theSec = new G4DynamicParticle; << 582 theSec->SetDefinition(tmpHadrons->operator[](i)->GetDefinition());
553 theSec->SetDefinition(tmpHadrons->operator << 583 theSec->SetMomentum(tmpHadrons->operator[](i)->GetMomentum());
554 theSec->SetMomentum(tmpHadrons->operator[] << 584 theResult.Get()->AddSecondary(theSec);
555 theResult.Get()->AddSecondary(theSec, secI << 585 #ifdef G4PHPDEBUG
556 #ifdef G4PHPDEBUG << 586 if( getenv("G4ParticleHPDebug")) G4cout << this << " G4ParticleHPInelasticBaseFS::BaseApply add secondary2 " << theSec->GetParticleDefinition()->GetParticleName() << " E= " << theSec->GetKineticEnergy() << " NSECO " << theResult.Get()->GetNumberOfSecondaries() << G4endl;
557 if (fManager->GetDEBUG()) << 587 #endif
558 G4cout << " G4ParticleHPInelasticBaseFS: << 588 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( << 589 delete tmpHadrons->operator[](i);
560 << " E=" << theSec->GetKineticEne << 590 }
561 << theResult.Get()->GetNumberOfSe << 591 #ifdef G4PHPDEBUG
562 #endif << 592 if( getenv("G4ParticleHPDebug")) G4cout << " G4ParticleHPInelasticBaseFS::BaseApply N photons " << nPhotons << G4endl;
563 delete (*tmpHadrons)[i]; << 593 #endif
564 } << 594 if(thePhotons != 0)
565 #ifdef G4PHPDEBUG << 595 {
566 if (fManager->GetDEBUG()) << 596 for(i=0; i<nPhotons; i++)
567 G4cout << " G4ParticleHPInelasticBaseFS::B << 597 {
568 #endif << 598 theSec = new G4DynamicParticle;
569 if (thePhotons != nullptr) { << 599 theSec->SetDefinition(thePhotons->operator[](i)->GetDefinition());
570 for (i = 0; i < nPhotons; ++i) { << 600 theSec->SetMomentum(thePhotons->operator[](i)->GetMomentum());
571 theSec = new G4DynamicParticle; << 601 theResult.Get()->AddSecondary(theSec);
572 theSec->SetDefinition((*thePhotons)[i]-> << 602 #ifdef G4PHPDEBUG
573 theSec->SetMomentum((*thePhotons)[i]->Ge << 603 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 604 #endif
582 delete thePhotons->operator[](i); 605 delete thePhotons->operator[](i);
583 } 606 }
584 } 607 }
585 << 608
586 // some garbage collection << 609 // some garbage collection
587 delete thePhotons; 610 delete thePhotons;
588 delete tmpHadrons; 611 delete tmpHadrons;
589 612
590 G4ParticleDefinition* targ_pd = ionTable->Ge << 613 //080721
591 G4LorentzVector targ_4p_lab( << 614 G4ParticleDefinition* targ_pd = G4IonTable::GetIonTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
592 theTarget.GetMomentum(), << 615 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- << 616 G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
594 G4LorentzVector proj_4p_lab = theTrack.Get4M << 617 G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
595 G4LorentzVector init_4p_lab = proj_4p_lab + << 618 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 619
603 // clean up the primary neutron << 620 // clean up the primary neutron
604 theResult.Get()->SetStatusChange(stopAndKill 621 theResult.Get()->SetStatusChange(stopAndKill);
605 } 622 }
606 623