Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/GlueXPhotonBeamGenerator.cc
Location:line 107, column 11
Description:Value stored to 'raddz' during its initialization is never read

Annotated Source Code

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
29int GlueXPhotonBeamGenerator::fGenerateNotSimulate = 0;
30int GlueXPhotonBeamGenerator::fBeamBackgroundTagOnly = 0;
31
32double GlueXPhotonBeamGenerator::fBeamBucketPeriod = 4. * ns;
33double GlueXPhotonBeamGenerator::fBeamStartZ = -24 * m;
34double GlueXPhotonBeamGenerator::fBeamDiameter = 0.5 * cm;
35double GlueXPhotonBeamGenerator::fBeamVelocity = 2.99792e8 * m/s;
36double GlueXPhotonBeamGenerator::fBeamOffset[2] = {0,0};
37
38int GlueXPhotonBeamGenerator::fForceFixedPolarization = false;
39double GlueXPhotonBeamGenerator::fFixedPolarization = 0;
40double 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
45int 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
56GlueXPhotonBeamGenerator::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
128GlueXPhotonBeamGenerator::~GlueXPhotonBeamGenerator()
129{
130 delete fMessenger;
131}
132
133void 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
220void GlueXPhotonBeamGenerator::GeneratePrimaryVertex(G4Event* anEvent)
221{
222 GenerateBeamPhoton(anEvent, 0);
223}
224
225void 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
587double 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
610int 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
628void 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
646double 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
676double 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}