| 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 // neutron_hp -- source file 26 // neutron_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 // 25-08-06 New Final State type (refFlag==3 , << 30 // 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) )
31 // is added by T. KOI 31 // is added by T. KOI
32 // 080904 Add Protection for negative energy r << 32 // 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T. Koi
33 // Koi <<
34 // 33 //
35 // P. Arce, June-2014 Conversion neutron_hp to 34 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
36 // 35 //
37 #include "G4ParticleHPElasticFS.hh" 36 #include "G4ParticleHPElasticFS.hh"
>> 37 #include "G4ParticleHPManager.hh"
38 38
39 #include "G4Alpha.hh" << 39 #include "G4PhysicalConstants.hh"
>> 40 #include "G4SystemOfUnits.hh"
>> 41 #include "G4ReactionProduct.hh"
>> 42 #include "G4Nucleus.hh"
>> 43 #include "G4Proton.hh"
40 #include "G4Deuteron.hh" 44 #include "G4Deuteron.hh"
41 #include "G4HadronicParameters.hh" << 45 #include "G4Triton.hh"
42 #include "G4IonTable.hh" << 46 #include "G4Alpha.hh"
>> 47 #include "G4ThreeVector.hh"
43 #include "G4LorentzVector.hh" 48 #include "G4LorentzVector.hh"
44 #include "G4Nucleus.hh" << 49 #include "G4IonTable.hh"
45 #include "G4ParticleHPDataUsed.hh" 50 #include "G4ParticleHPDataUsed.hh"
46 #include "G4ParticleHPManager.hh" <<
47 #include "G4PhysicalConstants.hh" <<
48 #include "G4PhysicsModelCatalog.hh" <<
49 #include "G4Pow.hh" 51 #include "G4Pow.hh"
50 #include "G4Proton.hh" <<
51 #include "G4ReactionProduct.hh" <<
52 #include "G4SystemOfUnits.hh" <<
53 #include "G4ThreeVector.hh" <<
54 #include "G4Triton.hh" <<
55 <<
56 #include "zlib.h" 52 #include "zlib.h"
>> 53 #include "G4PhysicsModelCatalog.hh"
>> 54
57 55
58 G4ParticleHPElasticFS::G4ParticleHPElasticFS() 56 G4ParticleHPElasticFS::G4ParticleHPElasticFS()
59 { 57 {
60 svtEmax = 0.0; << 58 secID = G4PhysicsModelCatalog::GetModelID( "model_NeutronHPElastic" );
61 dbrcEmax = 0.0; << 59
62 dbrcEmin = 0.0; << 60 hasXsec = false;
63 dbrcAmin = 0.0; << 61 theCoefficients = 0;
64 dbrcUse = false; << 62 theProbArray = 0;
65 xsForDBRC = nullptr; << 63
66 <<
67 secID = G4PhysicsModelCatalog::GetModelID("m <<
68 <<
69 hasXsec = false; <<
70 theCoefficients = nullptr; <<
71 theProbArray = nullptr; <<
72 <<
73 repFlag = 0; 64 repFlag = 0;
74 tE_of_repFlag3 = 0.0; 65 tE_of_repFlag3 = 0.0;
75 targetMass = 0.0; 66 targetMass = 0.0;
76 frameFlag = 0; 67 frameFlag = 0;
77 } 68 }
78 69
79 void G4ParticleHPElasticFS::Init(G4double A, G 70 void G4ParticleHPElasticFS::Init(G4double A, G4double Z, G4int M,
80 const G4Strin << 71 G4String& dirName, G4String&,
81 G4ParticleDef << 72 G4ParticleDefinition* )
82 { 73 {
83 G4String tString = "/FS"; 74 G4String tString = "/FS";
84 G4bool dbool = true; << 75 G4bool dbool;
85 SetA_Z(A, Z, M); << 76 G4ParticleHPDataUsed aFile =
86 const G4ParticleHPDataUsed& aFile = << 77 theNames.GetName(static_cast<G4int>(A), static_cast<G4int>(Z), M, dirName, tString, dbool);
87 theNames.GetName(theBaseA, theBaseZ, M, di << 78 G4String filename = aFile.GetName();
88 const G4String& filename = aFile.GetName(); << 79 SetAZMs( A, Z, M, aFile );
89 SetAZMs(aFile); << 80 //theBaseA = aFile.GetA();
>> 81 //theBaseZ = aFile.GetZ();
90 if (!dbool) { 82 if (!dbool) {
91 hasAnyData = false; 83 hasAnyData = false;
92 hasFSData = false; << 84 hasFSData = false;
93 hasXsec = false; 85 hasXsec = false;
94 return; 86 return;
95 } 87 }
96 88
97 // 130205 For compressed data files << 89 //130205 For compressed data files
98 std::istringstream theData(std::ios::in); 90 std::istringstream theData(std::ios::in);
99 G4ParticleHPManager::GetInstance()->GetDataS << 91 G4ParticleHPManager::GetInstance()->GetDataStream(filename,theData);
100 // 130205 END << 92 //130205 END
101 theData >> repFlag >> targetMass >> frameFla 93 theData >> repFlag >> targetMass >> frameFlag;
102 94
103 if (repFlag == 1) { 95 if (repFlag == 1) {
104 G4int nEnergy; 96 G4int nEnergy;
105 theData >> nEnergy; << 97 theData >> nEnergy;
106 theCoefficients = new G4ParticleHPLegendre 98 theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
107 theCoefficients->InitInterpolation(theData 99 theCoefficients->InitInterpolation(theData);
108 G4double temp, energy; 100 G4double temp, energy;
109 G4int tempdep, nLegendre; 101 G4int tempdep, nLegendre;
110 G4int i, ii; 102 G4int i, ii;
111 for (i = 0; i < nEnergy; i++) { << 103 for (i=0; i < nEnergy; i++) {
112 theData >> temp >> energy >> tempdep >> 104 theData >> temp >> energy >> tempdep >> nLegendre;
113 energy *= eV; << 105 energy *=eV;
114 theCoefficients->Init(i, energy, nLegend 106 theCoefficients->Init(i, energy, nLegendre);
115 theCoefficients->SetTemperature(i, temp) 107 theCoefficients->SetTemperature(i, temp);
116 G4double coeff = 0; 108 G4double coeff = 0;
117 for (ii = 0; ii < nLegendre; ii++) { 109 for (ii = 0; ii < nLegendre; ii++) {
118 // load legendre coefficients. 110 // load legendre coefficients.
119 theData >> coeff; 111 theData >> coeff;
120 theCoefficients->SetCoeff(i, ii + 1, c << 112 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
121 } 113 }
122 } 114 }
123 } << 115
124 else if (repFlag == 2) { << 116 } else if (repFlag == 2) {
125 G4int nEnergy; 117 G4int nEnergy;
126 theData >> nEnergy; 118 theData >> nEnergy;
127 theProbArray = new G4ParticleHPPartial(nEn 119 theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
128 theProbArray->InitInterpolation(theData); 120 theProbArray->InitInterpolation(theData);
129 G4double temp, energy; 121 G4double temp, energy;
130 G4int tempdep, nPoints; 122 G4int tempdep, nPoints;
131 for (G4int i = 0; i < nEnergy; i++) { 123 for (G4int i = 0; i < nEnergy; i++) {
132 theData >> temp >> energy >> tempdep >> 124 theData >> temp >> energy >> tempdep >> nPoints;
133 energy *= eV; 125 energy *= eV;
134 theProbArray->InitInterpolation(i, theDa 126 theProbArray->InitInterpolation(i, theData);
135 theProbArray->SetT(i, temp); 127 theProbArray->SetT(i, temp);
136 theProbArray->SetX(i, energy); 128 theProbArray->SetX(i, energy);
137 G4double prob, costh; 129 G4double prob, costh;
138 for (G4int ii = 0; ii < nPoints; ii++) { 130 for (G4int ii = 0; ii < nPoints; ii++) {
139 // fill probability arrays. 131 // fill probability arrays.
140 theData >> costh >> prob; 132 theData >> costh >> prob;
141 theProbArray->SetX(i, ii, costh); 133 theProbArray->SetX(i, ii, costh);
142 theProbArray->SetY(i, ii, prob); 134 theProbArray->SetY(i, ii, prob);
143 } 135 }
144 theProbArray->DoneSetXY(i); << 136 theProbArray->DoneSetXY( i );
145 } 137 }
146 } << 138
147 else if (repFlag == 3) { << 139 } else if (repFlag == 3) {
148 G4int nEnergy_Legendre; 140 G4int nEnergy_Legendre;
149 theData >> nEnergy_Legendre; << 141 theData >> nEnergy_Legendre;
150 if (nEnergy_Legendre <= 0) { << 142 if (nEnergy_Legendre <= 0 ) {
151 std::stringstream iss; 143 std::stringstream iss;
152 iss << "G4ParticleHPElasticFS::Init Data 144 iss << "G4ParticleHPElasticFS::Init Data Error repFlag is 3 but nEnergy_Legendre <= 0";
153 iss << "Z, A and M of problematic file i << 145 iss << "Z, A and M of problematic file is " << theNDLDataZ << ", "
154 << theNDLDataM << " respectively."; << 146 << theNDLDataA << " and " << theNDLDataM << " respectively.";
155 throw G4HadronicException(__FILE__, __LI << 147 throw G4HadronicException(__FILE__, __LINE__, iss.str() );
156 } 148 }
157 theCoefficients = new G4ParticleHPLegendre << 149 theCoefficients = new G4ParticleHPLegendreStore( nEnergy_Legendre );
158 theCoefficients->InitInterpolation(theData << 150 theCoefficients->InitInterpolation( theData );
159 G4double temp, energy; 151 G4double temp, energy;
160 G4int tempdep, nLegendre; 152 G4int tempdep, nLegendre;
161 153
162 for (G4int i = 0; i < nEnergy_Legendre; i+ 154 for (G4int i = 0; i < nEnergy_Legendre; i++) {
163 theData >> temp >> energy >> tempdep >> 155 theData >> temp >> energy >> tempdep >> nLegendre;
164 energy *= eV; << 156 energy *=eV;
165 theCoefficients->Init(i, energy, nLegend << 157 theCoefficients->Init( i , energy , nLegendre );
166 theCoefficients->SetTemperature(i, temp) << 158 theCoefficients->SetTemperature( i , temp );
167 G4double coeff = 0; 159 G4double coeff = 0;
168 for (G4int ii = 0; ii < nLegendre; ii++) 160 for (G4int ii = 0; ii < nLegendre; ii++) {
169 // load legendre coefficients. 161 // load legendre coefficients.
170 theData >> coeff; 162 theData >> coeff;
171 theCoefficients->SetCoeff(i, ii + 1, c << 163 theCoefficients->SetCoeff(i, ii+1, coeff); // @@@HPW@@@
172 } 164 }
173 } << 165 }
174 166
175 tE_of_repFlag3 = energy; << 167 tE_of_repFlag3 = energy;
176 168
177 G4int nEnergy_Prob; 169 G4int nEnergy_Prob;
178 theData >> nEnergy_Prob; 170 theData >> nEnergy_Prob;
179 theProbArray = new G4ParticleHPPartial(nEn << 171 theProbArray = new G4ParticleHPPartial( nEnergy_Prob , nEnergy_Prob );
180 theProbArray->InitInterpolation(theData); << 172 theProbArray->InitInterpolation( theData );
181 G4int nPoints; 173 G4int nPoints;
182 for (G4int i = 0; i < nEnergy_Prob; i++) { 174 for (G4int i = 0; i < nEnergy_Prob; i++) {
183 theData >> temp >> energy >> tempdep >> 175 theData >> temp >> energy >> tempdep >> nPoints;
184 energy *= eV; 176 energy *= eV;
185 177
186 // consistency check 178 // consistency check
187 if (i == 0) 179 if (i == 0)
188 // if ( energy != tE_of_repFlag3 ) //1 << 180 //if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines
189 if (std::abs(energy - tE_of_repFlag3) 181 if (std::abs(energy - tE_of_repFlag3) / tE_of_repFlag3 > 1.0e-15)
190 G4cout << "Warning Transition Energy << 182 G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
191 183
192 theProbArray->InitInterpolation(i, theDa << 184 theProbArray->InitInterpolation( i , theData );
193 theProbArray->SetT(i, temp); << 185 theProbArray->SetT( i , temp );
194 theProbArray->SetX(i, energy); << 186 theProbArray->SetX( i , energy );
195 G4double prob, costh; 187 G4double prob, costh;
196 for (G4int ii = 0; ii < nPoints; ii++) { 188 for (G4int ii = 0; ii < nPoints; ii++) {
197 // fill probability arrays. 189 // fill probability arrays.
198 theData >> costh >> prob; 190 theData >> costh >> prob;
199 theProbArray->SetX(i, ii, costh); << 191 theProbArray->SetX( i , ii , costh );
200 theProbArray->SetY(i, ii, prob); << 192 theProbArray->SetY( i , ii , prob );
201 } 193 }
202 theProbArray->DoneSetXY(i); << 194 theProbArray->DoneSetXY( i );
203 } 195 }
>> 196
>> 197 } else if (repFlag==0) {
>> 198 theData >> frameFlag;
>> 199
>> 200 } else {
>> 201 G4cout << "unusable number for repFlag: repFlag="<<repFlag<<G4endl;
>> 202 throw G4HadronicException(__FILE__, __LINE__, "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
204 } 203 }
205 else if (repFlag == 0) { << 204 //130205 For compressed data files(theData changed from ifstream to istringstream)
206 theData >> frameFlag; << 205 //theData.close();
207 } <<
208 else { <<
209 G4cout << "unusable number for repFlag: re <<
210 throw G4HadronicException(__FILE__, __LINE <<
211 "G4ParticleHPEla <<
212 } <<
213 // 130205 For compressed data files(theData <<
214 // theData.close(); <<
215 } 206 }
216 207
217 G4HadFinalState* G4ParticleHPElasticFS::ApplyY << 208
218 { << 209 G4HadFinalState*
219 if (theResult.Get() == nullptr) theResult.Pu << 210 G4ParticleHPElasticFS::ApplyYourself(const G4HadProjectile& theTrack)
>> 211 {
>> 212 if (theResult.Get() == NULL) theResult.Put(new G4HadFinalState);
220 theResult.Get()->Clear(); 213 theResult.Get()->Clear();
221 G4double eKinetic = theTrack.GetKineticEnerg 214 G4double eKinetic = theTrack.GetKineticEnergy();
222 const G4HadProjectile* incidentParticle = &t << 215 const G4HadProjectile *incidentParticle = &theTrack;
223 G4ReactionProduct theNeutron( << 216 G4ReactionProduct theNeutron(const_cast<G4ParticleDefinition*>(incidentParticle->GetDefinition() ));
224 const_cast<G4ParticleDefinition*>(incident << 217 theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect() );
225 theNeutron.SetMomentum(incidentParticle->Get <<
226 theNeutron.SetKineticEnergy(eKinetic); 218 theNeutron.SetKineticEnergy(eKinetic);
227 219
>> 220 G4ReactionProduct theTarget;
>> 221 G4Nucleus aNucleus;
228 G4ThreeVector neuVelo = 222 G4ThreeVector neuVelo =
229 (1. / incidentParticle->GetDefinition()->G << 223 (1./incidentParticle->GetDefinition()->GetPDGMass())*theNeutron.GetMomentum();
230 G4ReactionProduct theTarget = << 224 theTarget =
231 GetBiasedThermalNucleus(targetMass, neuVel << 225 aNucleus.GetBiasedThermalNucleus(targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
232 226
233 // Neutron and target defined as G4ReactionP 227 // Neutron and target defined as G4ReactionProducts
234 // Prepare Lorentz transformation to lab 228 // Prepare Lorentz transformation to lab
235 229
236 G4ThreeVector the3Neutron = theNeutron.GetMo 230 G4ThreeVector the3Neutron = theNeutron.GetMomentum();
237 G4double nEnergy = theNeutron.GetTotalEnergy 231 G4double nEnergy = theNeutron.GetTotalEnergy();
238 G4ThreeVector the3Target = theTarget.GetMome 232 G4ThreeVector the3Target = theTarget.GetMomentum();
239 G4double tEnergy = theTarget.GetTotalEnergy( 233 G4double tEnergy = theTarget.GetTotalEnergy();
240 G4ReactionProduct theCMS; 234 G4ReactionProduct theCMS;
241 G4double totE = nEnergy + tEnergy; << 235 G4double totE = nEnergy+tEnergy;
242 G4ThreeVector the3CMS = the3Target + the3Neu << 236 G4ThreeVector the3CMS = the3Target+the3Neutron;
243 theCMS.SetMomentum(the3CMS); 237 theCMS.SetMomentum(the3CMS);
244 G4double cmsMom = std::sqrt(the3CMS * the3CM << 238 G4double cmsMom = std::sqrt(the3CMS*the3CMS);
245 G4double sqrts = std::sqrt((totE - cmsMom) * << 239 G4double sqrts = std::sqrt((totE-cmsMom)*(totE+cmsMom));
246 theCMS.SetMass(sqrts); 240 theCMS.SetMass(sqrts);
247 theCMS.SetTotalEnergy(totE); 241 theCMS.SetTotalEnergy(totE);
248 << 242
249 // Data come as function of n-energy in nucl 243 // Data come as function of n-energy in nuclear rest frame
250 G4ReactionProduct boosted; 244 G4ReactionProduct boosted;
251 boosted.Lorentz(theNeutron, theTarget); 245 boosted.Lorentz(theNeutron, theTarget);
252 eKinetic = boosted.GetKineticEnergy(); // g << 246 eKinetic = boosted.GetKineticEnergy(); // get kinetic energy for scattering
253 G4double cosTh = -2; 247 G4double cosTh = -2;
254 248
255 if (repFlag == 1) { 249 if (repFlag == 1) {
256 cosTh = theCoefficients->SampleElastic(eKi 250 cosTh = theCoefficients->SampleElastic(eKinetic);
257 } << 251
258 else if (repFlag == 2) { << 252 } else if (repFlag == 2) {
259 cosTh = theProbArray->Sample(eKinetic); 253 cosTh = theProbArray->Sample(eKinetic);
260 } << 254
261 else if (repFlag == 3) { << 255 } else if (repFlag == 3) {
262 if (eKinetic <= tE_of_repFlag3) { 256 if (eKinetic <= tE_of_repFlag3) {
263 cosTh = theCoefficients->SampleElastic(e 257 cosTh = theCoefficients->SampleElastic(eKinetic);
264 } << 258 } else {
265 else { <<
266 cosTh = theProbArray->Sample(eKinetic); 259 cosTh = theProbArray->Sample(eKinetic);
267 } 260 }
268 } << 261
269 else if (repFlag == 0) { << 262 } else if (repFlag == 0) {
270 cosTh = 2. * G4UniformRand() - 1.; << 263 cosTh = 2.*G4UniformRand() - 1.;
271 } << 264
272 else { << 265 } else {
273 G4cout << "Unusable number for repFlag: re 266 G4cout << "Unusable number for repFlag: repFlag=" << repFlag << G4endl;
274 throw G4HadronicException(__FILE__, __LINE 267 throw G4HadronicException(__FILE__, __LINE__,
275 "G4ParticleHPEla 268 "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
276 } 269 }
277 270
278 if (cosTh < -1.1) { << 271 if (cosTh < -1.1) { return 0; }
279 return nullptr; <<
280 } <<
281 272
282 G4double phi = twopi * G4UniformRand(); << 273 G4double phi = twopi*G4UniformRand();
283 G4double cosPhi = std::cos(phi); 274 G4double cosPhi = std::cos(phi);
284 G4double sinPhi = std::sin(phi); 275 G4double sinPhi = std::sin(phi);
285 G4double theta = std::acos(cosTh); 276 G4double theta = std::acos(cosTh);
286 G4double sinth = std::sin(theta); 277 G4double sinth = std::sin(theta);
287 278
288 if (frameFlag == 1) { 279 if (frameFlag == 1) {
289 // Projectile scattering values cosTh are << 280 // Projectile scattering values cosTh are in target rest frame
290 // In this frame, do relativistic calculat 281 // In this frame, do relativistic calculation of scattered projectile and
291 // target 4-momenta 282 // target 4-momenta
292 283
293 theNeutron.Lorentz(theNeutron, theTarget); 284 theNeutron.Lorentz(theNeutron, theTarget);
294 G4double mN = theNeutron.GetMass(); 285 G4double mN = theNeutron.GetMass();
295 G4double Pinit = theNeutron.GetTotalMoment << 286 G4double Pinit = theNeutron.GetTotalMomentum(); // Incident momentum
296 G4double Einit = theNeutron.GetTotalEnergy << 287 G4double Einit = theNeutron.GetTotalEnergy(); // Incident energy
297 G4double mT = theTarget.GetMass(); 288 G4double mT = theTarget.GetMass();
298 289
299 G4double ratio = mT / mN; << 290 G4double ratio = mT/mN;
300 G4double sqt = std::sqrt(ratio * ratio - 1 << 291 G4double sqt = std::sqrt(ratio*ratio - 1.0 + cosTh*cosTh);
301 G4double beta = Pinit / (mT + Einit); // << 292 G4double beta = Pinit/(mT + Einit); // CMS beta
302 G4double denom = 1. - beta * beta * cosTh << 293 G4double denom = 1. - beta*beta*cosTh*cosTh;
303 G4double term1 = cosTh * (Einit * ratio + << 294 G4double term1 = cosTh*(Einit*ratio + mN)/(mN*ratio + Einit);
304 G4double pN = beta * mN * (term1 + sqt) / << 295 G4double pN = beta*mN*(term1 + sqt)/denom;
305 296
306 // Get the scattered momentum and rotate i 297 // Get the scattered momentum and rotate it in theta and phi
307 G4ThreeVector pDir = theNeutron.GetMomentu << 298 G4ThreeVector pDir = theNeutron.GetMomentum()/Pinit;
308 G4double px = pN * pDir.x(); << 299 G4double px = pN*pDir.x();
309 G4double py = pN * pDir.y(); << 300 G4double py = pN*pDir.y();
310 G4double pz = pN * pDir.z(); << 301 G4double pz = pN*pDir.z();
311 302
312 G4ThreeVector pcmRot; 303 G4ThreeVector pcmRot;
313 pcmRot.setX(px * cosTh * cosPhi - py * sin << 304 pcmRot.setX(px*cosTh*cosPhi - py*sinPhi + pz*sinth*cosPhi);
314 pcmRot.setY(px * cosTh * sinPhi + py * cos << 305 pcmRot.setY(px*cosTh*sinPhi + py*cosPhi + pz*sinth*sinPhi);
315 pcmRot.setZ(-px * sinth + pz * cosTh); << 306 pcmRot.setZ(-px*sinth + pz*cosTh);
316 theNeutron.SetMomentum(pcmRot); 307 theNeutron.SetMomentum(pcmRot);
317 G4double eN = std::sqrt(pN * pN + mN * mN) << 308 G4double eN = std::sqrt(pN*pN + mN*mN); // Scattered neutron energy
318 theNeutron.SetTotalEnergy(eN); 309 theNeutron.SetTotalEnergy(eN);
319 310
320 // Get the scattered target momentum << 311 // Get the scattered target momentum
321 G4ReactionProduct toLab(-1. * theTarget); << 312 G4ReactionProduct toLab(-1.*theTarget);
322 theTarget.SetMomentum(pDir * Pinit - pcmRo << 313 theTarget.SetMomentum(pDir*Pinit - pcmRot);
323 G4double eT = Einit - eN + mT; 314 G4double eT = Einit - eN + mT;
324 theTarget.SetTotalEnergy(eT); 315 theTarget.SetTotalEnergy(eT);
325 316
326 // Now back to lab frame << 317 // Now back to lab frame
327 theNeutron.Lorentz(theNeutron, toLab); 318 theNeutron.Lorentz(theNeutron, toLab);
328 theTarget.Lorentz(theTarget, toLab); 319 theTarget.Lorentz(theTarget, toLab);
329 320
330 // 111005 Protection for not producing 0 k << 321 //111005 Protection for not producing 0 kinetic energy target
331 if (theNeutron.GetKineticEnergy() <= 0) 322 if (theNeutron.GetKineticEnergy() <= 0)
332 theNeutron.SetTotalEnergy(theNeutron.Get << 323 theNeutron.SetTotalEnergy(theNeutron.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
333 * (1. + G4Pow: << 324 if (theTarget.GetKineticEnergy() <= 0)
334 if (theTarget.GetKineticEnergy() <= 0) << 325 theTarget.SetTotalEnergy(theTarget.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
335 theTarget.SetTotalEnergy(theTarget.GetMa << 326
336 } << 327 } else if (frameFlag == 2) {
337 else if (frameFlag == 2) { <<
338 // Projectile scattering values cosTh take 328 // Projectile scattering values cosTh taken from center of mass tabulation
339 329
340 G4LorentzVector proj(nEnergy, the3Neutron) 330 G4LorentzVector proj(nEnergy, the3Neutron);
341 G4LorentzVector targ(tEnergy, the3Target); 331 G4LorentzVector targ(tEnergy, the3Target);
342 G4ThreeVector boostToCM = proj.findBoostTo 332 G4ThreeVector boostToCM = proj.findBoostToCM(targ);
343 proj.boost(boostToCM); 333 proj.boost(boostToCM);
344 targ.boost(boostToCM); 334 targ.boost(boostToCM);
345 335
346 // Rotate projectile and target momenta by 336 // Rotate projectile and target momenta by CM scattering angle
347 // Note: at this point collision axis is n << 337 // Note: at this point collision axis is not along z axis, due to
348 // momentum given target nucleus by 338 // momentum given target nucleus by thermal process
349 G4double px = proj.px(); 339 G4double px = proj.px();
350 G4double py = proj.py(); 340 G4double py = proj.py();
351 G4double pz = proj.pz(); 341 G4double pz = proj.pz();
352 342
353 G4ThreeVector pcmRot; 343 G4ThreeVector pcmRot;
354 pcmRot.setX(px * cosTh * cosPhi - py * sin << 344 pcmRot.setX(px*cosTh*cosPhi - py*sinPhi + pz*sinth*cosPhi);
355 pcmRot.setY(px * cosTh * sinPhi + py * cos << 345 pcmRot.setY(px*cosTh*sinPhi + py*cosPhi + pz*sinth*sinPhi);
356 pcmRot.setZ(-px * sinth + pz * cosTh); << 346 pcmRot.setZ(-px*sinth + pz*cosTh);
357 proj.setVect(pcmRot); 347 proj.setVect(pcmRot);
358 targ.setVect(-pcmRot); 348 targ.setVect(-pcmRot);
359 349
360 // Back to lab frame 350 // Back to lab frame
361 proj.boost(-boostToCM); 351 proj.boost(-boostToCM);
362 targ.boost(-boostToCM); 352 targ.boost(-boostToCM);
363 353
364 theNeutron.SetMomentum(proj.vect()); << 354 theNeutron.SetMomentum(proj.vect() );
365 theNeutron.SetTotalEnergy(proj.e()); << 355 theNeutron.SetTotalEnergy(proj.e() );
366 356
367 theTarget.SetMomentum(targ.vect()); << 357 theTarget.SetMomentum(targ.vect() );
368 theTarget.SetTotalEnergy(targ.e()); << 358 theTarget.SetTotalEnergy(targ.e() );
369 359
370 // 080904 Add Protection for very low ener << 360 //080904 Add Protection for very low energy (1e-6eV) scattering
371 if (theNeutron.GetKineticEnergy() <= 0) { 361 if (theNeutron.GetKineticEnergy() <= 0) {
372 theNeutron.SetTotalEnergy(theNeutron.Get << 362 theNeutron.SetTotalEnergy(theNeutron.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
373 * (1. + G4Pow: <<
374 } 363 }
375 364
376 // 080904 Add Protection for very low ener << 365 //080904 Add Protection for very low energy (1e-6eV) scattering
377 if (theTarget.GetKineticEnergy() <= 0) { 366 if (theTarget.GetKineticEnergy() <= 0) {
378 theTarget.SetTotalEnergy(theTarget.GetMa << 367 theTarget.SetTotalEnergy(theTarget.GetMass()*(1. + G4Pow::GetInstance()->powA(10, -15.65) ) );
379 } 368 }
380 } << 369
381 else { << 370 } else {
382 G4cout << "Value of frameFlag (1=LAB, 2=CM 371 G4cout << "Value of frameFlag (1=LAB, 2=CMS): " << frameFlag;
383 throw G4HadronicException(__FILE__, __LINE 372 throw G4HadronicException(__FILE__, __LINE__,
384 "G4ParticleHPEla << 373 "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect");
385 } 374 }
386 375
387 // Everything is now in the lab frame 376 // Everything is now in the lab frame
388 // Set energy change and momentum change 377 // Set energy change and momentum change
389 theResult.Get()->SetEnergyChange(theNeutron. 378 theResult.Get()->SetEnergyChange(theNeutron.GetKineticEnergy());
390 theResult.Get()->SetMomentumChange(theNeutro 379 theResult.Get()->SetMomentumChange(theNeutron.GetMomentum().unit());
391 380
392 // Make recoil a G4DynamicParticle 381 // Make recoil a G4DynamicParticle
393 auto theRecoil = new G4DynamicParticle; << 382 G4DynamicParticle* theRecoil = new G4DynamicParticle;
394 theRecoil->SetDefinition(G4IonTable::GetIonT 383 theRecoil->SetDefinition(G4IonTable::GetIonTable()->GetIon(static_cast<G4int>(theBaseZ),
395 << 384 static_cast<G4int>(theBaseA), 0) );
396 theRecoil->SetMomentum(theTarget.GetMomentum 385 theRecoil->SetMomentum(theTarget.GetMomentum());
397 theResult.Get()->AddSecondary(theRecoil, sec 386 theResult.Get()->AddSecondary(theRecoil, secID);
398 387
399 // Postpone the tracking of the primary neut 388 // Postpone the tracking of the primary neutron
400 theResult.Get()->SetStatusChange(suspend); 389 theResult.Get()->SetStatusChange(suspend);
401 return theResult.Get(); 390 return theResult.Get();
402 } 391 }
403 392
404 void G4ParticleHPElasticFS::InitializeScatteri <<
405 { <<
406 // Initialize DBRC variables <<
407 svtEmax = G4HadronicParameters::Instance()-> <<
408 G4ParticleHPManager* manager = G4ParticleHPM <<
409 dbrcUse = manager->GetUseDBRC(); <<
410 dbrcEmax = manager->GetMaxEnergyDBRC(); <<
411 dbrcEmin = manager->GetMinEnergyDBRC(); <<
412 dbrcAmin = manager->GetMinADBRC(); <<
413 } <<
414 <<
415 G4ReactionProduct G4ParticleHPElasticFS::GetBi <<
416 <<
417 <<
418 { <<
419 // This new method implements the DBRC (Dopp <<
420 // on top of the SVT (Sampling of the Veloci <<
421 // The SVT algorithm was written by Loic Thu <<
422 // the method G4Nucleus::GetBiasedThermalNuc <<
423 // implemented the DBRC algorithm on top of <<
424 // While the SVT algorithm is still present <<
425 // the DBRC algorithm on top of the SVT one <<
426 // order to avoid a cycle dependency between <<
427 <<
428 InitializeScatteringKernelParameters(); <<
429 <<
430 // Set threshold for SVT algorithm <<
431 G4double E_threshold = svtEmax; <<
432 if (svtEmax == -1.) { <<
433 // If E_neutron <= 400*kB*T (400 is a comm <<
434 // then apply the Sampling ot the Velocity <<
435 // else consider the target nucleus being <<
436 E_threshold = 400.0 * 8.617333262E-11 * te <<
437 } <<
438 <<
439 // If DBRC is enabled and the nucleus is hea <<
440 if (dbrcUse && aMass >= dbrcAmin) { <<
441 E_threshold = std::max(svtEmax, dbrcEmax); <<
442 } <<
443 <<
444 G4double E_neutron = 0.5 * aVelocity.mag2() <<
445 <<
446 G4bool dbrcIsOn = dbrcUse && E_neutron >= db <<
447 <<
448 G4Nucleus aNucleus; <<
449 if (E_neutron > E_threshold || !dbrcIsOn) { <<
450 // Apply only the SVT algorithm, not the D <<
451 return aNucleus.GetBiasedThermalNucleus(ta <<
452 } <<
453 <<
454 G4ReactionProduct result; <<
455 result.SetMass(aMass * G4Neutron::Neutron()- <<
456 <<
457 // Beta = sqrt(m/2kT) <<
458 G4double beta = <<
459 std::sqrt(result.GetMass() <<
460 / (2. * 8.617333262E-11 * temp)) <<
461 <<
462 // Neutron speed vn <<
463 G4double vN_norm = aVelocity.mag(); <<
464 G4double vN_norm2 = vN_norm * vN_norm; <<
465 G4double y = beta * vN_norm; <<
466 <<
467 // Normalize neutron velocity <<
468 aVelocity = (1. / vN_norm) * aVelocity; <<
469 <<
470 // Variables for sampling of target speed an <<
471 G4double x2; <<
472 G4double randThresholdSVT; <<
473 G4double vT_norm, vT_norm2, mu; <<
474 G4double acceptThresholdSVT; <<
475 G4double vRelativeSpeed; <<
476 G4double cdf0 = 2. / (2. + std::sqrt(CLHEP:: <<
477 <<
478 // DBRC variables <<
479 G4double xsRelative = -99.; <<
480 G4double randThresholdDBRC = 0.; <<
481 // Calculate max cross-section in interval f <<
482 G4double eMin = <<
483 0.5 * G4Neutron::Neutron()->GetPDGMass() * <<
484 G4double eMax = <<
485 0.5 * G4Neutron::Neutron()->GetPDGMass() * <<
486 G4double xsMax = xsForDBRC->GetMaxY(eMin, eM <<
487 <<
488 do { <<
489 do { <<
490 // Sample the target velocity vT in the <<
491 if (G4UniformRand() < cdf0) { <<
492 // Sample in C45 from https://laws.lan <<
493 x2 = -std::log(G4UniformRand() * G4Uni <<
494 } <<
495 else { <<
496 // Sample in C61 from https://laws.lan <<
497 G4double ampl = std::cos(CLHEP::pi / 2 <<
498 x2 = -std::log(G4UniformRand()) - std: <<
499 } <<
500 <<
501 vT_norm = std::sqrt(x2) / beta; <<
502 vT_norm2 = vT_norm * vT_norm; <<
503 <<
504 // Sample cosine between the incident ne <<
505 mu = 2. * G4UniformRand() - 1.; <<
506 <<
507 // Define acceptance threshold for SVT <<
508 vRelativeSpeed = std::sqrt(vN_norm2 + vT <<
509 acceptThresholdSVT = vRelativeSpeed / (v <<
510 randThresholdSVT = G4UniformRand(); <<
511 } while (randThresholdSVT >= acceptThresho <<
512 <<
513 // Apply DBRC rejection <<
514 xsRelative = xsForDBRC->GetXsec(0.5 * G4Ne <<
515 * vRelativ <<
516 randThresholdDBRC = G4UniformRand(); <<
517 <<
518 } while (randThresholdDBRC >= xsRelative / x <<
519 <<
520 aNucleus.DoKinematicsOfThermalNucleus(mu, vT <<
521 <<
522 return result; <<
523 } <<
524 393