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 ]

Diff markup

Differences between /processes/hadronic/models/particle_hp/src/G4ParticleHPElasticFS.cc (Version 11.3.0) and /processes/hadronic/models/particle_hp/src/G4ParticleHPElasticFS.cc (Version 8.2.p1)


  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