Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/PairConversionGeneration.cc
Location:line 100, column 7
Description:Value stored to 'ff' is never read

Annotated Source Code

1//
2// PairConversionGeneration class implementation
3//
4// author: richard.t.jones at uconn.edu
5// version: december 17, 2016
6//
7// notes:
8//
9// This class computes differential pair production rates with full
10// polarization degrees of freedom taken into account for both the
11// photon and the leptons. Both incoherent pair production with a
12// recoil electron and nuclear/coherent atomic scattering without a
13// free recoil electron are included. Calculations are performed by
14// computing directly the Feynman amplitudes for all 8 leading-order
15// graphs and summing them together. No approximations are used.
16//
17// dependencies:
18//
19// 1. This class depends on the Dirac++ package by the same author.
20// 2. The Dirac++ package uses extended precision floating point arithmetic
21// called "long double" in c++. On Intel cpu architectures, this is
22// implemented using 80-bit floating-point registers known as "x87".
23// Extended precision is needed to keep rounding errors in the sum over
24// cancelling Feynman amplitudes from producing excessive rounding
25// errors in the computation of the differential cross section.
26//
27// units:
28// Any length is in m; energy,momentum,mass in GeV (c=1); angles in
29// radians; time in seconds; cross section in barns.
30
31#ifdef USING_DIRACXX1
32#define BOOST_PYTHON_WRAPPING1 1
33
34#include <PairConversionGeneration.hh>
35
36#include "Complex.h"
37#include "TCrossSection.h"
38#include "TLorentzBoost.h"
39
40#include "G4ios.hh"
41
42#include <iostream>
43
44inline unsigned int sqr(unsigned int x) { return x*x; }
45inline Int_t sqr(Int_t x) { return x*x; }
46inline Double_t sqr(Double_t x) { return x*x; }
47inline LDouble_t sqr(LDouble_t x) { return x*x; }
48inline Complex_t sqr(Complex_t x) { return x*x; }
49
50const TThreeVectorReal zeroVector(0,0,0);
51const TThreeVectorReal posXhat(1,0,0);
52const TThreeVectorReal negXhat(-1,0,0);
53const TThreeVectorReal posYhat(0,1,0);
54const TThreeVectorReal negYhat(0,-1,0);
55const TThreeVectorReal posZhat(0,0,1);
56const TThreeVectorReal negZhat(0,0,-1);
57
58PairConversionGeneration::PairConversionGeneration(std::vector<double> Z,
59 std::vector<double> A,
60 std::vector<double> w)
61{
62 if (Z.size() == 1) {
63 fConverterZ = Z[0];
64 }
65 else {
66 std::cerr << "PairConversionGeneration constructor error: "
67 << "converter material is a compound or mixture, " << std::endl
68 << "only pure elements are supported by the "
69 << "present algorithm, cannot continue." << std::endl;
70 exit(13);
71 }
72}
73
74PairConversionGeneration::~PairConversionGeneration()
75{}
76
77LDouble_t PairConversionGeneration::FFatomic(LDouble_t qRecoil)
78{
79 // return the atomic form factor of the pair converter
80 // normalized to unity at zero momentum transfer qRecoil (GeV/c).
81 // Lengths are in Angstroms in this function.
82
83 LDouble_t ff(0);
84 if (fConverterZ == 1) {
85 LDouble_t a0Bohr = 0.529177 / 1.97327e-6;
86 ff = 1 / pow(1 + pow(a0Bohr * qRecoil / 2, 2), 2);
87 }
88 else if (fConverterZ == 4) {
89 // parameterization for 4Be given by online database at
90 // http://lampx.tugraz.at/~hadley/ss1/crystaldiffraction
91 // /atomicformfactors/formfactors.php
92 LDouble_t acoeff[] = {1.5919, 1.1278, 0.5391, 0.7029};
93 LDouble_t bcoeff[] = {43.6427, 1.8623, 103.483, 0.5420};
94 LDouble_t ccoeff[] = {0.0385};
95 LDouble_t q_invA = qRecoil / 1.97327e-6;
96 LDouble_t ff = ccoeff[0];
97 for (int i=0; i < 4; ++i) {
98 ff += acoeff[i] * exp(-bcoeff[i] * pow(q_invA / (4 * M_PI3.14159265358979323846), 2));
99 }
100 ff /= fConverterZ;
Value stored to 'ff' is never read
101 }
102 else if (fConverterZ > 1 && fConverterZ < 93) {
103 // parameterization implemented by Bernard et al
104 // in Geant4 class G4BetheHeitler5DModel.
105 const LDouble_t beta = 2.17e5 / pow(fConverterZ, 1/3.);
106 ff = 1 / (1 + sqr(beta * qRecoil));
107 }
108 else {
109 std::cerr << "PairConversionGeneration::FFatomic error: "
110 << "no model currently implemented for element "
111 << "Z=" << fConverterZ << std::endl;
112 exit(13);
113 }
114 return ff;
115}
116
117LDouble_t PairConversionGeneration::DiffXS_pair(const TPhoton &gIn,
118 const TLepton &pOut,
119 const TLepton &eOut)
120{
121 // Calculates the lepton pair production cross section for a
122 // gamma ray off an atom at a particular recoil momentum vector q.
123 // The cross section is returned as d(sigma)/(dE dphi d^3q) where E is
124 // the energy of the final-state lepton and phi is its azimuthal angle.
125 // The polar angles of the pair are fixed by momentum conservation.
126 // It is assumed that gIn.Mom()[0] = eOut.Mom()[0]+pOut.Mom()[0], that
127 // the energy carried away by the recoil is zero in the laboratory frame,
128 // but it is not checked. The calculation is performed in the lab frame.
129 // This cross section is only a partial result, because it does not
130 // include the integral d^3 q over the form factor of the target. This
131 // depends on the crystal structure of the target atom, and so is left to
132 // be carried out by more specialized code. Units are barns/GeV^4.
133
134 TPhoton g0(gIn);
135 TLepton p1(pOut);
136 TLepton e2(eOut);
137
138 // Set the initial,final polarizations
139 p1.AllPol();
140 e2.AllPol();
141
142 // Multiply the basic cross section by the converter atomic form factor
143 LDouble_t result = TCrossSection::PairProduction(g0, e2, p1);
144 TFourVectorReal qR(gIn.Mom() - eOut.Mom() - pOut.Mom());
145 result *= sqr(fConverterZ * (1 - FFatomic(qR.Length())));
146 return result * 1e-6;
147
148 // The unpolarized Bethe-Heitler cross section is given here for comparison
149 LDouble_t kin = gIn.Mom()[0];
150 LDouble_t Epos = pOut.Mom()[0];
151 LDouble_t Eneg = eOut.Mom()[0];
152 LDouble_t mLepton = eOut.Mass();
153 LDouble_t delta = 136 * mLepton / pow(fConverterZ, 0.33333) *
154 kin / (Eneg * Epos);
155 LDouble_t aCoul = sqr(alphaQED * fConverterZ);
156 LDouble_t fCoul = aCoul * (1 / (1 + aCoul) + 0.20206 - 0.0369 * aCoul +
157 0.0083 * pow(aCoul, 2) - 0.002 * pow(aCoul, 3));
158 LDouble_t xsi = log(1440 / pow(fConverterZ, 0.66667)) /
159 (log(183 / pow(fConverterZ, 0.33333) - fCoul));
160 LDouble_t FofZ = (8./3.) * log(fConverterZ) + ((kin < 0.05)? 0 : 8 * fCoul);
161 LDouble_t Phi1 = 20.867 - 3.242 * delta + 0.625 * sqr(delta);
162 LDouble_t Phi2 = 20.209 - 1.930 * delta - 0.086 * sqr(delta);
163 LDouble_t Phi0 = 21.12 - 4.184 * log(delta + 0.952);
164 if (delta > 1) {
165 Phi1 = Phi2 = Phi0;
166 }
167 result = hbarcSqr / sqr(mLepton) * pow(alphaQED, 3) / kin
168 * fConverterZ * (fConverterZ + xsi)
169 * (
170 (sqr(Eneg) + sqr(Epos)) / sqr(kin) * (Phi1 - FofZ/2) +
171 (2./3.) * (Eneg * Epos) / sqr(kin) * (Phi2 - FofZ/2)
172 );
173 return result * 1e-6;
174}
175
176LDouble_t PairConversionGeneration::DiffXS_triplet(const TPhoton &gIn,
177 const TLepton &pOut,
178 const TLepton &eOut2,
179 const TLepton &eOut3)
180{
181 // Calculates the e+e- pair production rate on a free electron target,
182 // including incident photon polarization effects, for a given set of
183 // kinematics. The kinematics are specified by the initial photon
184 // energy kin, the mass of the e+e- pair M, the recoil momentum vector
185 // qR, the azimuthal angle of the plane containing the e+e- pair phi+,
186 // and the energy of the pair positron E+. The returned value is the
187 // differential cross section measured in barns/GeV^4 per atom, with
188 // fConverterZ electrons per atom, differential in
189 // (d^3 qR dphi+ dE+) = (M / 2 kin) (dM dqR^2 dphiR dphi+ dE+).
190
191 TPhoton g0(gIn);
192 TLepton e0(zeroVector, mElectron);
193 TLepton p1(pOut);
194 TLepton e2(eOut2);
195 TLepton e3(eOut3);
196
197 // Avoid double-counting due to identical fs electrons by requiring
198 // that e3 (recoil e-) have a lower momentum magnitude than e2.
199 if (e2.Mass() == e3.Mass() && e3.Mom().Length() > e2.Mom().Length())
200 return 0;
201
202 // Set the initial,final polarizations
203 e0.SetPol(zeroVector);
204 p1.AllPol();
205 e2.AllPol();
206 e3.AllPol();
207
208 // Correct the basic cross section by the converter atomic form factor
209 LDouble_t result = TCrossSection::TripletProduction(g0, e0, p1, e2, e3);
210 result *= fConverterZ * (1 - sqr(FFatomic(e3.Mom().Length())));
211 return result * 1e-6;
212}
213
214#endif