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 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 35 36 #include "globals.hh" 37 38 /** \file G4INCLDeuteronDensity.cc 39 * \brief Deuteron density in r and p accordin 40 * 41 * \date 6 March 2012 42 * \author Davide Mancusi 43 */ 44 45 #include "G4INCLDeuteronDensity.hh" 46 #include "G4INCLGlobals.hh" 47 // #include <cassert> 48 #include <algorithm> 49 50 namespace G4INCL { 51 52 namespace DeuteronDensity { 53 54 namespace { 55 56 const G4int coeffTableSize = 13; 57 58 /// \brief Coefficients for the deuteron 59 const G4double coeff1[coeffTableSize] = 60 0.88688076e+0, 61 -0.34717093e+0, 62 -0.30502380e+1, 63 0.56207766e+2, 64 -0.74957334e+3, 65 0.53365279e+4, 66 -0.22706863e+5, 67 0.60434469e+5, 68 -0.10292058e+6, 69 0.11223357e+6, 70 -0.75925226e+5, 71 0.29059715e+5, 72 -0.48157368e+4 73 }; 74 75 /// \brief Coefficients for the deuteron 76 const G4double coeff2[coeffTableSize] = 77 0.23135193e-1, 78 -0.85604572e+0, 79 0.56068193e+1, 80 -0.69462922e+2, 81 0.41631118e+3, 82 -0.12546621e+4, 83 0.12387830e+4, 84 0.33739172e+4, 85 -0.13041151e+5, 86 0.19512524e+5, 87 -0.15634324e+5, 88 0.66231089e+4, 89 -0.11698185e+4 90 }; 91 92 /// \brief Normalisation coefficient for 93 const G4double normalisationR = std::sqr 94 95 /// \brief Normalisation coefficient for 96 const G4double normalisationP = normalis 97 98 /// \brief Mysterious coefficient that a 99 const G4double al = 0.23162461; 100 101 } 102 103 G4double densityR(const G4double r) { 104 const G4double sWave = wavefunctionR(0, 105 const G4double dWave = wavefunctionR(2, 106 return r*r*(sWave*sWave + dWave*dWave); 107 } 108 109 G4double derivDensityR(const G4double r) { 110 const G4double sWave = wavefunctionR(0, 111 const G4double dWave = wavefunctionR(2, 112 const G4double sWaveDeriv = derivWavefun 113 const G4double dWaveDeriv = derivWavefun 114 return (sWave*sWaveDeriv + dWave*dWaveDe 115 } 116 117 G4double densityP(const G4double p) { 118 const G4double sWave = wavefunctionP(0, 119 const G4double dWave = wavefunctionP(2, 120 return p*p*(sWave*sWave + dWave*dWave); 121 } 122 123 G4double wavefunctionR(const G4int l, cons 124 // assert(l==0 || l==2); // only s- and d-wave 125 const G4double r = 2. * std::max(theR, 1 126 127 G4double result = 0.; 128 G4double fmr; 129 130 for(G4int i=0; i<coeffTableSize; ++i) { 131 fmr = r * (al+i); 132 if(l==0) { // s-wave 133 result += coeff1[i] * std::exp(-fmr) 134 } else { // d-wave 135 result += coeff2[i] * std::exp(-fmr) 136 } 137 } 138 139 result *= normalisationR/r; 140 return result; 141 } 142 143 G4double derivWavefunctionR(const G4int l, 144 // assert(l==0 || l==2); // only s- and d-wave 145 const G4double r = 2. * std::max(theR, 1 146 147 G4double result = 0.; 148 G4double fmr; 149 150 for(G4int i=0; i<coeffTableSize; ++i) { 151 fmr = r * (al+i); 152 if(l==0) { // s-wave 153 result += coeff1[i] * std::exp(-fmr) 154 } else { // d-wave 155 result += coeff2[i] * std::exp(-fmr) 156 } 157 } 158 159 result *= -normalisationR/(r*r); 160 return result; 161 } 162 163 G4double wavefunctionP(const G4int l, cons 164 // assert(l==0 || l==2); // only s- and d-wave 165 const G4double q = theQ / PhysicalConsta 166 const G4double q2 = q*q; 167 G4double result=0.; 168 G4double fmq, alPlusI; 169 for(G4int i=0; i<coeffTableSize; ++i) { 170 alPlusI = al+i; 171 fmq = q2 + alPlusI*alPlusI; 172 if(l==0) { // s-wave 173 result += coeff1[i] / fmq; 174 } else { // d-wave 175 result += coeff2[i] / fmq; 176 } 177 } 178 179 result *= normalisationP; 180 return result; 181 } 182 183 } 184 185 } 186