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);
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 {
103 std::cerr << "PairConversionGeneration::FFatomic error: "
104 << "no model currently implemented for element "
105 << "Z=" << fConverterZ << std::endl;
106 exit(13);
107 }
108 return ff;
109}
110
111LDouble_t PairConversionGeneration::DiffXS_pair(const TPhoton &gIn,
112 const TLepton &pOut,
113 const TLepton &eOut)
114{
115 // Calculates the e+e- pair production cross section for a
116 // gamma ray off an atom at a particular recoil momentum vector q.
117 // The cross section is returned as d(sigma)/(dE dphi d^3q) where E is
118 // the energy of the final-state electron and phi is its azimuthal angle.
119 // The polar angles of the pair are fixed by momentum conservation.
120 // It is assumed that gIn.Mom()[0] = eOut.Mom()[0]+pOut.Mom()[0], that
121 // the energy carried away by the recoil is zero in the laboratory frame,
122 // but it is not checked. The calculation is performed in the lab frame.
123 // This cross section is only a partial result, because it does not
124 // include the integral d^3 q over the form factor of the target. This
125 // depends on the crystal structure of the target atom, and so is left to
126 // be carried out by more specialized code. Units are barns/GeV^4.
127
128 TPhoton g0(gIn);
129 TLepton p1(pOut);
130 TLepton e2(eOut);
131
132 // Set the initial,final polarizations
133 p1.AllPol();
134 e2.AllPol();
135
136 // Multiply the basic cross section by the converter atomic form factor
137 LDouble_t result = TCrossSection::PairProduction(g0, e2, p1);
138 TFourVectorReal qR(gIn.Mom() - eOut.Mom() - pOut.Mom());
139 result *= sqr(fConverterZ * (1 - FFatomic(qR.Length())));
140 return result * 1e-6;
141
142 // The unpolarized Bethe-Heitler cross section is given here for comparison
143 LDouble_t kin = gIn.Mom()[0];
144 LDouble_t Epos = pOut.Mom()[0];
145 LDouble_t Eele = eOut.Mom()[0];
146 LDouble_t delta = 136 * mElectron / pow(fConverterZ, 0.33333) *
147 kin / (Eele * Epos);
148 LDouble_t aCoul = sqr(alphaQED * fConverterZ);
149 LDouble_t fCoul = aCoul * (1 / (1 + aCoul) + 0.20206 - 0.0369 * aCoul +
150 0.0083 * pow(aCoul, 2) - 0.002 * pow(aCoul, 3));
151 LDouble_t xsi = log(1440 / pow(fConverterZ, 0.66667)) /
152 (log(183 / pow(fConverterZ, 0.33333) - fCoul));
153 LDouble_t FofZ = (8./3.) * log(fConverterZ) + ((kin < 0.05)? 0 : 8 * fCoul);
154 LDouble_t Phi1 = 20.867 - 3.242 * delta + 0.625 * sqr(delta);
155 LDouble_t Phi2 = 20.209 - 1.930 * delta - 0.086 * sqr(delta);
156 LDouble_t Phi0 = 21.12 - 4.184 * log(delta + 0.952);
157 if (delta > 1) {
158 Phi1 = Phi2 = Phi0;
159 }
160 result = hbarcSqr / sqr(mElectron) * pow(alphaQED, 3) / kin
161 * fConverterZ * (fConverterZ + xsi)
162 * (
163 (sqr(Eele) + sqr(Epos)) / sqr(kin) * (Phi1 - FofZ/2) +
164 (2./3.) * (Eele * Epos) / sqr(kin) * (Phi2 - FofZ/2)
165 );
166 return result * 1e-6;
167}
168
169LDouble_t PairConversionGeneration::DiffXS_triplet(const TPhoton &gIn,
170 const TLepton &pOut,
171 const TLepton &eOut2,
172 const TLepton &eOut3)
173{
174 // Calculates the e+e- pair production rate on a free electron target,
175 // including incident photon polarization effects, for a given set of
176 // kinematics. The kinematics are specified by the initial photon
177 // energy kin, the mass of the e+e- pair M, the recoil momentum vector
178 // qR, the azimuthal angle of the plane containing the e+e- pair phi+,
179 // and the energy of the pair positron E+. The returned value is the
180 // differential cross section measured in barns/GeV^4 per atom, with
181 // fConverterZ electrons per atom, differential in
182 // (d^3 qR dphi+ dE+) = (M / 2 kin) (dM dqR^2 dphiR dphi+ dE+).
183
184 TPhoton g0(gIn);
185 TLepton e0(zeroVector, mElectron);
186 TLepton p1(pOut);
187 TLepton e2(eOut2);
188 TLepton e3(eOut3);
189
190 // Avoid double-counting due to identical fs electrons by requiring
191 // that e3 (recoil e-) have a lower momentum magnitude than e2.
192 if (e3.Mom().Length() > e2.Mom().Length())
193 return 0;
194
195 // Set the initial,final polarizations
196 e0.SetPol(zeroVector);
197 p1.AllPol();
198 e2.AllPol();
199 e3.AllPol();
200
201 // Correct the basic cross section by the converter atomic form factor
202 LDouble_t result = TCrossSection::TripletProduction(g0, e0, p1, e2, e3);
203 result *= fConverterZ * (1 - sqr(FFatomic(e3.Mom().Length())));
204 return result * 1e-6;
205}
206
207#endif