1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
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 | |
44 | inline unsigned int sqr(unsigned int x) { return x*x; } |
45 | inline Int_t sqr(Int_t x) { return x*x; } |
46 | inline Double_t sqr(Double_t x) { return x*x; } |
47 | inline LDouble_t sqr(LDouble_t x) { return x*x; } |
48 | inline Complex_t sqr(Complex_t x) { return x*x; } |
49 | |
50 | const TThreeVectorReal zeroVector(0,0,0); |
51 | const TThreeVectorReal posXhat(1,0,0); |
52 | const TThreeVectorReal negXhat(-1,0,0); |
53 | const TThreeVectorReal posYhat(0,1,0); |
54 | const TThreeVectorReal negYhat(0,-1,0); |
55 | const TThreeVectorReal posZhat(0,0,1); |
56 | const TThreeVectorReal negZhat(0,0,-1); |
57 | |
58 | PairConversionGeneration::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 | |
74 | PairConversionGeneration::~PairConversionGeneration() |
75 | {} |
76 | |
77 | LDouble_t PairConversionGeneration::FFatomic(LDouble_t qRecoil) |
78 | { |
79 | |
80 | |
81 | |
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 | |
90 | |
91 | |
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 | |
104 | |
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 | |
117 | LDouble_t PairConversionGeneration::DiffXS_pair(const TPhoton &gIn, |
118 | const TLepton &pOut, |
119 | const TLepton &eOut) |
120 | { |
121 | |
122 | |
123 | |
124 | |
125 | |
126 | |
127 | |
128 | |
129 | |
130 | |
131 | |
132 | |
133 | |
134 | TPhoton g0(gIn); |
135 | TLepton p1(pOut); |
136 | TLepton e2(eOut); |
137 | |
138 | |
139 | p1.AllPol(); |
140 | e2.AllPol(); |
141 | |
142 | |
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 | |
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 | |
176 | LDouble_t PairConversionGeneration::DiffXS_triplet(const TPhoton &gIn, |
177 | const TLepton &pOut, |
178 | const TLepton &eOut2, |
179 | const TLepton &eOut3) |
180 | { |
181 | |
182 | |
183 | |
184 | |
185 | |
186 | |
187 | |
188 | |
189 | |
190 | |
191 | TPhoton g0(gIn); |
192 | TLepton e0(zeroVector, mElectron); |
193 | TLepton p1(pOut); |
194 | TLepton e2(eOut2); |
195 | TLepton e3(eOut3); |
196 | |
197 | |
198 | |
199 | if (e2.Mass() == e3.Mass() && e3.Mom().Length() > e2.Mom().Length()) |
200 | return 0; |
201 | |
202 | |
203 | e0.SetPol(zeroVector); |
204 | p1.AllPol(); |
205 | e2.AllPol(); |
206 | e3.AllPol(); |
207 | |
208 | |
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 |