include/Geant4/G4ParamExpTwoBodyAngDst.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // Author:  Michael Kelsey (SLAC)
0027 // Date:    20 February 2013
0028 //
0029 // Description: implementation base for exponential parametrization of
0030 //               two-body angular distributions in Bertini-style cascade
0031 //
0032 // 20130227  Renamed from "Barashenkov" to "ParamExp" to fix misattribution.
0033 // 20130620  Fix Coverity #50296, recursive #include.
0034 // 20130924  M. Kelsey -- Use G4Log, G4Exp for CPU speedup
0035
0036 #ifndef G4ParamExpTwoBodyAngDst_icc
0037 #define G4ParamExpTwoBodyAngDst_icc 1
0038
0039 #include "G4Log.hh"
0040 #include "G4Exp.hh"
0041 #include "Randomize.hh"
0042 #include <cfloat>
0043
0044
0045 template <G4int NKEBINS> G4double
0046 G4ParamExpTwoBodyAngDst<NKEBINS>::
0047 GetCosTheta(const G4double& ekin, const G4double& pcm) const
0048 {
0049   if (verboseLevel>3) {
0050     G4cout << theName << "::GetCosTheta: ekin " << ekin << " pcm " << pcm
0051        << G4endl;
0052   }
0053
0054   // Get parameter values for interaction energy
0055   G4double pA = interpolator.interpolate(ekin, smallScale);
0056   G4double pC = interpolator.interpolate(ekin, largeScale);
0057   G4double pCos = interpolator.interpolate(ekin, cosScale);
0058   G4double pFrac = interpolator.interpolate(ekin, angleCut);
0059
0060   // Bound parameters by their physical ranges
0061   pCos = std::max(-1., std::min(pCos,1.));
0062   pFrac = std::max(0., std::min(pFrac,1.));
0063
0064   if (verboseLevel>3) {
0065     G4cout << " pFrac " << pFrac << " pA " << pA << " pC " << pC
0066        << " pCos " << pCos << G4endl;
0067   }
0068
0069   G4bool smallAngle = (G4UniformRand() < pFrac);    // 0 < theta < theta0
0070
0071   G4double term1 = 2.0 * pcm*pcm * (smallAngle ? pA : pC);
0072
0073   if (std::abs(term1) < 1e-7) return 1.0;   // No actual scattering here!
0074   if (term1 > DBL_MAX_EXP) return 1.0;
0075
0076   G4double term2 = G4Exp(-2.0*term1);
0077   G4double randScale = (G4Exp(-term1*(1.0 - pCos)) - term2)/(1.0 - term2);
0078
0079   G4double randVal;
0080   if (smallAngle) randVal = (1.0 - randScale)*G4UniformRand() + randScale;
0081   else randVal = randScale*G4UniformRand();
0082
0083   G4double costh = 1.0 + G4Log(randVal*(1.0 - term2) + term2)/term1;
0084
0085   if (verboseLevel>3) {
0086     G4cout << " term1 " << term1 << " term2 " << term2 << " randVal "
0087        << randVal << " => costheta " << costh << G4endl;
0088   }
0089
0090   return costh;
0091 }
0092
0093 #endif  /* G4ParamExpTwoBodyAngDst_icc */