| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron trans
29 //
30 // 25-08-06 New Final State type (refFlag==3 ,
31 // is added by T. KOI
32 // 080904 Add Protection for negative energy r
33 // Koi
34 //
35 // P. Arce, June-2014 Conversion neutron_hp to
36 //
37 #include "G4ParticleHPElasticFS.hh"
38
39 #include "G4Alpha.hh"
40 #include "G4Deuteron.hh"
41 #include "G4HadronicParameters.hh"
42 #include "G4IonTable.hh"
43 #include "G4LorentzVector.hh"
44 #include "G4Nucleus.hh"
45 #include "G4ParticleHPDataUsed.hh"
46 #include "G4ParticleHPManager.hh"
47 #include "G4PhysicalConstants.hh"
48 #include "G4PhysicsModelCatalog.hh"
49 #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"
57
58 G4ParticleHPElasticFS::G4ParticleHPElasticFS()
59 {
60 svtEmax = 0.0;
61 dbrcEmax = 0.0;
62 dbrcEmin = 0.0;
63 dbrcAmin = 0.0;
64 dbrcUse = false;
65 xsForDBRC = nullptr;
66
67 secID = G4PhysicsModelCatalog::GetModelID("m
68
69 hasXsec = false;
70 theCoefficients = nullptr;
71 theProbArray = nullptr;
72
73 repFlag = 0;
74 tE_of_repFlag3 = 0.0;
75 targetMass = 0.0;
76 frameFlag = 0;
77 }
78
79 void G4ParticleHPElasticFS::Init(G4double A, G
80 const G4Strin
81 G4ParticleDef
82 {
83 G4String tString = "/FS";
84 G4bool dbool = true;
85 SetA_Z(A, Z, M);
86 const G4ParticleHPDataUsed& aFile =
87 theNames.GetName(theBaseA, theBaseZ, M, di
88 const G4String& filename = aFile.GetName();
89 SetAZMs(aFile);
90 if (!dbool) {
91 hasAnyData = false;
92 hasFSData = false;
93 hasXsec = false;
94 return;
95 }
96
97 // 130205 For compressed data files
98 std::istringstream theData(std::ios::in);
99 G4ParticleHPManager::GetInstance()->GetDataS
100 // 130205 END
101 theData >> repFlag >> targetMass >> frameFla
102
103 if (repFlag == 1) {
104 G4int nEnergy;
105 theData >> nEnergy;
106 theCoefficients = new G4ParticleHPLegendre
107 theCoefficients->InitInterpolation(theData
108 G4double temp, energy;
109 G4int tempdep, nLegendre;
110 G4int i, ii;
111 for (i = 0; i < nEnergy; i++) {
112 theData >> temp >> energy >> tempdep >>
113 energy *= eV;
114 theCoefficients->Init(i, energy, nLegend
115 theCoefficients->SetTemperature(i, temp)
116 G4double coeff = 0;
117 for (ii = 0; ii < nLegendre; ii++) {
118 // load legendre coefficients.
119 theData >> coeff;
120 theCoefficients->SetCoeff(i, ii + 1, c
121 }
122 }
123 }
124 else if (repFlag == 2) {
125 G4int nEnergy;
126 theData >> nEnergy;
127 theProbArray = new G4ParticleHPPartial(nEn
128 theProbArray->InitInterpolation(theData);
129 G4double temp, energy;
130 G4int tempdep, nPoints;
131 for (G4int i = 0; i < nEnergy; i++) {
132 theData >> temp >> energy >> tempdep >>
133 energy *= eV;
134 theProbArray->InitInterpolation(i, theDa
135 theProbArray->SetT(i, temp);
136 theProbArray->SetX(i, energy);
137 G4double prob, costh;
138 for (G4int ii = 0; ii < nPoints; ii++) {
139 // fill probability arrays.
140 theData >> costh >> prob;
141 theProbArray->SetX(i, ii, costh);
142 theProbArray->SetY(i, ii, prob);
143 }
144 theProbArray->DoneSetXY(i);
145 }
146 }
147 else if (repFlag == 3) {
148 G4int nEnergy_Legendre;
149 theData >> nEnergy_Legendre;
150 if (nEnergy_Legendre <= 0) {
151 std::stringstream iss;
152 iss << "G4ParticleHPElasticFS::Init Data
153 iss << "Z, A and M of problematic file i
154 << theNDLDataM << " respectively.";
155 throw G4HadronicException(__FILE__, __LI
156 }
157 theCoefficients = new G4ParticleHPLegendre
158 theCoefficients->InitInterpolation(theData
159 G4double temp, energy;
160 G4int tempdep, nLegendre;
161
162 for (G4int i = 0; i < nEnergy_Legendre; i+
163 theData >> temp >> energy >> tempdep >>
164 energy *= eV;
165 theCoefficients->Init(i, energy, nLegend
166 theCoefficients->SetTemperature(i, temp)
167 G4double coeff = 0;
168 for (G4int ii = 0; ii < nLegendre; ii++)
169 // load legendre coefficients.
170 theData >> coeff;
171 theCoefficients->SetCoeff(i, ii + 1, c
172 }
173 }
174
175 tE_of_repFlag3 = energy;
176
177 G4int nEnergy_Prob;
178 theData >> nEnergy_Prob;
179 theProbArray = new G4ParticleHPPartial(nEn
180 theProbArray->InitInterpolation(theData);
181 G4int nPoints;
182 for (G4int i = 0; i < nEnergy_Prob; i++) {
183 theData >> temp >> energy >> tempdep >>
184 energy *= eV;
185
186 // consistency check
187 if (i == 0)
188 // if ( energy != tE_of_repFlag3 ) //1
189 if (std::abs(energy - tE_of_repFlag3)
190 G4cout << "Warning Transition Energy
191
192 theProbArray->InitInterpolation(i, theDa
193 theProbArray->SetT(i, temp);
194 theProbArray->SetX(i, energy);
195 G4double prob, costh;
196 for (G4int ii = 0; ii < nPoints; ii++) {
197 // fill probability arrays.
198 theData >> costh >> prob;
199 theProbArray->SetX(i, ii, costh);
200 theProbArray->SetY(i, ii, prob);
201 }
202 theProbArray->DoneSetXY(i);
203 }
204 }
205 else if (repFlag == 0) {
206 theData >> frameFlag;
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 }
216
217 G4HadFinalState* G4ParticleHPElasticFS::ApplyY
218 {
219 if (theResult.Get() == nullptr) theResult.Pu
220 theResult.Get()->Clear();
221 G4double eKinetic = theTrack.GetKineticEnerg
222 const G4HadProjectile* incidentParticle = &t
223 G4ReactionProduct theNeutron(
224 const_cast<G4ParticleDefinition*>(incident
225 theNeutron.SetMomentum(incidentParticle->Get
226 theNeutron.SetKineticEnergy(eKinetic);
227
228 G4ThreeVector neuVelo =
229 (1. / incidentParticle->GetDefinition()->G
230 G4ReactionProduct theTarget =
231 GetBiasedThermalNucleus(targetMass, neuVel
232
233 // Neutron and target defined as G4ReactionP
234 // Prepare Lorentz transformation to lab
235
236 G4ThreeVector the3Neutron = theNeutron.GetMo
237 G4double nEnergy = theNeutron.GetTotalEnergy
238 G4ThreeVector the3Target = theTarget.GetMome
239 G4double tEnergy = theTarget.GetTotalEnergy(
240 G4ReactionProduct theCMS;
241 G4double totE = nEnergy + tEnergy;
242 G4ThreeVector the3CMS = the3Target + the3Neu
243 theCMS.SetMomentum(the3CMS);
244 G4double cmsMom = std::sqrt(the3CMS * the3CM
245 G4double sqrts = std::sqrt((totE - cmsMom) *
246 theCMS.SetMass(sqrts);
247 theCMS.SetTotalEnergy(totE);
248
249 // Data come as function of n-energy in nucl
250 G4ReactionProduct boosted;
251 boosted.Lorentz(theNeutron, theTarget);
252 eKinetic = boosted.GetKineticEnergy(); // g
253 G4double cosTh = -2;
254
255 if (repFlag == 1) {
256 cosTh = theCoefficients->SampleElastic(eKi
257 }
258 else if (repFlag == 2) {
259 cosTh = theProbArray->Sample(eKinetic);
260 }
261 else if (repFlag == 3) {
262 if (eKinetic <= tE_of_repFlag3) {
263 cosTh = theCoefficients->SampleElastic(e
264 }
265 else {
266 cosTh = theProbArray->Sample(eKinetic);
267 }
268 }
269 else if (repFlag == 0) {
270 cosTh = 2. * G4UniformRand() - 1.;
271 }
272 else {
273 G4cout << "Unusable number for repFlag: re
274 throw G4HadronicException(__FILE__, __LINE
275 "G4ParticleHPEla
276 }
277
278 if (cosTh < -1.1) {
279 return nullptr;
280 }
281
282 G4double phi = twopi * G4UniformRand();
283 G4double cosPhi = std::cos(phi);
284 G4double sinPhi = std::sin(phi);
285 G4double theta = std::acos(cosTh);
286 G4double sinth = std::sin(theta);
287
288 if (frameFlag == 1) {
289 // Projectile scattering values cosTh are
290 // In this frame, do relativistic calculat
291 // target 4-momenta
292
293 theNeutron.Lorentz(theNeutron, theTarget);
294 G4double mN = theNeutron.GetMass();
295 G4double Pinit = theNeutron.GetTotalMoment
296 G4double Einit = theNeutron.GetTotalEnergy
297 G4double mT = theTarget.GetMass();
298
299 G4double ratio = mT / mN;
300 G4double sqt = std::sqrt(ratio * ratio - 1
301 G4double beta = Pinit / (mT + Einit); //
302 G4double denom = 1. - beta * beta * cosTh
303 G4double term1 = cosTh * (Einit * ratio +
304 G4double pN = beta * mN * (term1 + sqt) /
305
306 // Get the scattered momentum and rotate i
307 G4ThreeVector pDir = theNeutron.GetMomentu
308 G4double px = pN * pDir.x();
309 G4double py = pN * pDir.y();
310 G4double pz = pN * pDir.z();
311
312 G4ThreeVector pcmRot;
313 pcmRot.setX(px * cosTh * cosPhi - py * sin
314 pcmRot.setY(px * cosTh * sinPhi + py * cos
315 pcmRot.setZ(-px * sinth + pz * cosTh);
316 theNeutron.SetMomentum(pcmRot);
317 G4double eN = std::sqrt(pN * pN + mN * mN)
318 theNeutron.SetTotalEnergy(eN);
319
320 // Get the scattered target momentum
321 G4ReactionProduct toLab(-1. * theTarget);
322 theTarget.SetMomentum(pDir * Pinit - pcmRo
323 G4double eT = Einit - eN + mT;
324 theTarget.SetTotalEnergy(eT);
325
326 // Now back to lab frame
327 theNeutron.Lorentz(theNeutron, toLab);
328 theTarget.Lorentz(theTarget, toLab);
329
330 // 111005 Protection for not producing 0 k
331 if (theNeutron.GetKineticEnergy() <= 0)
332 theNeutron.SetTotalEnergy(theNeutron.Get
333 * (1. + G4Pow:
334 if (theTarget.GetKineticEnergy() <= 0)
335 theTarget.SetTotalEnergy(theTarget.GetMa
336 }
337 else if (frameFlag == 2) {
338 // Projectile scattering values cosTh take
339
340 G4LorentzVector proj(nEnergy, the3Neutron)
341 G4LorentzVector targ(tEnergy, the3Target);
342 G4ThreeVector boostToCM = proj.findBoostTo
343 proj.boost(boostToCM);
344 targ.boost(boostToCM);
345
346 // Rotate projectile and target momenta by
347 // Note: at this point collision axis is n
348 // momentum given target nucleus by
349 G4double px = proj.px();
350 G4double py = proj.py();
351 G4double pz = proj.pz();
352
353 G4ThreeVector pcmRot;
354 pcmRot.setX(px * cosTh * cosPhi - py * sin
355 pcmRot.setY(px * cosTh * sinPhi + py * cos
356 pcmRot.setZ(-px * sinth + pz * cosTh);
357 proj.setVect(pcmRot);
358 targ.setVect(-pcmRot);
359
360 // Back to lab frame
361 proj.boost(-boostToCM);
362 targ.boost(-boostToCM);
363
364 theNeutron.SetMomentum(proj.vect());
365 theNeutron.SetTotalEnergy(proj.e());
366
367 theTarget.SetMomentum(targ.vect());
368 theTarget.SetTotalEnergy(targ.e());
369
370 // 080904 Add Protection for very low ener
371 if (theNeutron.GetKineticEnergy() <= 0) {
372 theNeutron.SetTotalEnergy(theNeutron.Get
373 * (1. + G4Pow:
374 }
375
376 // 080904 Add Protection for very low ener
377 if (theTarget.GetKineticEnergy() <= 0) {
378 theTarget.SetTotalEnergy(theTarget.GetMa
379 }
380 }
381 else {
382 G4cout << "Value of frameFlag (1=LAB, 2=CM
383 throw G4HadronicException(__FILE__, __LINE
384 "G4ParticleHPEla
385 }
386
387 // Everything is now in the lab frame
388 // Set energy change and momentum change
389 theResult.Get()->SetEnergyChange(theNeutron.
390 theResult.Get()->SetMomentumChange(theNeutro
391
392 // Make recoil a G4DynamicParticle
393 auto theRecoil = new G4DynamicParticle;
394 theRecoil->SetDefinition(G4IonTable::GetIonT
395
396 theRecoil->SetMomentum(theTarget.GetMomentum
397 theResult.Get()->AddSecondary(theRecoil, sec
398
399 // Postpone the tracking of the primary neut
400 theResult.Get()->SetStatusChange(suspend);
401 return theResult.Get();
402 }
403
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