File: | scratch/gluex/scan-build-work/hdgeant4^hdr4251/src/GlueXPhotonBeamGenerator.cc |
Location: | line 107, column 11 |
Description: | Value stored to 'raddz' during its initialization is never read |
1 | // |
2 | // class implementation for GlueXPhotonBeamGenerator |
3 | // |
4 | // author: richard.t.jones at uconn.edu |
5 | // version: december 24, 2016 |
6 | |
7 | #include "GlueXPhotonBeamGenerator.hh" |
8 | #include "GlueXPrimaryGeneratorAction.hh" |
9 | #include "GlueXUserEventInformation.hh" |
10 | #include "GlueXUserOptions.hh" |
11 | #include "HddmOutput.hh" |
12 | #include <G4PhysicalConstants.hh> |
13 | #include <G4SystemOfUnits.hh> |
14 | #include <Randomize.hh> |
15 | |
16 | #include <JANA/JApplication.h> |
17 | #include <JANA/JCalibration.h> |
18 | |
19 | // This flag allows HDGeant4 to be used as a Monte Carlo event |
20 | // generator of a coherent bremsstrahlung beam. When running in |
21 | // "generate not simulate" mode, the output hddm file is filled |
22 | // with MC information describing photons incident on the |
23 | // primary collimator, ie. before collimation. This mode is |
24 | // enabled when the flag GENBEAM 'precol' is present together |
25 | // with the BEAM card in the control.in file. See comments at |
26 | // the head of GlueXBeamConversionProcess.cc for other options |
27 | // related to the GENBEAM control.in flag. |
28 | |
29 | int GlueXPhotonBeamGenerator::fGenerateNotSimulate = 0; |
30 | int GlueXPhotonBeamGenerator::fBeamBackgroundTagOnly = 0; |
31 | |
32 | double GlueXPhotonBeamGenerator::fBeamBucketPeriod = 4. * ns; |
33 | double GlueXPhotonBeamGenerator::fBeamStartZ = -24 * m; |
34 | double GlueXPhotonBeamGenerator::fBeamDiameter = 0.5 * cm; |
35 | double GlueXPhotonBeamGenerator::fBeamVelocity = 2.99792e8 * m/s; |
36 | double GlueXPhotonBeamGenerator::fBeamOffset[2] = {0,0}; |
37 | |
38 | int GlueXPhotonBeamGenerator::fForceFixedPolarization = false; |
39 | double GlueXPhotonBeamGenerator::fFixedPolarization = 0; |
40 | double GlueXPhotonBeamGenerator::fFixedPolarization_phi = 0; |
41 | |
42 | // This utility function is useful in the debugger, |
43 | // but do not use it for actual simulation operatons. |
44 | |
45 | int g4random_seeds(CLHEP::HepRandomEngine *eng=0) { |
46 | const long int *seed = G4RandomG4MTHepRandom::getTheSeeds(); |
47 | G4cout(*G4cout_p) << "current seeds are " << seed[0] << ", " << seed[1] << G4endlstd::endl; |
48 | printf("seeds array pointer is %lx\n", (unsigned long)seed); |
49 | CLHEP::HepRandomEngine *engine = CLHEP::HepRandom::getTheEngine(); |
50 | printf("randoms generator engine pointer is %lx\n", (unsigned long)engine); |
51 | if (eng != 0) |
52 | CLHEP::HepRandom::setTheEngine(eng); |
53 | return 0; |
54 | } |
55 | |
56 | GlueXPhotonBeamGenerator::GlueXPhotonBeamGenerator(CobremsGeneration *gen) |
57 | : fCobrems(gen), |
58 | fTagger(0) |
59 | { |
60 | GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance(); |
61 | if (user_opts == 0) { |
62 | G4cerr(*G4cerr_p) << "Error in GlueXPhotonBeamGenerator constructor - " |
63 | << "GlueXUserOptions::GetInstance() returns null, " |
64 | << "cannot continue." << G4endlstd::endl; |
65 | exit(-1); |
66 | } |
67 | |
68 | std::map<int,std::string> infile; |
69 | std::map<int,double> beampars; |
70 | std::map<int,std::string> genbeampars; |
71 | if (!user_opts->Find("INFILE", infile) && |
72 | user_opts->Find("BEAM", beampars) && |
73 | user_opts->Find("GENBEAM", genbeampars)) |
74 | { |
75 | if (genbeampars.find(1) != genbeampars.end() && |
76 | (genbeampars[1] == "precol" || |
77 | genbeampars[1] == "PRECOL" || |
78 | genbeampars[1] == "Precol" || |
79 | genbeampars[1] == "PreCol" )) |
80 | { |
81 | fGenerateNotSimulate = 1; |
82 | } |
83 | else { |
84 | fGenerateNotSimulate = -1; |
85 | } |
86 | GlueXUserEventInformation::setWriteNoHitEvents(1); |
87 | } |
88 | std::map<int, int> bgratepars; |
89 | std::map<int, int> bggatepars; |
90 | std::map<int, int> bgtagonlypars; |
91 | if (user_opts->Find("BEAM", beampars) && |
92 | user_opts->Find("BGRATE", bgratepars) && |
93 | user_opts->Find("BGGATE", bggatepars)) |
94 | { |
95 | if (user_opts->Find("BGTAGONLY", bgtagonlypars)) { |
96 | fBeamBackgroundTagOnly = bgtagonlypars[1]; |
97 | } |
98 | else { |
99 | fBeamBackgroundTagOnly = 0; |
100 | } |
101 | } |
102 | |
103 | // These cutoffs should be set empirically, as low as possible |
104 | // for good efficiency, but not too low so as to avoid excessive |
105 | // warnings about Pcut violations. |
106 | |
107 | double raddz = fCobrems->getTargetThickness() * m; |
Value stored to 'raddz' during its initialization is never read | |
108 | |
109 | prepareImportanceSamplingPDFs(); |
110 | |
111 | // Create interface for interactive commands |
112 | |
113 | fMessenger = new G4GenericMessenger(this, "/PhotonBeam/", |
114 | "Photon beam generator control"); |
115 | fMessenger->DeclareMethod("enableFixedPolarization", |
116 | &GlueXPhotonBeamGenerator::enableFixedPolarization, |
117 | "Tell the photon beam generator to force a fixed polarization\n" |
118 | " on all beam photons created by the simulation.\n" |
119 | " Arguments are:\n" |
120 | " polarization: value in the range [0,1]\n" |
121 | " phi: azimuthal angle (degrees) [0, 180]"); |
122 | fMessenger->DeclareMethod("disableFixedPolarization", |
123 | &GlueXPhotonBeamGenerator::disableFixedPolarization, |
124 | "Tell the photon beam generator not to force a fixed polarization\n" |
125 | " on beam photons created by the simulation."); |
126 | } |
127 | |
128 | GlueXPhotonBeamGenerator::~GlueXPhotonBeamGenerator() |
129 | { |
130 | delete fMessenger; |
131 | } |
132 | |
133 | void GlueXPhotonBeamGenerator::prepareImportanceSamplingPDFs() |
134 | { |
135 | // Construct lookup tables representing the PDFs used for |
136 | // importance-sampling the coherent bremsstrahlung kinematics. |
137 | |
138 | const int Ndim = 500; |
139 | double Emin = fCobrems->getPhotonEnergyMin() * GeV; |
140 | double Emax = fCobrems->getBeamEnergy() * GeV; |
141 | |
142 | // Compute approximate PDF for dNc/dx |
143 | double xmin = Emin / Emax; |
144 | double dx = (1 - xmin) / Ndim; |
145 | double xarr[Ndim], yarr[Ndim]; |
146 | double ymax = 0; |
147 | for (int i=0; i < Ndim; ++i) { |
148 | xarr[i] = xmin + (i + 0.5) * dx; |
149 | yarr[i] = fCobrems->Rate_dNcdx(xarr[i]); |
150 | yarr[i] = (yarr[i] > 0)? yarr[i] : 0; |
151 | ymax = (yarr[i] > ymax)? yarr[i] : ymax; |
152 | } |
153 | fCobrems->applyBeamCrystalConvolution(Ndim, xarr, yarr); |
154 | double pcut = 5; |
155 | double psum = 0; |
156 | for (int i=0; i < Ndim; ++i) { |
157 | double yplusbl = yarr[i] + ymax / pcut; |
158 | psum += yplusbl; |
159 | fCoherentPDFx.randvar.push_back(xarr[i]); |
160 | fCoherentPDFx.density.push_back(yplusbl); |
161 | fCoherentPDFx.integral.push_back(psum); |
162 | } |
163 | for (int i=0; i < Ndim; ++i) { |
164 | fCoherentPDFx.density[i] /= psum * dx; |
165 | fCoherentPDFx.integral[i] /= psum; |
166 | } |
167 | fCoherentPDFx.Pcut = 3 * psum / Ndim; |
168 | |
169 | // Compute approximate PDF for dNi/dlogx |
170 | double logxmin = log(xmin); |
171 | double dlogx = -logxmin / Ndim; |
172 | double dNidlogx; |
173 | double qsum = 0; |
174 | for (int i=0; i < Ndim; ++i) { |
175 | double logx = logxmin + (i + 0.5) * dlogx; |
176 | double x = exp(logx); |
177 | dNidlogx = fCobrems->Rate_dNidxdt2(x, 0) * x; |
178 | dNidlogx = (dNidlogx > 0)? dNidlogx : 0; |
179 | qsum += dNidlogx; |
180 | fIncoherentPDFlogx.randvar.push_back(logx); |
181 | fIncoherentPDFlogx.density.push_back(dNidlogx); |
182 | fIncoherentPDFlogx.integral.push_back(qsum); |
183 | } |
184 | for (int i=0; i < Ndim; ++i) { |
185 | fIncoherentPDFlogx.density[i] /= qsum * dlogx; |
186 | fIncoherentPDFlogx.integral[i] /= qsum; |
187 | } |
188 | fIncoherentPDFlogx.Pcut = qsum / Ndim; |
189 | |
190 | // Compute approximate PDF for dNi/dy |
191 | fIncoherentPDFtheta02 = 1.8; |
192 | double ymin = 1e-3; |
193 | double dy = (1 - ymin) / Ndim; |
194 | double dNidxdy; |
195 | double tsum = 0; |
196 | for (int i=0; i < Ndim; ++i) { |
197 | double y = ymin + (i + 0.5) * dy; |
198 | double theta2 = fIncoherentPDFtheta02 * (1 / y - 1); |
199 | dNidxdy = fCobrems->Rate_dNidxdt2(0.5, theta2) * |
200 | fIncoherentPDFtheta02 / (y*y); |
201 | dNidxdy = (dNidxdy > 0)? dNidxdy : 0; |
202 | tsum += dNidxdy; |
203 | fIncoherentPDFy.randvar.push_back(y); |
204 | fIncoherentPDFy.density.push_back(dNidxdy); |
205 | fIncoherentPDFy.integral.push_back(tsum); |
206 | } |
207 | for (int i=0; i < Ndim; ++i) { |
208 | fIncoherentPDFy.density[i] /= tsum * dy; |
209 | fIncoherentPDFy.integral[i] /= tsum; |
210 | } |
211 | fIncoherentPDFy.Pcut = 1.5 * tsum / Ndim; |
212 | fCoherentPDFx.Pmax = 0; |
213 | fCoherentPDFx.Psum = 0; |
214 | fIncoherentPDFlogx.Pmax = 0; |
215 | fIncoherentPDFlogx.Psum = 0; |
216 | fIncoherentPDFy.Pmax = 0; |
217 | fIncoherentPDFy.Psum = 0; |
218 | } |
219 | |
220 | void GlueXPhotonBeamGenerator::GeneratePrimaryVertex(G4Event* anEvent) |
221 | { |
222 | GenerateBeamPhoton(anEvent, 0); |
223 | } |
224 | |
225 | void GlueXPhotonBeamGenerator::GenerateBeamPhoton(G4Event* anEvent, double t0) |
226 | { |
227 | // Generates a single beam photon according to the coherent bremsstrahlung |
228 | // model defined by class CobremsGeneration. The photon begins its lifetime |
229 | // just upstream of the primary collimator (WARNING: position is hard-wired |
230 | // in the code below) and is tracked by the simulation from there forward. |
231 | // Its time t0 should identify its beam bucket, ie. the time the photon |
232 | // would reach the midplane of the target. To enable beam motion spreading, |
233 | // define the beam box size below. |
234 | |
235 | // The algorithm below generates coherent bremsstrahlung photons using a |
236 | // importance-sampling technique. This algorithm requires that we prepare |
237 | // an approximate probability density function for the generated photons. |
238 | // The function is not in general equal to the true physical PDF, which |
239 | // varies from event to event depending on the direction of the incident |
240 | // electron, and also the exact angle of the crystal planes at the point |
241 | // of scattering which moves because of the mosaic spread of the crystal. |
242 | // The important thing is that the approximate PDF be reasonably close to |
243 | // the average over all beam particles and the crystal mosaic, and that |
244 | // deviations from event to event are sufficiently small that rejection |
245 | // sampling can be used to take them into account with high efficiency. |
246 | // |
247 | // The kinematics of bremsstrahlung are described by three independent |
248 | // variables (x, theta, phi) where x is the photon energy in units of the |
249 | // incident electron energy, and theta,phi are the polar,azimuthal angles |
250 | // of the photon in a lab frame tilted so that the incident electron comes |
251 | // in along the z axis. Polar angle theta is represented by dimensionless |
252 | // variable y = theta0^2 / (theta^2 + theta0^2) where contant theta0 is |
253 | // chosen to optimize the uniformity of the PDF in y. On each event, |
254 | // a new random tuple (x, phi, y) is generated on the interval x:[0,1], |
255 | // phi:[0,2pi], y:[0,1] using a split-and-recombine strategy. One side |
256 | // of the split covers the coherent process and the other side covers the |
257 | // incoherent one. |
258 | // |
259 | // 1) coherent process - the PDF here is continuous in x,phi according |
260 | // the dNc/(dx dphi), and the dependence on y is a sequence of delta |
261 | // functions corresponding to the different planes that contribute to |
262 | // the scattering at the given value of x. Here we take advantage of |
263 | // the fact that the marginal distribution dNc/dx is proportional to |
264 | // dNc/(dx dphi) at phi=pi/4. This allows us to decompose the generation |
265 | // into two stages, first generating x from dNc/dx and then generating |
266 | // phi from dNc/(dx dphi) at fixed x. The x generation step is performed |
267 | // using importance sampling based on the average PDF stored in table |
268 | // fCoherentPDF, followed by rejection sampling based on the value of |
269 | // dNc/(dx dphi) computed for the particular kinematics of each event. |
270 | // The y value is obtained by sampling the weighted list of q2 values |
271 | // that contributed the to q-sum in the calculation of dNc/(dx dphi). |
272 | // |
273 | // 2) incoherent process - the PDF here is continuous in x,phi,y but it |
274 | // is uniform in phi, so it is effectively a 2D distribution. Here we |
275 | // take advantage of the fact that x and y are independent variables |
276 | // to a good approximation, which allows us to generate x using |
277 | // importance sampling from an approximation to dNi/(dx dtheta^2) at |
278 | // theta=0 and y ~ uniform [0,1], then employ rejection sampling based |
279 | // on the exact PDF dNi/(dx dtheta2) to get a true sample. |
280 | // |
281 | // Recombination after the split is very simple. First we integrate over |
282 | // phi in both cases to obtain values dNc/dx and dNi/(dx dy). It turns |
283 | // out that in spite of the fact that the y-dependence is discrete in the |
284 | // coherent case and continuous in the incoherent case, the sum over the |
285 | // probabilities for all values of y in dNc/dx is always normalized to 1 |
286 | // independently for all values of x. Hence we can treat y as a psuedo |
287 | // coordinate y' ~ Unif[0,1] and form a 2D PDF dNc/(dx dy') which is |
288 | // numerically equal to dNc/dx, do the rejection sampling in combination |
289 | // with that applied to dNi/(dx dy) and then replace the fake variable y' |
290 | // with the true y that was sampled as described above. |
291 | |
292 | assert (fCoherentPDFx.integral.size() != 0)((fCoherentPDFx.integral.size() != 0) ? static_cast<void> (0) : __assert_fail ("fCoherentPDFx.integral.size() != 0", "src/GlueXPhotonBeamGenerator.cc" , 292, __PRETTY_FUNCTION__)); |
293 | |
294 | // The GlueXUserEventInformation constructor can set the random number |
295 | // seeds for this event, so this must happen at here at the top. |
296 | GlueXUserEventInformation *event_info; |
297 | if (t0 == 0) { |
298 | event_info = new GlueXUserEventInformation(); |
299 | anEvent->SetUserInformation(event_info); |
300 | } |
301 | else { |
302 | event_info = (GlueXUserEventInformation*)anEvent->GetUserInformation(); |
303 | } |
304 | assert (event_info != 0)((event_info != 0) ? static_cast<void> (0) : __assert_fail ("event_info != 0", "src/GlueXPhotonBeamGenerator.cc", 304, __PRETTY_FUNCTION__ )); |
305 | |
306 | double phiMosaic = 2*M_PI3.14159265358979323846 * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
307 | double rhoMosaic = sqrt(-2 * log(G4UniformRand()G4MTHepRandom::getTheEngine()->flat())); |
308 | rhoMosaic *= fCobrems->getTargetCrystalMosaicSpread() * radian; |
309 | double thxMosaic = rhoMosaic * cos(phiMosaic); |
310 | double thyMosaic = rhoMosaic * sin(phiMosaic); |
311 | |
312 | double xemittance = fCobrems->getBeamEmittance() * m*radian; |
313 | double yemittance = xemittance / 2.5; // nominal, should be checked |
314 | double xspotsize = fCobrems->getCollimatorSpotrms() * m; |
315 | double yspotsize = xspotsize; // nominal, should be checked |
316 | double phiBeam = 2*M_PI3.14159265358979323846 * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
317 | double rhoBeam = sqrt(-2 * log(G4UniformRand()G4MTHepRandom::getTheEngine()->flat())); |
318 | double thxBeam = (xemittance / xspotsize) * rhoBeam * cos(phiBeam); |
319 | double thyBeam = (yemittance / yspotsize) * rhoBeam * sin(phiBeam); |
320 | |
321 | double raddz = fCobrems->getTargetThickness() * m; |
322 | double varMS = fCobrems->Sigma2MS(raddz/m * G4UniformRand()G4MTHepRandom::getTheEngine()->flat()); |
323 | double rhoMS = sqrt(-2 * varMS * log(G4UniformRand()G4MTHepRandom::getTheEngine()->flat())); |
324 | double phiMS = 2*M_PI3.14159265358979323846 * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
325 | double thxMS = rhoMS * cos(phiMS); |
326 | double thyMS = rhoMS * sin(phiMS); |
327 | |
328 | double targetThetax = fCobrems->getTargetThetax() * radian; |
329 | double targetThetay = fCobrems->getTargetThetay() * radian; |
330 | double targetThetaz = fCobrems->getTargetThetaz() * radian; |
331 | double thetax = thxBeam + thxMS - targetThetax - thxMosaic; |
332 | double thetay = thyBeam + thyMS - targetThetay - thyMosaic; |
333 | double thetaz = -targetThetaz; |
334 | fCobrems->setTargetOrientation(thetax/radian, thetay/radian, thetaz/radian); |
335 | |
336 | // Generate with importance sampling |
337 | double x = 0; |
338 | double phi = 0; |
339 | double theta2 = 0; |
340 | double polarization = 0; |
341 | double polarization_phi = 0; |
342 | double coherent_fraction = fCoherentPDFx.Pcut; |
343 | coherent_fraction /= fCoherentPDFx.Pcut + fIncoherentPDFy.Pcut; |
344 | while (true) { |
345 | double splitrand = G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
346 | if (splitrand < coherent_fraction) { // try coherent generation |
347 | ++fCoherentPDFx.Ntested; |
348 | |
349 | double u = G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
350 | int i = fCoherentPDFx.search(u); |
351 | double ui = fCoherentPDFx.integral[i]; |
352 | double u_i = (i > 0)? fCoherentPDFx.integral[i-1] : 0; |
353 | double xi = fCoherentPDFx.randvar[i]; |
354 | double dx = (i > 0)? xi - fCoherentPDFx.randvar[i-1]: |
355 | fCoherentPDFx.randvar[i+1] - xi; |
356 | x = xi + dx * (0.5 - (ui - u) / (ui - u_i)); |
357 | assert (x > 0)((x > 0) ? static_cast<void> (0) : __assert_fail ("x > 0" , "src/GlueXPhotonBeamGenerator.cc", 357, __PRETTY_FUNCTION__ )); |
358 | double dNcdxPDF = (ui - u_i) / dx; |
359 | double dNcdx = fCobrems->Rate_dNcdx(x); |
360 | double Pfactor = dNcdx / dNcdxPDF; |
361 | if (Pfactor > fCoherentPDFx.Pmax) |
362 | fCoherentPDFx.Pmax = Pfactor; |
363 | if (Pfactor > fCoherentPDFx.Pcut) { |
364 | G4cout(*G4cout_p) << "Warning in GenerateBeamPhoton - Pfactor " << Pfactor |
365 | << " exceeds fCoherentPDFx.Pcut = " << fCoherentPDFx.Pcut |
366 | << G4endlstd::endl |
367 | << " present x = " << x << G4endlstd::endl |
368 | << " present maximum Pfactor = " |
369 | << fCoherentPDFx.Pmax << G4endlstd::endl |
370 | << " current generator efficiency = " |
371 | << fCoherentPDFx.Npassed / |
372 | (fCoherentPDFx.Ntested + 1e-99) |
373 | << G4endlstd::endl; |
374 | } |
375 | ++fCoherentPDFx.Psum += Pfactor; |
376 | if (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() * fCoherentPDFx.Pcut > Pfactor) { |
377 | continue; |
378 | } |
379 | ++fCoherentPDFx.Npassed; |
380 | |
381 | double freq; |
382 | double fmax = dNcdx / M_PI3.14159265358979323846; |
383 | while (true) { |
384 | phi = 2*M_PI3.14159265358979323846 * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
385 | freq = fCobrems->Rate_dNcdxdp(x, phi); |
386 | assert (freq < fmax)((freq < fmax) ? static_cast<void> (0) : __assert_fail ("freq < fmax", "src/GlueXPhotonBeamGenerator.cc", 386, __PRETTY_FUNCTION__ )); |
387 | if (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() * fmax < freq) |
388 | break; |
389 | } |
390 | double uq = freq * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
391 | int j = ImportanceSampler::search(uq, fCobrems->fQ2weight); |
392 | theta2 = fCobrems->fQ2theta2[j]; |
393 | polarization = fCobrems->Polarization(x, theta2, phi); |
394 | polarization_phi = M_PI3.14159265358979323846 / 2; |
395 | break; |
396 | } |
397 | else { // try incoherent generation |
398 | ++fIncoherentPDFlogx.Ntested; |
399 | |
400 | double ux = G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
401 | int i = fIncoherentPDFlogx.search(ux); |
402 | double ui = fIncoherentPDFlogx.integral[i]; |
403 | double u_i = (i > 0)? fIncoherentPDFlogx.integral[i-1] : 0; |
404 | double logxi = fIncoherentPDFlogx.randvar[i]; |
405 | double dlogx = (i > 0)? logxi - fIncoherentPDFlogx.randvar[i-1]: |
406 | fIncoherentPDFlogx.randvar[i+1] - logxi; |
407 | double logx = logxi + dlogx * (0.5 - (ui - ux) / (ui - u_i)); |
408 | assert (logx < 0)((logx < 0) ? static_cast<void> (0) : __assert_fail ( "logx < 0", "src/GlueXPhotonBeamGenerator.cc", 408, __PRETTY_FUNCTION__ )); |
409 | x = exp(logx); |
410 | double dNidxdyPDF = (ui - u_i) / x / dlogx; |
411 | double uy = G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
412 | int j = fIncoherentPDFy.search(uy); |
413 | double uj = fIncoherentPDFy.integral[j]; |
414 | double u_j = (j > 0)? fIncoherentPDFy.integral[j-1] : 0; |
415 | double yj = fIncoherentPDFy.randvar[j]; |
416 | double dy = (j > 0)? yj - fIncoherentPDFy.randvar[j-1]: |
417 | fIncoherentPDFy.randvar[j+1] - yj; |
418 | double y = yj + dy * (0.5 - (uj - uy) / (uj - u_j)); |
419 | assert (y > 0)((y > 0) ? static_cast<void> (0) : __assert_fail ("y > 0" , "src/GlueXPhotonBeamGenerator.cc", 419, __PRETTY_FUNCTION__ )); |
420 | dNidxdyPDF *= (uj - u_j) / dy; |
421 | theta2 = fIncoherentPDFtheta02 * (1 / (y + 1e-99) - 1); |
422 | double dNidxdy = fCobrems->Rate_dNidxdt2(x, theta2) * |
423 | fIncoherentPDFtheta02 / (y*y + 1e-99); |
424 | double Pfactor = dNidxdy / dNidxdyPDF; |
425 | if (Pfactor > fIncoherentPDFlogx.Pmax) |
426 | fIncoherentPDFlogx.Pmax = Pfactor; |
427 | if (Pfactor > fIncoherentPDFy.Pcut) { |
428 | G4cout(*G4cout_p) << "Warning in GenerateBeamPhoton - Pfactor " << Pfactor |
429 | << " exceeds fIncoherentPDFy.Pcut = " |
430 | << fIncoherentPDFy.Pcut |
431 | << G4endlstd::endl |
432 | << " present x = " << x << G4endlstd::endl |
433 | << " present y = " << y << G4endlstd::endl |
434 | << " present maximum Pfactor = " |
435 | << fIncoherentPDFlogx.Pmax << G4endlstd::endl |
436 | << " current generator efficiency = " |
437 | << fIncoherentPDFlogx.Npassed / |
438 | (fIncoherentPDFlogx.Ntested + 1e-99) |
439 | << G4endlstd::endl; |
440 | } |
441 | fIncoherentPDFlogx.Psum += Pfactor; |
442 | if (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() * fIncoherentPDFy.Pcut > Pfactor) { |
443 | continue; |
444 | } |
445 | ++fIncoherentPDFlogx.Npassed; |
446 | |
447 | phi = 2*M_PI3.14159265358979323846 * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
448 | polarization = fCobrems->AbremsPolarization(x, theta2, phi); |
449 | polarization_phi = phi - M_PI3.14159265358979323846 / 2; |
450 | break; |
451 | } |
452 | } |
453 | |
454 | #if VERBOSE_COBREMS_SPLITTING |
455 | if (fIncoherentPDFlogx.Npassed / 100 * 100 == fIncoherentPDFlogx.Npassed) { |
456 | G4cout(*G4cout_p) << "coherent rate is " |
457 | << fCoherentPDFx.Psum / (fCoherentPDFx.Ntested + 1e-99) |
458 | << ", efficiency is " |
459 | << fCoherentPDFx.Npassed / (fCoherentPDFx.Ntested + 1e-99) |
460 | << G4endlstd::endl |
461 | << "incoherent rate is " |
462 | << fIncoherentPDFlogx.Psum / (fIncoherentPDFlogx.Ntested + 1e-99) |
463 | << ", efficiency is " |
464 | << fIncoherentPDFlogx.Npassed / (fIncoherentPDFlogx.Ntested + 1e-99) |
465 | << G4endlstd::endl |
466 | << "counts are " |
467 | << fCoherentPDFx.Npassed << " / " << fIncoherentPDFlogx.Npassed |
468 | << " = " |
469 | << fCoherentPDFx.Npassed / (fIncoherentPDFlogx.Npassed + 1e-99) |
470 | << G4endlstd::endl; |
471 | } |
472 | #endif |
473 | |
474 | // Put the radiator back the way you found it |
475 | fCobrems->setTargetOrientation(targetThetax/radian, |
476 | targetThetay/radian, |
477 | targetThetaz/radian); |
478 | |
479 | // Define the particle kinematics and polarization in lab coordinates |
480 | G4ParticleDefinition *part = GlueXPrimaryGeneratorAction::GetParticle("gamma"); |
481 | double Emax = fCobrems->getBeamEnergy() * GeV; |
482 | double Erms = fCobrems->getBeamErms() * GeV; |
483 | double Ebeam = Emax + Erms * G4RandGaussG4MTRandGaussQ::shoot(); |
484 | double theta = sqrt(theta2) * electron_mass_c2 / Emax; |
485 | double alphax = thxBeam + thxMS + theta * cos(phi); |
486 | double alphay = thyBeam + thyMS + theta * sin(phi); |
487 | double pabs = Ebeam * x; |
488 | double px = pabs * alphax; |
489 | double py = pabs * alphay; |
490 | double pz = sqrt(pabs*pabs - px*px - py*py); |
491 | double colphi = 2*M_PI3.14159265358979323846 * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
492 | double vspotrms = fCobrems->getCollimatorSpotrms() * m; |
493 | double colrho = vspotrms * sqrt(-2 * log(G4UniformRand()G4MTHepRandom::getTheEngine()->flat())); |
494 | double colDist = fCobrems->getCollimatorDistance() * m; |
495 | double radx = colrho * cos(colphi) - colDist * thxBeam; |
496 | double rady = colrho * sin(colphi) - colDist * thyBeam; |
497 | double colx = radx + colDist * alphax; |
498 | double coly = rady + colDist * alphay; |
499 | #if defined BEAM_BOX_SIZE |
500 | colx += BEAM_BOX_SIZE * (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5); |
501 | coly += BEAM_BOX_SIZE * (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5); |
502 | #endif |
503 | colx += fBeamOffset[0]; |
504 | coly += fBeamOffset[1]; |
505 | G4ThreeVector vtx(colx, coly, fBeamStartZ); |
506 | if (fForceFixedPolarization) { |
507 | polarization = fFixedPolarization; |
508 | polarization_phi = fFixedPolarization_phi; |
509 | } |
510 | G4ThreeVector pol(polarization * cos(polarization_phi), |
511 | polarization * sin(polarization_phi), |
512 | -(px * polarization * cos(polarization_phi) + |
513 | py * polarization * sin(polarization_phi)) / pz); |
514 | // use upper half-space to define the polarization plane |
515 | if (pol[2] < 0) |
516 | pol = -pol; |
517 | G4ThreeVector mom(px, py, pz); |
518 | |
519 | // If beam photon is primary particle, use it to initialize event info |
520 | double targetCenterZ = GlueXPrimaryGeneratorAction::getTargetCenterZ(); |
521 | int bg = 1; |
522 | double tvtx; |
523 | if (t0 == 0) { |
524 | tvtx = (vtx[2] - targetCenterZ) / fBeamVelocity; |
525 | tvtx -= GenerateTriggerTime(anEvent); |
526 | if (fGenerateNotSimulate == 0) { |
527 | event_info->AddBeamParticle(1, tvtx, vtx, mom, pol); |
528 | } |
529 | bg = 0; |
530 | } |
531 | else { |
532 | tvtx = fBeamBucketPeriod * floor(t0 / fBeamBucketPeriod + 0.5); |
533 | tvtx += (vtx[2] - targetCenterZ) / fBeamVelocity; |
534 | if (fBeamBackgroundTagOnly) { |
535 | double ttag = tvtx + (targetCenterZ - vtx[2]) / fBeamVelocity; |
536 | GenerateTaggerHit(anEvent, pabs, ttag, bg); |
537 | return; |
538 | } |
539 | } |
540 | |
541 | // Generate new primary for the beam photon |
542 | G4PrimaryVertex* vertex = new G4PrimaryVertex(vtx, tvtx); |
543 | G4PrimaryParticle* photon = new G4PrimaryParticle(part, px, py, pz); |
544 | photon->SetPolarization(pol); |
545 | vertex->SetPrimary(photon); |
546 | if (fGenerateNotSimulate < 1) { |
547 | anEvent->AddPrimaryVertex(vertex); |
548 | } |
549 | |
550 | // Include information about the radiating electron, but do not track it |
551 | G4ParticleDefinition *ve = GlueXPrimaryGeneratorAction::GetParticle("e-"); |
552 | double vpx = Ebeam * thxBeam; |
553 | double vpy = Ebeam * thyBeam; |
554 | double vpz = sqrt(Ebeam*Ebeam - vpx*vpx - vpy*vpy); |
555 | G4PrimaryParticle *velectron = new G4PrimaryParticle(ve, vpx, vpy, vpz); |
556 | double colvx = radx + colDist * thxBeam; |
557 | double colvy = rady + colDist * thyBeam; |
558 | colvx += fBeamOffset[0]; |
559 | colvy += fBeamOffset[1]; |
560 | G4ThreeVector vvtx(colvx, colvy, fBeamStartZ); |
561 | G4PrimaryVertex *vvertex = new G4PrimaryVertex(vvtx, tvtx); |
562 | vvertex->SetPrimary(velectron); |
563 | event_info->AddPrimaryVertex(*vvertex); |
564 | delete vvertex; |
565 | |
566 | // If running in event generation only mode, default is not to |
567 | // save the event to the output file. This will be set back to |
568 | // true if the beam particle makes it to the reference plane. |
569 | if (fGenerateNotSimulate == -1) { |
570 | event_info->SetKeepEvent(0); |
571 | } |
572 | |
573 | // If bg beam particle, append to MC record |
574 | if (bg || fGenerateNotSimulate == 1) { |
575 | event_info->AddPrimaryVertex(*vertex); |
576 | } |
577 | |
578 | if (fGenerateNotSimulate == 0) { |
579 | |
580 | // Register a tagger hit for each beam photon |
581 | |
582 | double ttag = tvtx + (targetCenterZ - vtx[2]) / fBeamVelocity; |
583 | GenerateTaggerHit(anEvent, pabs, ttag, bg); |
584 | } |
585 | } |
586 | |
587 | double GlueXPhotonBeamGenerator::GenerateTriggerTime(const G4Event *event) |
588 | { |
589 | // The primary interaction vertex time is referenced to a clock |
590 | // whose t=0 is synchronized to the crossing of a beam bunch |
591 | // through the clock reference plane. This beam bunch may not |
592 | // contain the beam particle whose interaction generated the |
593 | // vertex, but it represents best-guess based on the arrival time |
594 | // of the L1 trigger signal. The spread in the L1 relative to the |
595 | // interacting bunch time is parameterized as a Gaussian. |
596 | |
597 | static int last_run_number = 0; |
598 | int run_number = HddmOutput::getRunNo(); |
599 | if (run_number != last_run_number) { |
600 | fBeamBucketPeriod = getBeamBucketPeriod(run_number); |
601 | double refZ = getRFreferencePlaneZ(run_number); |
602 | GlueXPrimaryGeneratorAction::setRFreferencePlaneZ(refZ); |
603 | last_run_number = run_number; |
604 | } |
605 | double L1sigmat = GlueXPrimaryGeneratorAction::getL1triggerTimeSigma(); |
606 | double t0 = L1sigmat * G4RandGaussG4MTRandGaussQ::shoot(); |
607 | return fBeamBucketPeriod * floor(t0 / fBeamBucketPeriod + 0.5); |
608 | } |
609 | |
610 | int GlueXPhotonBeamGenerator::GenerateTaggerHit(const G4Event *event, |
611 | double energy, |
612 | double time, |
613 | int bg) |
614 | { |
615 | // Register a tagger hit |
616 | |
617 | int runNo = HddmOutput::getRunNo(); |
618 | if (fTagger == 0) { |
619 | fTagger = new GlueXPseudoDetectorTAG(runNo); |
620 | } |
621 | else if (fTagger->getRunNo() != runNo) { |
622 | delete fTagger; |
623 | fTagger = new GlueXPseudoDetectorTAG(runNo); |
624 | } |
625 | return fTagger->addTaggerPhoton(event, energy, time, bg); |
626 | } |
627 | |
628 | void GlueXPhotonBeamGenerator::GenerateRFsync(const G4Event *event) |
629 | { |
630 | // Append a RF sync element to the output record, assuming that |
631 | // GenerateTriggerTime has already been called for this event. |
632 | |
633 | int runNo = HddmOutput::getRunNo(); |
634 | if (fTagger == 0) { |
635 | fTagger = new GlueXPseudoDetectorTAG(runNo); |
636 | } |
637 | else if (fTagger->getRunNo() != runNo) { |
638 | delete fTagger; |
639 | fTagger = new GlueXPseudoDetectorTAG(runNo); |
640 | } |
641 | double tsync = 512*ns * G4UniformRand()G4MTHepRandom::getTheEngine()->flat(); |
642 | tsync = fBeamBucketPeriod * floor(tsync / fBeamBucketPeriod + 0.5); |
643 | fTagger->addRFsync(event, tsync); |
644 | } |
645 | |
646 | double GlueXPhotonBeamGenerator::getBeamBucketPeriod(int runno) |
647 | { |
648 | // Look up the beam bucket period for this run in ccdb |
649 | // unless the user has already set the value by hand. |
650 | |
651 | if (runno > 0) { |
652 | jana::JCalibration *jcalib = japp->GetJCalibration(runno); |
653 | G4cout(*G4cout_p) << "JCalibration context: " << jcalib->GetContext() |
654 | << G4endlstd::endl; |
655 | std::map<std::string, double> result; |
656 | std::string map_key("/PHOTON_BEAM/RF/beam_period"); |
657 | if (jcalib->Get(map_key, result)) { |
658 | G4cerr(*G4cerr_p) << "Error in GlueXPhotonBeamGenerator::getBeamBucketPeriod" |
659 | << " - error fetching " << map_key << " from ccdb, " |
660 | << "cannot continue." << G4endlstd::endl; |
661 | exit(-1); |
662 | } |
663 | else if (result.find("beam_period") != result.end()) { |
664 | fBeamBucketPeriod = result["beam_period"] * ns; |
665 | } |
666 | else { |
667 | G4cerr(*G4cerr_p) << "Error in GlueXPhotonBeamGenerator::getBeamBucketPeriod" |
668 | << " - error finding value for " << map_key |
669 | << " in ccdb, cannot continue." << G4endlstd::endl; |
670 | exit(-1); |
671 | } |
672 | } |
673 | return fBeamBucketPeriod; |
674 | } |
675 | |
676 | double GlueXPhotonBeamGenerator::getRFreferencePlaneZ(int runno) |
677 | { |
678 | // Look up the reference plane Z for this run in ccdb |
679 | // unless the user has already set the value by hand. |
680 | |
681 | double refZ = 65 * cm; |
682 | |
683 | if (runno > 0) { |
684 | jana::JCalibration *jcalib = japp->GetJCalibration(runno); |
685 | G4cout(*G4cout_p) << "JCalibration context: " << jcalib->GetContext() |
686 | << G4endlstd::endl; |
687 | std::map<std::string, double> result; |
688 | std::string map_key("/PHOTON_BEAM/RF/reference_plane_z"); |
689 | if (jcalib->Get(map_key, result)) { |
690 | G4cerr(*G4cerr_p) << "Error in GlueXPhotonBeamGenerator::getRFreferencePlaneZ" |
691 | << " - error fetching " << map_key << " from ccdb, " |
692 | << "keeping default value " << refZ / cm << " cm." << G4endlstd::endl; |
693 | } |
694 | else if (result.find("z_position") != result.end()) { |
695 | refZ = result["z_position"] * cm; |
696 | G4cout(*G4cout_p) << "Info: RF reference plane set to " << refZ / cm << " cm." |
697 | << G4endlstd::endl; |
698 | } |
699 | else { |
700 | G4cerr(*G4cerr_p) << "Error in GlueXPhotonBeamGenerator::getRFreferencePlaneZ" |
701 | << " - error finding value for " << map_key |
702 | << " in ccdb, cannot continue." << G4endlstd::endl; |
703 | exit(-1); |
704 | } |
705 | } |
706 | return refZ; |
707 | } |