Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLDeuteronDensity.cc
Go to the documentation of this file.
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 // 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 Helsinki Institute of Physics, Finland
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 
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 
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 
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 
93  const G4double normalisationR = std::sqrt(32. * Math::pi) * 0.28212;
94 
96  const G4double normalisationP = normalisationR / (std::sqrt(4. * Math::pi) * std::pow(PhysicalConstants::hc,1.5));
97 
99  const G4double al = 0.23162461;
100 
101  }
102 
104  const G4double sWave = wavefunctionR(0, r);
105  const G4double dWave = wavefunctionR(2, r);
106  return r*r*(sWave*sWave + dWave*dWave);
107  }
108 
110  const G4double sWave = wavefunctionR(0, r);
111  const G4double dWave = wavefunctionR(2, r);
112  const G4double sWaveDeriv = derivWavefunctionR(0, r);
113  const G4double dWaveDeriv = derivWavefunctionR(2, r);
114  return (sWave*sWaveDeriv + dWave*dWaveDeriv) / Math::twoPi;
115  }
116 
118  const G4double sWave = wavefunctionP(0, p);
119  const G4double dWave = wavefunctionP(2, p);
120  return p*p*(sWave*sWave + dWave*dWave);
121  }
122 
123  G4double wavefunctionR(const G4int l, const G4double theR) {
124 // assert(l==0 || l==2); // only s- and d-waves in a deuteron
125  const G4double r = 2. * std::max(theR, 1.e-4);
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) * (1.+3./fmr+3./(fmr*fmr));
136  }
137  }
138 
139  result *= normalisationR/r;
140  return result;
141  }
142 
143  G4double derivWavefunctionR(const G4int l, const G4double theR) {
144 // assert(l==0 || l==2); // only s- and d-waves in a deuteron
145  const G4double r = 2. * std::max(theR, 1.e-4);
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) * (fmr + 1.);
154  } else { // d-wave
155  result += coeff2[i] * std::exp(-fmr) * (fmr + 4. + 9./fmr + 9./(fmr*fmr));
156  }
157  }
158 
159  result *= -normalisationR/(r*r);
160  return result;
161  }
162 
163  G4double wavefunctionP(const G4int l, const G4double theQ) {
164 // assert(l==0 || l==2); // only s- and d-waves in a deuteron
165  const G4double q = theQ / PhysicalConstants::hc;
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 }
G4double G4ParticleHPJENDLHEData::G4double result
G4double derivWavefunctionR(const G4int l, const G4double r)
const char * p
Definition: xmltok.h:285
const G4double pi
const G4double hc
[MeV*fm]
int G4int
Definition: G4Types.hh:78
G4double derivDensityR(const G4double r)
First derivative of the r-space density function.
G4double densityR(const G4double r)
PDF for a nucleon in r space.
Deuteron density in r and p according to the Paris potential.
G4double densityP(const G4double p)
PDF for a nucleon in p space.
G4double wavefunctionP(const G4int l, const G4double p)
G4double wavefunctionR(const G4int l, const G4double r)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4double twoPi
double G4double
Definition: G4Types.hh:76