Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/particle_hp/src/G4ParticleHPElasticFS.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 // neutron_hp -- source file
 27 // J.P. Wellisch, Nov-1996
 28 // A prototype of the low energy neutron transport model.
 29 //
 30 // 25-08-06 New Final State type (refFlag==3 , Legendre (Low Energy) + Probability (High Energy) )
 31 //          is added by T. KOI
 32 // 080904 Add Protection for negative energy results in very low energy ( 1E-6 eV ) scattering by T.
 33 // Koi
 34 //
 35 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
 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("model_NeutronHPElastic");
 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, G4double Z, G4int M,
 80                                  const G4String& dirName, const G4String&,
 81                                  G4ParticleDefinition*)
 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, dirName, tString, dbool);
 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()->GetDataStream(filename, theData);
100   // 130205 END
101   theData >> repFlag >> targetMass >> frameFlag;
102 
103   if (repFlag == 1) {
104     G4int nEnergy;
105     theData >> nEnergy;
106     theCoefficients = new G4ParticleHPLegendreStore(nEnergy);
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 >> nLegendre;
113       energy *= eV;
114       theCoefficients->Init(i, energy, nLegendre);
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, coeff);  // @@@HPW@@@
121       }
122     }
123   }
124   else if (repFlag == 2) {
125     G4int nEnergy;
126     theData >> nEnergy;
127     theProbArray = new G4ParticleHPPartial(nEnergy, nEnergy);
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 >> nPoints;
133       energy *= eV;
134       theProbArray->InitInterpolation(i, theData);
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 Error repFlag is 3 but nEnergy_Legendre <= 0";
153       iss << "Z, A and M of problematic file is " << theNDLDataZ << ", " << theNDLDataA << " and "
154           << theNDLDataM << " respectively.";
155       throw G4HadronicException(__FILE__, __LINE__, iss.str());
156     }
157     theCoefficients = new G4ParticleHPLegendreStore(nEnergy_Legendre);
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 >> nLegendre;
164       energy *= eV;
165       theCoefficients->Init(i, energy, nLegendre);
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, coeff);  // @@@HPW@@@
172       }
173     }
174 
175     tE_of_repFlag3 = energy;
176 
177     G4int nEnergy_Prob;
178     theData >> nEnergy_Prob;
179     theProbArray = new G4ParticleHPPartial(nEnergy_Prob, nEnergy_Prob);
180     theProbArray->InitInterpolation(theData);
181     G4int nPoints;
182     for (G4int i = 0; i < nEnergy_Prob; i++) {
183       theData >> temp >> energy >> tempdep >> nPoints;
184       energy *= eV;
185 
186       // consistency check
187       if (i == 0)
188         // if ( energy != tE_of_repFlag3 ) //110620TK This is too tight for 32bit machines
189         if (std::abs(energy - tE_of_repFlag3) / tE_of_repFlag3 > 1.0e-15)
190           G4cout << "Warning Transition Energy of repFlag3 is not consistent." << G4endl;
191 
192       theProbArray->InitInterpolation(i, theData);
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: repFlag=" << repFlag << G4endl;
210     throw G4HadronicException(__FILE__, __LINE__,
211                               "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
212   }
213   // 130205 For compressed data files(theData changed from ifstream to istringstream)
214   // theData.close();
215 }
216 
217 G4HadFinalState* G4ParticleHPElasticFS::ApplyYourself(const G4HadProjectile& theTrack)
218 {
219   if (theResult.Get() == nullptr) theResult.Put(new G4HadFinalState);
220   theResult.Get()->Clear();
221   G4double eKinetic = theTrack.GetKineticEnergy();
222   const G4HadProjectile* incidentParticle = &theTrack;
223   G4ReactionProduct theNeutron(
224     const_cast<G4ParticleDefinition*>(incidentParticle->GetDefinition()));
225   theNeutron.SetMomentum(incidentParticle->Get4Momentum().vect());
226   theNeutron.SetKineticEnergy(eKinetic);
227 
228   G4ThreeVector neuVelo =
229     (1. / incidentParticle->GetDefinition()->GetPDGMass()) * theNeutron.GetMomentum();
230   G4ReactionProduct theTarget =
231     GetBiasedThermalNucleus(targetMass, neuVelo, theTrack.GetMaterial()->GetTemperature());
232 
233   // Neutron and target defined as G4ReactionProducts
234   // Prepare Lorentz transformation to lab
235 
236   G4ThreeVector the3Neutron = theNeutron.GetMomentum();
237   G4double nEnergy = theNeutron.GetTotalEnergy();
238   G4ThreeVector the3Target = theTarget.GetMomentum();
239   G4double tEnergy = theTarget.GetTotalEnergy();
240   G4ReactionProduct theCMS;
241   G4double totE = nEnergy + tEnergy;
242   G4ThreeVector the3CMS = the3Target + the3Neutron;
243   theCMS.SetMomentum(the3CMS);
244   G4double cmsMom = std::sqrt(the3CMS * the3CMS);
245   G4double sqrts = std::sqrt((totE - cmsMom) * (totE + cmsMom));
246   theCMS.SetMass(sqrts);
247   theCMS.SetTotalEnergy(totE);
248 
249   // Data come as function of n-energy in nuclear rest frame
250   G4ReactionProduct boosted;
251   boosted.Lorentz(theNeutron, theTarget);
252   eKinetic = boosted.GetKineticEnergy();  // get kinetic energy for scattering
253   G4double cosTh = -2;
254 
255   if (repFlag == 1) {
256     cosTh = theCoefficients->SampleElastic(eKinetic);
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(eKinetic);
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: repFlag=" << repFlag << G4endl;
274     throw G4HadronicException(__FILE__, __LINE__,
275                               "G4ParticleHPElasticFS::Init -- unusable number for repFlag");
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 in target rest frame
290     // In this frame, do relativistic calculation of scattered projectile and
291     // target 4-momenta
292 
293     theNeutron.Lorentz(theNeutron, theTarget);
294     G4double mN = theNeutron.GetMass();
295     G4double Pinit = theNeutron.GetTotalMomentum();  // Incident momentum
296     G4double Einit = theNeutron.GetTotalEnergy();  // Incident energy
297     G4double mT = theTarget.GetMass();
298 
299     G4double ratio = mT / mN;
300     G4double sqt = std::sqrt(ratio * ratio - 1.0 + cosTh * cosTh);
301     G4double beta = Pinit / (mT + Einit);  // CMS beta
302     G4double denom = 1. - beta * beta * cosTh * cosTh;
303     G4double term1 = cosTh * (Einit * ratio + mN) / (mN * ratio + Einit);
304     G4double pN = beta * mN * (term1 + sqt) / denom;
305 
306     // Get the scattered momentum and rotate it in theta and phi
307     G4ThreeVector pDir = theNeutron.GetMomentum() / Pinit;
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 * sinPhi + pz * sinth * cosPhi);
314     pcmRot.setY(px * cosTh * sinPhi + py * cosPhi + pz * sinth * sinPhi);
315     pcmRot.setZ(-px * sinth + pz * cosTh);
316     theNeutron.SetMomentum(pcmRot);
317     G4double eN = std::sqrt(pN * pN + mN * mN);  // Scattered neutron energy
318     theNeutron.SetTotalEnergy(eN);
319 
320     // Get the scattered target momentum
321     G4ReactionProduct toLab(-1. * theTarget);
322     theTarget.SetMomentum(pDir * Pinit - pcmRot);
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 kinetic energy target
331     if (theNeutron.GetKineticEnergy() <= 0)
332       theNeutron.SetTotalEnergy(theNeutron.GetMass()
333                                 * (1. + G4Pow::GetInstance()->powA(10, -15.65)));
334     if (theTarget.GetKineticEnergy() <= 0)
335       theTarget.SetTotalEnergy(theTarget.GetMass() * (1. + G4Pow::GetInstance()->powA(10, -15.65)));
336   }
337   else if (frameFlag == 2) {
338     // Projectile scattering values cosTh taken from center of mass tabulation
339 
340     G4LorentzVector proj(nEnergy, the3Neutron);
341     G4LorentzVector targ(tEnergy, the3Target);
342     G4ThreeVector boostToCM = proj.findBoostToCM(targ);
343     proj.boost(boostToCM);
344     targ.boost(boostToCM);
345 
346     // Rotate projectile and target momenta by CM scattering angle
347     // Note: at this point collision axis is not along z axis, due to
348     //       momentum given target nucleus by thermal process
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 * sinPhi + pz * sinth * cosPhi);
355     pcmRot.setY(px * cosTh * sinPhi + py * cosPhi + pz * sinth * sinPhi);
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 energy (1e-6eV) scattering
371     if (theNeutron.GetKineticEnergy() <= 0) {
372       theNeutron.SetTotalEnergy(theNeutron.GetMass()
373                                 * (1. + G4Pow::GetInstance()->powA(10, -15.65)));
374     }
375 
376     // 080904 Add Protection for very low energy (1e-6eV) scattering
377     if (theTarget.GetKineticEnergy() <= 0) {
378       theTarget.SetTotalEnergy(theTarget.GetMass() * (1. + G4Pow::GetInstance()->powA(10, -15.65)));
379     }
380   }
381   else {
382     G4cout << "Value of frameFlag (1=LAB, 2=CMS): " << frameFlag;
383     throw G4HadronicException(__FILE__, __LINE__,
384                               "G4ParticleHPElasticFS::ApplyYourSelf frameflag incorrect");
385   }
386 
387   // Everything is now in the lab frame
388   // Set energy change and momentum change
389   theResult.Get()->SetEnergyChange(theNeutron.GetKineticEnergy());
390   theResult.Get()->SetMomentumChange(theNeutron.GetMomentum().unit());
391 
392   // Make recoil a G4DynamicParticle
393   auto theRecoil = new G4DynamicParticle;
394   theRecoil->SetDefinition(G4IonTable::GetIonTable()->GetIon(static_cast<G4int>(theBaseZ),
395                                                              static_cast<G4int>(theBaseA), 0));
396   theRecoil->SetMomentum(theTarget.GetMomentum());
397   theResult.Get()->AddSecondary(theRecoil, secID);
398 
399   // Postpone the tracking of the primary neutron
400   theResult.Get()->SetStatusChange(suspend);
401   return theResult.Get();
402 }
403 
404 void G4ParticleHPElasticFS::InitializeScatteringKernelParameters()
405 {
406   // Initialize DBRC variables
407   svtEmax = G4HadronicParameters::Instance()->GetNeutronKineticEnergyThresholdForSVT();
408   G4ParticleHPManager* manager = G4ParticleHPManager::GetInstance();
409   dbrcUse = manager->GetUseDBRC();
410   dbrcEmax = manager->GetMaxEnergyDBRC();
411   dbrcEmin = manager->GetMinEnergyDBRC();
412   dbrcAmin = manager->GetMinADBRC();
413 }
414 
415 G4ReactionProduct G4ParticleHPElasticFS::GetBiasedThermalNucleus(const G4double aMass,
416                                                                  G4ThreeVector aVelocity,
417                                                                  const G4double temp)
418 {
419   // This new method implements the DBRC (Doppler Broadening Rejection Correction) algorithm
420   // on top of the SVT (Sampling of the Velocity of the Target nucleus) algorithm.
421   // The SVT algorithm was written by Loic Thulliez (CEA-Saclay) on 2021/05/04 in
422   // the method G4Nucleus::GetBiasedThermalNucleus; Marek Zmeskal on 2022/11/30
423   // implemented the DBRC algorithm on top of the SVT one.
424   // While the SVT algorithm is still present also in G4Nucleus::GetBiasedThermalNucleus,
425   // the DBRC algorithm on top of the SVT one has been moved in this new method, in
426   // order to avoid a cycle dependency between hadronic/util and hadronic/model/particle_hp.
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 common value encounter in MC neutron transport code)
434     // then apply the Sampling ot the Velocity of the Target (SVT) method;
435     // else consider the target nucleus being without motion
436     E_threshold = 400.0 * 8.617333262E-11 * temp;
437   }
438 
439   // If DBRC is enabled and the nucleus is heavy enough, then update the energy threshold
440   if (dbrcUse && aMass >= dbrcAmin) {
441     E_threshold = std::max(svtEmax, dbrcEmax);
442   }
443 
444   G4double E_neutron = 0.5 * aVelocity.mag2() * G4Neutron::Neutron()->GetPDGMass();  // E=0.5*m*v2
445 
446   G4bool dbrcIsOn = dbrcUse && E_neutron >= dbrcEmin && aMass >= dbrcAmin && E_neutron <= dbrcEmax;
447 
448   G4Nucleus aNucleus;
449   if (E_neutron > E_threshold || !dbrcIsOn) {
450     // Apply only the SVT algorithm, not the DBRC one
451     return aNucleus.GetBiasedThermalNucleus(targetMass, aVelocity, temp);
452   }
453 
454   G4ReactionProduct result;
455   result.SetMass(aMass * G4Neutron::Neutron()->GetPDGMass());
456 
457   // Beta = sqrt(m/2kT)
458   G4double beta =
459     std::sqrt(result.GetMass()
460               / (2. * 8.617333262E-11 * temp));  // kT E-5[eV] mass E-11[MeV] => beta in [m/s]-1
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 and SVT rejection
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::pi) * y);
477 
478   // DBRC variables
479   G4double xsRelative = -99.;
480   G4double randThresholdDBRC = 0.;
481   // Calculate max cross-section in interval from  v - 4/beta  to  v + 4/beta  for rejection
482   G4double eMin =
483     0.5 * G4Neutron::Neutron()->GetPDGMass() * (vN_norm - 4. / beta) * (vN_norm - 4. / beta);
484   G4double eMax =
485     0.5 * G4Neutron::Neutron()->GetPDGMass() * (vN_norm + 4. / beta) * (vN_norm + 4. / beta);
486   G4double xsMax = xsForDBRC->GetMaxY(eMin, eMax);
487 
488   do {
489     do {
490       // Sample the target velocity vT in the laboratory frame
491       if (G4UniformRand() < cdf0) {
492         // Sample in C45 from https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf
493         x2 = -std::log(G4UniformRand() * G4UniformRand());
494       }
495       else {
496         // Sample in C61 from https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721.pdf
497         G4double ampl = std::cos(CLHEP::pi / 2.0 * G4UniformRand());
498         x2 = -std::log(G4UniformRand()) - std::log(G4UniformRand()) * ampl * ampl;
499       }
500 
501       vT_norm = std::sqrt(x2) / beta;
502       vT_norm2 = vT_norm * vT_norm;
503 
504       // Sample cosine between the incident neutron and the target in the laboratory frame
505       mu = 2. * G4UniformRand() - 1.;
506 
507       // Define acceptance threshold for SVT
508       vRelativeSpeed = std::sqrt(vN_norm2 + vT_norm2 - 2 * vN_norm * vT_norm * mu);
509       acceptThresholdSVT = vRelativeSpeed / (vN_norm + vT_norm);
510       randThresholdSVT = G4UniformRand();
511     } while (randThresholdSVT >= acceptThresholdSVT);
512 
513     // Apply DBRC rejection
514     xsRelative = xsForDBRC->GetXsec(0.5 * G4Neutron::Neutron()->GetPDGMass() * vRelativeSpeed
515                                     * vRelativeSpeed);
516     randThresholdDBRC = G4UniformRand();
517 
518   } while (randThresholdDBRC >= xsRelative / xsMax);
519 
520   aNucleus.DoKinematicsOfThermalNucleus(mu, vT_norm, aVelocity, result);
521 
522   return result;
523 }
524