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