Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/GlueXPrimaryGeneratorAction.cc
Location:line 967, column 4
Description:Called C++ object pointer is null

Annotated Source Code

1//
2// class implementation for GlueXPrimaryGeneratorAction
3//
4// author: richard.t.jones at uconn.edu
5// version: may 12, 2012
6
7#include "GlueXPrimaryGeneratorAction.hh"
8#include "GlueXPrimaryGenerator.hh"
9#include "GlueXUserEventInformation.hh"
10#include "GlueXUserOptions.hh"
11#include "G4OpticalPhoton.hh"
12#include "HddmOutput.hh"
13
14#include "G4Event.hh"
15#include "G4ParticleTable.hh"
16#include "G4ProcessManager.hh"
17#include "G4IonTable.hh"
18#include "G4SystemOfUnits.hh"
19#include "Randomize.hh"
20
21typedef GlueXPrimaryGeneratorAction::source_type_t source_type_t;
22typedef GlueXPrimaryGeneratorAction::single_particle_gun_t particle_gun_t;
23
24source_type_t GlueXPrimaryGeneratorAction::fSourceType = SOURCE_TYPE_NONE;
25
26std::ifstream *GlueXPrimaryGeneratorAction::fHDDMinfile = 0;
27hddm_s::istream *GlueXPrimaryGeneratorAction::fHDDMistream = 0;
28
29G4ParticleTable *GlueXPrimaryGeneratorAction::fParticleTable = 0;
30
31particle_gun_t GlueXPrimaryGeneratorAction::fGunParticle;
32
33double GlueXPrimaryGeneratorAction::fBeamBackgroundRate = 0;
34double GlueXPrimaryGeneratorAction::fBeamBackgroundGateStart = 0;
35double GlueXPrimaryGeneratorAction::fBeamBackgroundGateStop = 0;
36double GlueXPrimaryGeneratorAction::fL1triggerTimeSigma = 10 * ns;
37double GlueXPrimaryGeneratorAction::fRFreferencePlaneZ = 65 * cm;
38double GlueXPrimaryGeneratorAction::fTargetCenterZ = 65 * cm;
39double GlueXPrimaryGeneratorAction::fTargetLength = 29.9746 * cm;
40double GlueXPrimaryGeneratorAction::fBeamEndpointEnergy = 12 * GeV;
41double GlueXPrimaryGeneratorAction::fBeamPeakEnergy = 9 * GeV;
42
43G4Mutex GlueXPrimaryGeneratorAction::fMutex = G4MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } };
44std::list<GlueXPrimaryGeneratorAction*> GlueXPrimaryGeneratorAction::fInstance;
45
46double GlueXPrimaryGeneratorAction::DIRC_LUT_X[48] = {0};
47double GlueXPrimaryGeneratorAction::DIRC_BAR_Y[48] = {0};
48double GlueXPrimaryGeneratorAction::DIRC_LUT_Z = 0;
49double GlueXPrimaryGeneratorAction::DIRC_QZBL_DY = 0;
50double GlueXPrimaryGeneratorAction::DIRC_QZBL_DZ = 0;
51double GlueXPrimaryGeneratorAction::DIRC_OWDG_DZ = 0;
52double GlueXPrimaryGeneratorAction::DIRC_LED_OBCS_FDTH_X = 0;
53double GlueXPrimaryGeneratorAction::DIRC_LED_OBCS_FDTH_Z = 0;
54double GlueXPrimaryGeneratorAction::DIRC_LED_OBCN_FDTH_X = 0;
55double GlueXPrimaryGeneratorAction::DIRC_LED_OBCN_FDTH_Z = 0;
56double GlueXPrimaryGeneratorAction::DIRC_LED_OBCN_FDTH1_Y = 0;
57double GlueXPrimaryGeneratorAction::DIRC_LED_OBCN_FDTH2_Y = 0;
58double GlueXPrimaryGeneratorAction::DIRC_LED_OBCN_FDTH3_Y = 0;
59double GlueXPrimaryGeneratorAction::DIRC_LED_OBCS_FDTH4_Y = 0;
60double GlueXPrimaryGeneratorAction::DIRC_LED_OBCS_FDTH5_Y = 0;
61double GlueXPrimaryGeneratorAction::DIRC_LED_OBCS_FDTH6_Y = 0;
62
63//--------------------------------------------
64// GlueXPrimaryGeneratorAction (constructor)
65//--------------------------------------------
66
67GlueXPrimaryGeneratorAction::GlueXPrimaryGeneratorAction()
68 : fParticleGun(0),
69 fCobremsGeneration(0),
70 fPhotonBeamGenerator(0),
71 fPrimaryGenerator(0),
72 fBeamvertex_activated(0)
73{
74 G4AutoLock barrier(&fMutex);
75 fInstance.push_back(this);
76
77 // Initializaton is driven by the control.in file, which
78 // gets read and parsed only once, by the first constructor.
79
80 fParticleGun = new GlueXParticleGun();
81
82 if (fSourceType == SOURCE_TYPE_HDDM) {
83 clone_photon_beam_generator();
84 fPrimaryGenerator = new GlueXPrimaryGenerator(fHDDMistream);
85 return;
86 }
87 else if (fSourceType == SOURCE_TYPE_COBREMS_GEN) {
88 clone_photon_beam_generator();
89 return;
90 }
91 else if (fSourceType == SOURCE_TYPE_PARTICLE_GUN) {
92 fParticleGun->SetParticleDefinition(fGunParticle.partDef);
93 return;
94 }
95
96 fParticleTable = G4ParticleTable::GetParticleTable();
97
98 GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance();
99 if (user_opts == 0) {
100 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction constructor - "
101 << "GlueXUserOptions::GetInstance() returns null, "
102 << "cannot continue." << G4endlstd::endl;
103 exit(-1);
104 }
105
106 // get positions for LUT from XML geometry
107 std::map<int, int> dirclutpars;
108 std::map<int, int> dircledpars;
109
110 if (user_opts->Find("DIRCLUT", dirclutpars)) {
111
112 int runno = HddmOutput::getRunNo();
113 extern jana::JApplication *japp;
114 if (japp == 0) {
115 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction constructor - "
116 << "jana global DApplication object not set, "
117 << "cannot continue." << G4endlstd::endl;
118 exit(-1);
119 }
120 jana::JGeometry *jgeom = japp->GetJGeometry(runno);
121 if (japp == 0) { // dummy
122 jgeom = 0;
123 G4cout(*G4cout_p) << "DIRC: ALL parameters loaded from ccdb" << G4endlstd::endl;
124 }
125
126 std::vector<double> DIRC;
127 std::vector<double> DRCC;
128 std::vector<double> DCML00_XYZ;
129 std::vector<double> DCML01_XYZ;
130 std::vector<double> DCML10_XYZ;
131 std::vector<double> DCML11_XYZ;
132 std::vector<double> WNGL00_XYZ;
133 std::vector<double> WNGL01_XYZ;
134 std::vector<double> WNGL10_XYZ;
135 std::vector<double> WNGL11_XYZ;
136 std::vector<double> OWDG_XYZ;
137 jgeom->Get("//section/composition/posXYZ[@volume='DIRC']/@X_Y_Z", DIRC);
138 jgeom->Get("//composition[@name='DIRC']/posXYZ[@volume='DRCC']/@X_Y_Z", DRCC);
139 jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML00']/@X_Y_Z", DCML00_XYZ);
140 jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML01']/@X_Y_Z", DCML01_XYZ);
141 jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML10']/@X_Y_Z", DCML10_XYZ);
142 jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='DCML11']/@X_Y_Z", DCML11_XYZ);
143 jgeom->Get("//composition[@name='DCML00']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL00_XYZ);
144 jgeom->Get("//composition[@name='DCML01']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL01_XYZ);
145 jgeom->Get("//composition[@name='DCML10']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL10_XYZ);
146 jgeom->Get("//composition[@name='DCML11']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL11_XYZ);
147 jgeom->Get("//composition[@name='DCML11']/posXYZ[@volume='WNGL']/@X_Y_Z", WNGL11_XYZ);
148 jgeom->Get("//trd[@name='OWDG']/@Xmp_Ymp_Z", OWDG_XYZ);
149 DIRC_LUT_Z = (DIRC[2] + DRCC[2] + DCML01_XYZ[2] + 0.8625) * cm;
150 DIRC_QZBL_DY = 3.5 * cm; // nominal width to generate LUT
151 DIRC_QZBL_DZ = 1.725 * cm; // nominal thickness to generate LUT
152 DIRC_OWDG_DZ = OWDG_XYZ[4];
153
154 // set array of bar positions
155 for (int i=0; i<48; i++) {
156 std::vector<double> DCBR_XYZ;
157 if (i<12) {
158 std::stringstream geomDCML10;
159 geomDCML10 << "//composition[@name='DCML10']/posXYZ[@volume='DCBR"
160 << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z";
161 jgeom->Get(geomDCML10.str(), DCBR_XYZ);
162 DIRC_BAR_Y[i] = (DCML10_XYZ[1] - DCBR_XYZ[1]) * cm;
163 DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML10_XYZ[0] - WNGL10_XYZ[0] + DIRC_OWDG_DZ) * cm;
164 }
165 else if (i<24) {
166 std::stringstream geomDCML11;
167 geomDCML11 << "//composition[@name='DCML11']/posXYZ[@volume='DCBR"
168 << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z";
169 jgeom->Get(geomDCML11.str(), DCBR_XYZ);
170 DIRC_BAR_Y[i] = (DCML11_XYZ[1] - DCBR_XYZ[1]) * cm;
171 DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML11_XYZ[0] - WNGL11_XYZ[0] + DIRC_OWDG_DZ) * cm;
172 }
173 else if (i<36) {
174 std::stringstream geomDCML01;
175 geomDCML01 << "//composition[@name='DCML01']/posXYZ[@volume='DCBR"
176 << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z";
177 jgeom->Get(geomDCML01.str(), DCBR_XYZ);
178 DIRC_BAR_Y[i] = (DCML01_XYZ[1] + DCBR_XYZ[1]) * cm;
179 DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML01_XYZ[0] + WNGL01_XYZ[0] - DIRC_OWDG_DZ) * cm;
180 }
181 else if (i<48) {
182 std::stringstream geomDCML00;
183 geomDCML00 << "//composition[@name='DCML00']/posXYZ[@volume='DCBR"
184 << std::setfill('0') << std::setw(2) << i << "']/@X_Y_Z";
185 jgeom->Get(geomDCML00.str(), DCBR_XYZ);
186 DIRC_BAR_Y[i] = (DCML00_XYZ[1] + DCBR_XYZ[1]) * cm;
187 DIRC_LUT_X[i] = (DIRC[0] + DRCC[0] + DCML00_XYZ[0] + WNGL00_XYZ[0] - DIRC_OWDG_DZ) * cm;
188 }
189 }
190 }
191
192 if (user_opts->Find("DIRCLED", dircledpars)) {
193 int runno = HddmOutput::getRunNo();
194 extern jana::JApplication *japp;
195 if (japp == 0) {
196 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction constructor - "
197 << "jana global DApplication object not set, "
198 << "cannot continue." << G4endlstd::endl;
199 exit(-1);
200 }
201 jana::JGeometry *jgeom = japp->GetJGeometry(runno);
202 if (japp == 0) { // dummy
203 jgeom = 0;
204 G4cout(*G4cout_p) << "DIRC: ALL parameters loaded from ccdb" << G4endlstd::endl;
205 }
206 std::vector<double> DIRC;
207 std::vector<double> DRCC;
208 std::vector<double> OBCS_XYZ;
209 std::vector<double> OBCN_XYZ;
210 std::vector<double> MRAN_XYZ;
211 std::vector<double> MRAS_XYZ;
212 std::vector<double> WM1N_XYZ;
213 std::vector<double> WM1S_XYZ;
214 std::vector<double> WM1N_BOX_XYZ;
215 std::vector<double> WM1S_BOX_XYZ;
216
217 jgeom->Get("//section/composition/posXYZ[@volume='DIRC']/@X_Y_Z", DIRC);
218 jgeom->Get("//composition[@name='DIRC']/posXYZ[@volume='DRCC']/@X_Y_Z", DRCC);
219 jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='OBCN']/@X_Y_Z", OBCN_XYZ);
220 jgeom->Get("//composition[@name='DRCC']/posXYZ[@volume='OBCS']/@X_Y_Z", OBCS_XYZ);
221 jgeom->Get("//composition[@name='OBCN']/posXYZ[@volume='MRAN']/@X_Y_Z", MRAN_XYZ);
222 jgeom->Get("//composition[@name='OBCS']/posXYZ[@volume='MRAS']/@X_Y_Z", MRAS_XYZ);
223 jgeom->Get("//composition[@name='MRAN']/posXYZ[@volume='WM1N']/@X_Y_Z", WM1N_XYZ);
224 jgeom->Get("//composition[@name='MRAS']/posXYZ[@volume='WM1S']/@X_Y_Z", WM1S_XYZ);
225 jgeom->Get("//box[@name='WM1N']/@X_Y_Z", WM1N_BOX_XYZ);
226 jgeom->Get("//box[@name='WM1S']/@X_Y_Z", WM1S_BOX_XYZ);
227
228 DIRC_LED_OBCN_FDTH_X = (DIRC[0] + DRCC[0] + OBCN_XYZ[0] + MRAN_XYZ[0] +
229 WM1N_XYZ[0] + WM1N_BOX_XYZ[0]/2. + 1.27) * cm;
230 DIRC_LED_OBCS_FDTH_X = (DIRC[0] + DRCC[0] + OBCS_XYZ[0] + MRAS_XYZ[0] +
231 WM1S_XYZ[0] - WM1S_BOX_XYZ[0]/2. - 1.27) * cm;
232
233 DIRC_LED_OBCN_FDTH_Z = (DIRC[2] + DRCC[2] + OBCN_XYZ[2] + MRAN_XYZ[2] +
234 WM1N_XYZ[2] - WM1N_BOX_XYZ[2]/2. ) * cm;
235 DIRC_LED_OBCS_FDTH_Z = (DIRC[2] + DRCC[2] + OBCS_XYZ[2] + MRAS_XYZ[2] +
236 WM1S_XYZ[2] - WM1S_BOX_XYZ[2]/2. ) * cm;
237
238 DIRC_LED_OBCN_FDTH1_Y = (DIRC[1] + DRCC[1] + OBCN_XYZ[1] + MRAN_XYZ[1] +
239 WM1N_XYZ[1] - WM1N_BOX_XYZ[1]/2. + 17.235932) * cm;
240 DIRC_LED_OBCN_FDTH2_Y = (DIRC[1] + DRCC[1] + OBCN_XYZ[1] + MRAN_XYZ[1] +
241 WM1N_XYZ[1] - WM1N_BOX_XYZ[1]/2. + 17.235932 + 31.800038 ) * cm;
242 DIRC_LED_OBCN_FDTH3_Y = (DIRC[1] + DRCC[1] + OBCN_XYZ[1] + MRAN_XYZ[1] +
243 WM1N_XYZ[1] - WM1N_BOX_XYZ[1]/2. + 17.235932 + 31.800038 * 2. ) * cm;
244 DIRC_LED_OBCS_FDTH4_Y = (DIRC[1] + DRCC[1] + OBCS_XYZ[1] + MRAS_XYZ[1] +
245 WM1S_XYZ[1] + WM1S_BOX_XYZ[1]/2. - 17.235932) * cm;
246 DIRC_LED_OBCS_FDTH5_Y = (DIRC[1] + DRCC[1] + OBCS_XYZ[1] + MRAS_XYZ[1] +
247 WM1S_XYZ[1] + WM1S_BOX_XYZ[1]/2. - 17.235932 - 31.800038 ) * cm;
248 DIRC_LED_OBCS_FDTH6_Y = (DIRC[1] + DRCC[1] + OBCS_XYZ[1] + MRAS_XYZ[1] +
249 WM1S_XYZ[1] + WM1S_BOX_XYZ[1]/2. - 17.235932 - 31.800038 * 2. ) * cm;
250
251 }
252
253 std::map<int,std::string> infile;
254 std::map<int,double> beampars;
255 std::map<int,double> genbeampars;
256 std::map<int,double> kinepars;
257
258 // Three event source options are supported:
259 // 1) external generator, hddm input stream source
260 // 2) internal coherent bremsstrahlung beam generator
261 // 3) internal particle gun generator
262
263 if (user_opts->Find("INFILE", infile) ||
264 user_opts->Find("INFI", infile))
265 {
266 fHDDMinfile = new std::ifstream(infile[1].c_str());
267 if (!fHDDMinfile->is_open()) {
268 G4cerr(*G4cerr_p) << "GlueXPrimaryGeneratorAction error: "
269 << "Unable to open HDDM input file: " << infile[1]
270 << G4endlstd::endl;
271 exit(-1);
272 }
273 fHDDMistream = new hddm_s::istream(*fHDDMinfile);
274 G4cout(*G4cout_p) << "Opened input file: " << infile[1] << G4endlstd::endl;
275 std::map<int,int> skippars;
276 if (user_opts->Find("SKIP", skippars))
277 {
278 if (skippars[1] > 0)
279 {
280 fHDDMistream->skip(skippars[1]);
281 G4cout(*G4cout_p) << "skipped first " << skippars[1] << " input events." << G4endlstd::endl;
282 }
283 }
284 fPrimaryGenerator = new GlueXPrimaryGenerator(fHDDMistream);
285 fSourceType = SOURCE_TYPE_HDDM;
286 }
287
288 else if (user_opts->Find("BEAM", beampars) ||
289 user_opts->Find("GENBEAM", genbeampars))
290 {
291 fSourceType = SOURCE_TYPE_COBREMS_GEN;
292 }
293
294 else if (user_opts->Find("DIRCLUT", dirclutpars))
295 {
296 fGunParticle.geantType = 0;
297 fGunParticle.pdgType = 999999;
298 fGunParticle.partDef = fParticleTable->FindParticle("opticalphoton");
299 fGunParticle.deltaR = 0;
300 fGunParticle.deltaZ = 0;
301 fGunParticle.mom = 3.18 * eV;
302
303 fGunParticle.deltaMom = 0;
304 fGunParticle.deltaTheta = 0;
305 fGunParticle.deltaPhi = 0;
306 fParticleGun->SetParticleDefinition(fGunParticle.partDef);
307
308 fSourceType = SOURCE_TYPE_PARTICLE_GUN;
309 }
310
311 else if (user_opts->Find("DIRCLED", dircledpars))
312 {
313 fGunParticle.geantType = 0;
314 fGunParticle.pdgType = 999999;
315 fGunParticle.partDef = fParticleTable->FindParticle("opticalphoton");
316 fGunParticle.deltaR = 0;
317 fGunParticle.deltaZ = 0;
318 fGunParticle.mom = 3.0613 * eV;
319
320 fGunParticle.deltaMom = 0;
321 fGunParticle.deltaTheta = 0;
322 fGunParticle.deltaPhi = 0;
323 fParticleGun->SetParticleDefinition(fGunParticle.partDef);
324
325 fSourceType = SOURCE_TYPE_PARTICLE_GUN;
326 }
327
328 else if (user_opts->Find("KINE", kinepars))
329 {
330 if (kinepars[1] == 1000 || kinepars[1] == 148 || kinepars[1] == 48) {
331 fGunParticle.geantType = 0;
332 fGunParticle.pdgType = 999999;
333 fGunParticle.partDef = fParticleTable->FindParticle("geantino");
334 }
335 else if (kinepars[1] == 1001) {
336 fGunParticle.geantType = 0;
337 fGunParticle.pdgType = 999999;
338 fGunParticle.partDef = fParticleTable->FindParticle("chargedgeantino");
339 }
340 else if (kinepars[1] == 1050) {
341 fGunParticle.geantType = 0;
342 fGunParticle.pdgType = 999999;
343 fGunParticle.partDef = fParticleTable->FindParticle("opticalphoton");
344 }
345 else {
346 if (kinepars[1] > 100)
347 fGunParticle.geantType = kinepars[1] - 100;
348 else
349 fGunParticle.geantType = kinepars[1];
350 fGunParticle.pdgType = ConvertGeant3ToPdg(fGunParticle.geantType);
351 fGunParticle.partDef = fParticleTable->FindParticle(fGunParticle.pdgType);
352 }
353 if (fGunParticle.partDef == 0) {
354 G4cerr(*G4cerr_p) << "GlueXPrimaryGeneratorAction constructor error - "
355 << "Unknown GEANT particle type: " << kinepars[1]
356 << " was specified in the control.in file." << G4endlstd::endl;
357 exit(-1);
358 }
359 fParticleGun->SetParticleDefinition(fGunParticle.partDef);
360
361 double x(0), y(0), z(65 * cm);
362 std::map<int,double> scappars;
363 if (user_opts->Find("SCAP", scappars)) {
364 x = scappars[1] * cm;
365 y = scappars[2] * cm;
366 z = scappars[3] * cm;
367 }
368 fGunParticle.pos.set(x,y,z);
369 std::map<int,double> tgtwidthpars;
370 if (user_opts->Find("tgtwidth", tgtwidthpars)) {
371 fGunParticle.deltaR = tgtwidthpars[1] * cm;
372 fGunParticle.deltaZ = tgtwidthpars[2] * cm;
373 }
374 else {
375 fGunParticle.deltaR = 0;
376 fGunParticle.deltaZ = 0;
377 }
378 fGunParticle.plogOption = 0;
379 std::map<int,int> plogpars;
380 if (user_opts->Find("PLOG", plogpars)) {
381 fGunParticle.plogOption = plogpars[1];
382 }
383 fGunParticle.tlogOption = 0;
384 std::map<int,int> tlogpars;
385 if (user_opts->Find("TLOG", tlogpars)) {
386 fGunParticle.tlogOption = tlogpars[1];
387 }
388
389 fGunParticle.mom = kinepars[2] * GeV;
390 if (kinepars[1] > 100) {
391 fGunParticle.theta = kinepars[3] * degree;
392 fGunParticle.phi = kinepars[4] * degree;
393 fGunParticle.deltaMom = kinepars[5] * GeV;
394 fGunParticle.deltaTheta = kinepars[6] * degree;
395 fGunParticle.deltaPhi = kinepars[7] * degree;
396 }
397 else {
398 fGunParticle.deltaMom = 0;
399 fGunParticle.theta = 90 * degree;
400 fGunParticle.deltaTheta = 180 * degree;
401 fGunParticle.phi = 0;
402 fGunParticle.deltaPhi = 360 * degree;
403 }
404 fSourceType = SOURCE_TYPE_PARTICLE_GUN;
405 }
406
407 if (user_opts->Find("BEAM", beampars)) {
408 double beamE0 = beampars[1] * GeV;
409 double beamEpeak = beampars[2] * GeV;
410 double beamEmin = ((beampars[3] > 0)? beampars[3] : 0.120) * GeV;
411 double radColDist = ((beampars[4] > 0)? beampars[4] : 76.) * m;
412 double colDiam = ((beampars[5] > 0)? beampars[5] : 0.0034) * m;
413 double beamEmit = ((beampars[6] > 0)? beampars[6] : 2.5e-9) * m;
414 double radThick = ((beampars[7] > 0)? beampars[7] : 20e-6) * m;
415 double spotRMS = ((beampars[8] > 0)? beampars[8] : 5e-4) * m;
416 double spotX = ((beampars[9] != 0)? beampars[9] : 0) * m;
417 double spotY = ((beampars[10] != 0)? beampars[10] : 0) * m;
418
419 if (beamE0 == 0) {
420 G4cerr(*G4cerr_p) << "GlueXPrimaryGeneratorAction error: "
421 << "BEAM card specified in control.in but required values "
422 << "Emax and/or Epeak are missing, cannot continue."
423 << G4endlstd::endl;
424 exit(-1);
425 }
426
427 fBeamEndpointEnergy = beamE0;
428 fBeamPeakEnergy = beamEpeak;
429
430 // CobremsGeneration has its own standard units that it uses:
431 // length : m
432 // angles : radians
433 // energy : GeV
434 // time : s
435 // current: microAmps
436
437 fCobremsGeneration = new CobremsGeneration(beamE0/GeV, beamEpeak/GeV);
438 fCobremsGeneration->setPhotonEnergyMin(beamEmin/GeV);
439 fCobremsGeneration->setCollimatorDistance(radColDist/m);
440 fCobremsGeneration->setCollimatorDiameter(colDiam/m);
441 fCobremsGeneration->setBeamEmittance(beamEmit/(m*radian));
442 fCobremsGeneration->setTargetThickness(radThick/m);
443 fCobremsGeneration->setCollimatorSpotrms(spotRMS/m);
444 fPhotonBeamGenerator = new GlueXPhotonBeamGenerator(fCobremsGeneration);
445 fPhotonBeamGenerator->setBeamOffset(spotX, spotY);
446 }
447 else { // look up beam properties in the ccdb based on run number
448 int runno = HddmOutput::getRunNo();
449 extern jana::JApplication *japp;
450 if (japp == 0) {
451 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction constructor - "
452 << "jana global DApplication object not set, "
453 << "cannot continue." << G4endlstd::endl;
454 exit(-1);
455 }
456 std::map<string, float> beam_parms;
457 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
458 jcalib->Get("PHOTON_BEAM/endpoint_energy", beam_parms);
459 fBeamEndpointEnergy = beam_parms.at("PHOTON_BEAM_ENDPOINT_ENERGY")*GeV;
460 std::map<string, float> coherent_parms;
461 jcalib->Get("PHOTON_BEAM/coherent_energy", coherent_parms);
462 fBeamPeakEnergy = coherent_parms.at("cohedge_energy")*GeV;
463
464 // CobremsGeneration has its own standard units that it uses:
465 // length : m
466 // angles : radians
467 // energy : GeV
468 // time : s
469 // current: microAmps
470
471 fCobremsGeneration = new CobremsGeneration(fBeamEndpointEnergy/GeV,
472 fBeamPeakEnergy/GeV);
473 fCobremsGeneration->setPhotonEnergyMin(fBeamEndpointEnergy/GeV * 0.25);
474 fCobremsGeneration->setCollimatorDistance(76.);
475 fCobremsGeneration->setCollimatorDiameter(0.005);
476 fCobremsGeneration->setBeamEmittance(2e-9);
477 fCobremsGeneration->setTargetThickness(50e-6);
478 fCobremsGeneration->setCollimatorSpotrms(0.5e-3);
479 fPhotonBeamGenerator = new GlueXPhotonBeamGenerator(fCobremsGeneration);
480 std::map<string, float> beam_spot_params;
481 jcalib->Get("/PHOTON_BEAM/beam_spot", beam_spot_params);
482 double spotX = beam_spot_params.at("x") * cm;
483 double spotY = beam_spot_params.at("y") * cm;
484 fPhotonBeamGenerator->setBeamOffset(spotX, spotY);
485 std::cout << "no beam card found, looking up beam properties in ccdb, "
486 << "run number " << runno
487 << ": endpoint=" << fBeamEndpointEnergy/GeV << "GeV,"
488 << " peak=" << fBeamPeakEnergy/GeV << "GeV" << std::endl;
489 }
490
491 std::map<int, double> bgratepars;
492 std::map<int, double> bggatepars;
493 if (user_opts->Find("BGRATE", bgratepars) &&
494 user_opts->Find("BGGATE", bggatepars))
495 {
496 fBeamBackgroundRate = bgratepars[1] * 1/ns;
497 fBeamBackgroundGateStart = bggatepars[1] * ns;
498 fBeamBackgroundGateStop = bggatepars[2] * ns;
499 if (fBeamBackgroundRate > 0 &&
500 fBeamBackgroundGateStart >= fBeamBackgroundGateStop)
501 {
502 G4cerr(*G4cerr_p) << "GlueXPrimaryGeneratorAction error: "
503 << "BGRATE is non-zero, but the time window specified "
504 << "in BGGATE is invalid."
505 << G4endlstd::endl;
506 exit(-1);
507 }
508 }
509
510 std::map<int,double> trefsigma;
511 if (user_opts->Find("trefsigma", trefsigma)) {
512 fL1triggerTimeSigma = trefsigma[1] * ns;
513 }
514 else {
515 fL1triggerTimeSigma = 10 * ns;
516 }
517}
518
519GlueXPrimaryGeneratorAction::GlueXPrimaryGeneratorAction(const
520 GlueXPrimaryGeneratorAction &src)
521 : GlueXPrimaryGeneratorAction()
522{
523 fBeamvertex = src.fBeamvertex;
524 fBeamvertex_activated = src.fBeamvertex_activated;
525}
526
527GlueXPrimaryGeneratorAction &GlueXPrimaryGeneratorAction::operator=(const
528 GlueXPrimaryGeneratorAction &src)
529{
530 fParticleGun = new GlueXParticleGun();
531 fCobremsGeneration = 0;
532 fPhotonBeamGenerator = 0;
533 fPrimaryGenerator = 0;
534
535 if (fSourceType == SOURCE_TYPE_HDDM) {
536 fPrimaryGenerator = new GlueXPrimaryGenerator(fHDDMistream);
537 }
538 else if (fSourceType == SOURCE_TYPE_PARTICLE_GUN) {
539 fParticleGun->SetParticleDefinition(fGunParticle.partDef);
540 }
541 if (src.fPhotonBeamGenerator) {
542 clone_photon_beam_generator();
543 }
544 fBeamvertex = src.fBeamvertex;
545 fBeamvertex_activated = src.fBeamvertex_activated;
546 return *this;
547}
548
549//--------------------------------------------
550// ~GlueXPrimaryGeneratorAction (destructor)
551//--------------------------------------------
552
553GlueXPrimaryGeneratorAction::~GlueXPrimaryGeneratorAction()
554{
555 G4AutoLock barrier(&fMutex);
556 fInstance.remove(this);
557 if (fPrimaryGenerator)
558 delete fPrimaryGenerator;
559 if (fCobremsGeneration)
560 delete fCobremsGeneration;
561 if (fPhotonBeamGenerator)
562 delete fPhotonBeamGenerator;
563 delete fParticleGun;
564 if (fInstance.size() == 0) {
565 if (fHDDMistream)
566 delete fHDDMistream;
567 if (fHDDMinfile)
568 delete fHDDMinfile;
569 }
570}
571
572const CobremsGeneration *GlueXPrimaryGeneratorAction::GetCobremsGeneration() const
573{
574 return fCobremsGeneration;
575}
576
577//--------------------------------------------
578// GeneratePrimaries
579//--------------------------------------------
580
581void GlueXPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
582{
583 switch (fSourceType) {
584 case SOURCE_TYPE_HDDM:
585 GeneratePrimariesHDDM(anEvent);
586 break;
587 case SOURCE_TYPE_COBREMS_GEN:
588 GeneratePrimariesCobrems(anEvent);
589 break;
590 case SOURCE_TYPE_PARTICLE_GUN:
591 GeneratePrimariesParticleGun(anEvent);
592 break;
593 default:
594 G4cout(*G4cout_p) << "No event source selected, cannot continue!" << G4endlstd::endl;
595 exit(-1);
596 }
597}
598
599//--------------------------------------------
600// GeneratePrimariesParticleGun
601//--------------------------------------------
602
603void GlueXPrimaryGeneratorAction::GeneratePrimariesParticleGun(G4Event* anEvent)
604{
605 // GlueXUserEventInformation constructor can seed the random number
606 // generator for this event, so this must happen here at the top.
607 GlueXUserEventInformation *event_info = new GlueXUserEventInformation();
608 anEvent->SetUserInformation(event_info);
609 GlueXUserEventInformation::setWriteNoHitEvents(1);
610
611 // Unbelievably, GEANT4's G4ParticleGun class insists on printing
612 // a message whenever the momentum or energy is changed, unless
613 // the other is 0. Here, we reset the particle gun energy using
614 // our own derived class. (Sheesh!!)
615 fParticleGun->Reset();
616
617 // std::cout<<"GlueXPrimaryGeneratorAction:: GENERATE PRIMARIES PARTICLE GUN"<<std::endl;
618
619 GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance();
620 std::map<int,int> dirclutpars;
621 std::map<int,int> dircledpars;
622
623 // place and smear the particle gun origin
624 G4ThreeVector pos(fGunParticle.pos);
625 if (fGunParticle.deltaR > 0) {
626 double dx, dy;
627 while (true) {
628 double rx = G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5;
629 double ry = G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5;
630 if (rx*rx + ry*ry <= 0.25) {
631 dx = rx * 2 * fGunParticle.deltaR;
632 dy = ry * 2 * fGunParticle.deltaR;
633 break;
634 }
635 }
636 pos += G4ThreeVector(dx, dy, 0);
637 }
638 if (fGunParticle.deltaZ > 0) {
639 double dz = (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5) * fGunParticle.deltaZ;
640 pos += G4ThreeVector(0, 0, dz);
641 }
642 fParticleGun->SetParticlePosition(pos);
643
644 // Assign and optionally smear the particle momentum
645 double p = fGunParticle.mom;
646 double thetap = fGunParticle.theta;
647 double phip = fGunParticle.phi;
648 if (fGunParticle.deltaMom > 0) {
649 if (fGunParticle.plogOption) {
650 double pmin = p - fGunParticle.deltaMom / 2;
651 double pmax = p + fGunParticle.deltaMom / 2;
652 pmin = (pmin > 0)? pmin : 1e-6*GeV;
653 p = pmin * pow(pmax/pmin, G4UniformRand()G4MTHepRandom::getTheEngine()->flat());
654 }
655 else {
656 p += (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5) * fGunParticle.deltaMom;
657 }
658 }
659 if (fGunParticle.deltaTheta > 0) {
660 if (fGunParticle.plogOption) {
661 double thetamin = thetap - fGunParticle.deltaTheta / 2;
662 double thetamax = thetap + fGunParticle.deltaTheta / 2;
663 thetamin = (thetamin > 0)? thetamin : 1e-6*degree;
664 thetap = thetamin * pow(thetamax/thetamin, G4UniformRand()G4MTHepRandom::getTheEngine()->flat());
665 }
666 else {
667 thetap += (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5) * fGunParticle.deltaTheta;
668 }
669 }
670 if (fGunParticle.deltaPhi > 0)
671 phip += (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5) * fGunParticle.deltaPhi;
672
673 // Special case of Cherenkov photon gun for DIRC Look Up Tables (LUT)
674 if (user_opts->Find("DIRCLUT", dirclutpars)) {
675
676 // array of bar y-positions for LUT from JGeometry
677 double y = 0.; // no shift
678 double x = DIRC_LUT_X[dirclutpars[1]];
679 double z = DIRC_LUT_Z;
680
681 G4ThreeVector vec(0,0,1);
682 double rand1 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
683 double rand2 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
684 vec.setTheta(acos(rand1));
685 vec.setPhi(2*M_PI3.14159265358979323846*rand2);
686 vec.rotateY(M_PI3.14159265358979323846/2.);
687 y = DIRC_BAR_Y[dirclutpars[1]];
688 if (dirclutpars[1] < 24) {
689 vec.rotateY(M_PI3.14159265358979323846);
690 }
691
692 // spread over end of bar in y and z
693 y += DIRC_QZBL_DY/2.0 - DIRC_QZBL_DY*G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
694 z += DIRC_QZBL_DZ/2.0 - DIRC_QZBL_DZ*G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
695
696 thetap = vec.theta();
697 phip = vec.phi();
698 fParticleGun->SetParticlePosition(G4ThreeVector(x,y,z));
699 }
700
701
702
703 double DeltaT = 0.;
704 // Special case of Cherenkov photon gun for DIRC LED generator
705 if (user_opts->Find("DIRCLED", dircledpars)){
706
707 double x(0.),y(0.),z(0.);
708
709 int FDTH = -1;
710 vector<int> FDTHs = {};
711 for (int par_index = 1; par_index <= 6; par_index++)
712 {
713 int passed_FDTH = dircledpars[par_index];
714 if (passed_FDTH && 0 < passed_FDTH && passed_FDTH < 7)
715 FDTHs.push_back(dircledpars[par_index]);
716 }
717 int NumFDTHs = int(FDTHs.size());
718 double rand0 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
719 for (int FDTH_index = 0; FDTH_index < NumFDTHs; FDTH_index++)
720 {
721 if (rand0 <= (FDTH_index+1)*(1./NumFDTHs))
722 {
723 FDTH = FDTHs[FDTH_index]; break;
724 }
725 else
726 continue;
727 }
728
729 switch (FDTH)
730 {
731 case 1:
732 x = DIRC_LED_OBCN_FDTH_X;
733 y = DIRC_LED_OBCN_FDTH1_Y;
734 z = DIRC_LED_OBCN_FDTH_Z;
735 break;
736
737 case 2:
738 x = DIRC_LED_OBCN_FDTH_X;
739 y = DIRC_LED_OBCN_FDTH2_Y;
740 z = DIRC_LED_OBCN_FDTH_Z;
741 break;
742
743 case 3:
744 x = DIRC_LED_OBCN_FDTH_X;
745 y = DIRC_LED_OBCN_FDTH3_Y;
746 z = DIRC_LED_OBCN_FDTH_Z;
747 break;
748
749 case 4:
750 x = DIRC_LED_OBCS_FDTH_X;
751 y = DIRC_LED_OBCS_FDTH4_Y;
752 z = DIRC_LED_OBCS_FDTH_Z;
753 break;
754
755 case 5:
756 x = DIRC_LED_OBCS_FDTH_X;
757 y = DIRC_LED_OBCS_FDTH5_Y;
758 z = DIRC_LED_OBCS_FDTH_Z;
759 break;
760
761 case 6:
762 x = DIRC_LED_OBCS_FDTH_X;
763 y = DIRC_LED_OBCS_FDTH6_Y;
764 z = DIRC_LED_OBCS_FDTH_Z;
765 break;
766
767 default:
768 break;
769 }
770
771 //z -= 0.5*cm;
772 double theta_range = 12.5; // in degrees
773 double inclination_wrt_bars_deg = 0.;//negative->away from 3-segment mirror; positive->towards 3-segment mirror
774 double angle_towards_center_deg = 0.;//bending angle for the two feedthroughs on the sides towards the center
775
776 if (x > 0.)
777 inclination_wrt_bars_deg *= -1.;
778 if (FDTH == 3 || FDTH == 4)
779 angle_towards_center_deg *= -1.;
780 if (FDTH == 2 || FDTH == 5)
781 angle_towards_center_deg = 0.;
782
783 G4ThreeVector vec(0,0,1);
784 double rand1 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
785 double rand2 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
786
787 double costheta = -1. + rand1 * (std::cos((180.-theta_range) * M_PI3.14159265358979323846 / 180.) + 1.);
788
789 double rand3 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
790 double theta_to_set = acos(costheta);
791 if (rand3 < 0.5)
792 theta_to_set = 2*M_PI3.14159265358979323846 - theta_to_set;
793
794 vec.setTheta(theta_to_set);
795 vec.setPhi(2*M_PI3.14159265358979323846*rand2);
796 vec.rotateY(inclination_wrt_bars_deg*deg);
797 vec.rotateX(angle_towards_center_deg*deg);
798/*
799 //For square diffuser case
800 double theta_range = 25.; // in degrees
801 double diffuser_X = 1.697 * cm;
802 double diffuser_Y = 1.697 * cm;
803 double inclination_wrt_bars_deg = -6.;
804
805 G4ThreeVector vec(0,0,-1.);
806 double rand1 = G4UniformRand();
807 double rand2 = G4UniformRand();
808 double a = 2.* tan(theta_range*M_PI/180.);
809 vec.setX((rand1-0.5)*a);
810 vec.setY((rand2-0.5)*a);
811 vec.setZ(-1.);
812
813 double diffuser_offset_x = diffuser_X/2. - diffuser_X * G4UniformRand();
814 double diffuser_offset_y = diffuser_Y/2. - diffuser_Y * G4UniformRand();
815
816 G4ThreeVector diffuser_offset_vec(diffuser_offset_x,diffuser_offset_y,0.);
817 diffuser_offset_vec.rotateY(inclination_wrt_bars_deg*deg);
818
819 x += diffuser_offset_vec.x();
820 y += diffuser_offset_vec.y();
821 z += diffuser_offset_vec.z();
822
823*/
824
825 // time smearing from LED pulse shape
826 double t_rise = 0.84;
827 double t_FWHM = 1.5;
828
829 double t_range_LED = t_rise + t_FWHM;
830 double t_rand1 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
831 double t_rand2 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
832 double t_LEDShape = t_rand1 * t_range_LED;
833
834 double pulse_val = -1.;
835 while (t_rand2 > pulse_val)
836 {
837 t_rand1 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
838 t_rand2 = G4UniformRand()G4MTHepRandom::getTheEngine()->flat();
839 t_LEDShape = t_rand1 * t_range_LED;
840
841 if (0.<= t_LEDShape && t_LEDShape < t_rise)
842 pulse_val = 1./t_rise * t_LEDShape;
843 else if (t_rise <= t_LEDShape && t_LEDShape < t_FWHM)
844 pulse_val = 1.;
845 else if (t_FWHM <= t_LEDShape && t_LEDShape < t_range_LED)
846 pulse_val = -1./t_rise * t_LEDShape + 1./t_rise * t_range_LED;
847 else
848 pulse_val = 0.;
849
850 if (t_rand2 <= pulse_val)
851 break;
852 }
853 double DeltaT_LED = t_LEDShape;
854
855
856 // delay from LED feedthrough cables
857 double DeltaT_CableDelay = 0.;
858 if (FDTH == 2 || FDTH == 5)
859 DeltaT_CableDelay = 10.;
860 if (FDTH == 3 || FDTH == 6)
861 DeltaT_CableDelay = 20.;
862
863
864 DeltaT = DeltaT_LED + DeltaT_CableDelay;
865
866 G4ThreeVector pos_vec(x,y,z);
867 thetap = vec.theta();
868 phip = vec.phi();
869 fParticleGun->SetParticlePosition(pos_vec);
870
871 }
872
873 G4ThreeVector mom(p * sin(thetap) * cos(phip),
874 p * sin(thetap) * sin(phip),
875 p * cos(thetap));
876 fParticleGun->SetParticleMomentum(mom);
877 fParticleGun->SetParticleTime(0.+DeltaT);
878
879 // Fire the gun
880 fParticleGun->GeneratePrimaryVertex(anEvent);
881
882 // Store generated particle info so it can be written to output file
883 event_info->AddPrimaryVertex(*anEvent->GetPrimaryVertex());
884}
885
886//--------------------------------------------
887// GeneratePrimariesHDDM
888//--------------------------------------------
889
890void GlueXPrimaryGeneratorAction::GeneratePrimariesHDDM(G4Event* anEvent)
891{
892 if (fPrimaryGenerator != 0) {
1
Taking false branch
→
893 double beamDiameter = GlueXPhotonBeamGenerator::getBeamDiameter();
894 double beamVelocity = GlueXPhotonBeamGenerator::getBeamVelocity();
895 double x, y, z;
896 if (fBeamvertex_activated == 0)
897 configure_beam_vertex();
898 if (fBeamvertex_activated == 1) {
899
900 // generate a random beam spot according to a gaussian model
901
902 double ur = sqrt(-2 * log(G4UniformRand()G4MTHepRandom::getTheEngine()->flat()));
903 double phi = G4UniformRand()G4MTHepRandom::getTheEngine()->flat() * 2*M_PI3.14159265358979323846;
904 double uz = G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5;
905 double u[2] = {ur * cos(phi), ur * sin(phi)};
906 u[0] *= fBeamvertex.sigma[0];
907 u[1] *= fBeamvertex.sigma[1];
908 x = u[0] * cos(fBeamvertex.alpha) + u[1] * sin(fBeamvertex.alpha);
909 y =-u[0] * sin(fBeamvertex.alpha) + u[1] * cos(fBeamvertex.alpha);
910 z = uz * fBeamvertex.length;
911 x += fBeamvertex.x + fBeamvertex.dxdz * z;
912 y += fBeamvertex.y + fBeamvertex.dydz * z;
913 z += fBeamvertex.z;
914 }
915 else {
916
917 // fall back on the old hard-coded cylindrical beam spot
918
919 while (true) {
920 x = G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5;
921 y = G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5;
922 if (x*x + y*y <= 0.25) {
923 x *= beamDiameter;
924 y *= beamDiameter;
925 break;
926 }
927 }
928 z = fTargetCenterZ + (G4UniformRand()G4MTHepRandom::getTheEngine()->flat() - 0.5) * fTargetLength;
929 }
930 assert (fPrimaryGenerator != 0)((fPrimaryGenerator != 0) ? static_cast<void> (0) : __assert_fail
("fPrimaryGenerator != 0", "src/GlueXPrimaryGeneratorAction.cc"
, 930, __PRETTY_FUNCTION__));
931 assert (fPhotonBeamGenerator != 0)((fPhotonBeamGenerator != 0) ? static_cast<void> (0) : __assert_fail
("fPhotonBeamGenerator != 0", "src/GlueXPrimaryGeneratorAction.cc"
, 931, __PRETTY_FUNCTION__));
932 double ttag = fPhotonBeamGenerator->GenerateTriggerTime(anEvent);
933 double trel = (z - fRFreferencePlaneZ) / beamVelocity;
934 fPrimaryGenerator->SetParticlePosition(G4ThreeVector(x,y,z));
935 fPrimaryGenerator->SetParticleTime(trel + ttag);
936 fPrimaryGenerator->GeneratePrimaryVertex(anEvent);
937 GlueXUserEventInformation *eventinfo;
938 eventinfo = (GlueXUserEventInformation*)anEvent->GetUserInformation();
939 if (eventinfo) {
940 // The above-assigned vertex coordinates and time are advisory to
941 // GeneratePrimaryVertex, and may have been overridden by values
942 // read from the input MC event record.
943 G4PrimaryVertex *vertex = anEvent->GetPrimaryVertex();
944 trel = (vertex->GetZ0() - fRFreferencePlaneZ) / beamVelocity;
945 ttag = vertex->GetT0() - trel;
946 double E = eventinfo->GetBeamPhotonEnergy();
947 fPhotonBeamGenerator->GenerateTaggerHit(anEvent, E, ttag);
948 }
949 else {
950 return;
951 }
952 }
953
954 // Superimpose any requested background minimum-bias beam interactions
955
956 if (fBeamBackgroundRate > 0 && fPhotonBeamGenerator != 0) {
2
←
Assuming pointer value is null
→
3
←
Taking false branch
→
957 double t = fBeamBackgroundGateStart;
958 while (true) {
959 t += -log(G4UniformRand()G4MTHepRandom::getTheEngine()->flat()) / fBeamBackgroundRate;
960 if (t > fBeamBackgroundGateStop)
961 break;
962 fPhotonBeamGenerator->GenerateBeamPhoton(anEvent, t);
963 }
964 }
965
966 // Add the RF sync
967 fPhotonBeamGenerator->GenerateRFsync(anEvent);
4
←
Called C++ object pointer is null
968}
969
970void GlueXPrimaryGeneratorAction::GeneratePrimariesCobrems(G4Event* anEvent)
971{
972 if (fPhotonBeamGenerator != 0) {
973 fPhotonBeamGenerator->GeneratePrimaryVertex(anEvent);
974 }
975}
976
977// Convert particle types from Geant3 types to PDG scheme
978
979int GlueXPrimaryGeneratorAction::ConvertGeant3ToPdg(int Geant3type)
980{
981 // This method was imported from ROOT source file TDatabasePDG.cc
982
983 switch (Geant3type) {
984 case 1 : return 22; // photon
985 case 2 : return -11; // e+
986 case 3 : return 11; // e-
987 case 4 : return 12; // e-neutrino (NB: flavour undefined by Geant)
988 case 5 : return -13; // mu+
989 case 6 : return 13; // mu-
990 case 7 : return 111; // pi0
991 case 8 : return 211; // pi+
992 case 9 : return -211; // pi-
993 case 10 : return 130; // K long
994 case 11 : return 321; // K+
995 case 12 : return -321; // K-
996 case 13 : return 2112; // n
997 case 14 : return 2212; // p
998 case 15 : return -2212; // anti-proton
999 case 16 : return 310; // K short
1000 case 17 : return 221; // eta
1001 case 18 : return 3122; // Lambda
1002 case 19 : return 3222; // Sigma+
1003 case 20 : return 3212; // Sigma0
1004 case 21 : return 3112; // Sigma-
1005 case 22 : return 3322; // Xi0
1006 case 23 : return 3312; // Xi-
1007 case 24 : return 3334; // Omega- (PB)
1008 case 25 : return -2112; // anti-neutron
1009 case 26 : return -3122; // anti-Lambda
1010 case 27 : return -3222; // anti-Sigma-
1011 case 28 : return -3212; // anti-Sigma0
1012 case 29 : return -3112; // anti-Sigma+
1013 case 30 : return -3322; // anti-Xi0
1014 case 31 : return -3312; // anti-Xi+
1015 case 32 : return -3334; // anti-Omega+
1016 case 45 : return 1000010020; // deuteron
1017 case 46 : return 1000010030; // triton
1018 case 47 : return 1000020040; // alpha
1019 case 48 : return 0; // geantino (no PDG type)
1020 case 49 : return 1000020030; // He3 ion
1021 case 50 : return 0; // Cerenkov photon (no PDG type)
1022 case 61 : return 1000030060; // Li6
1023 case 62 : return 1000030070; // Li7
1024 case 63 : return 1000040070; // Be7
1025 case 64 : return 1000040090; // Be9
1026 case 65 : return 1000050100; // B10
1027 case 66 : return 1000050110; // B11
1028 case 67 : return 1000060120; // C12
1029 case 68 : return 1000070140; // N14
1030 case 69 : return 1000080160; // O16
1031 case 70 : return 1000090190; // F19
1032 case 71 : return 1000100200; // Ne20
1033 case 72 : return 1000110230; // Na23
1034 case 73 : return 1000120240; // Mg24
1035 case 74 : return 1000130270; // Al27
1036 case 75 : return 1000140280; // Si28
1037 case 76 : return 1000150310; // P31
1038 case 77 : return 1000160320; // S32
1039 case 78 : return 1000170350; // Cl35
1040 case 79 : return 1000180360; // Ar36
1041 case 80 : return 1000190390; // K39
1042 case 81 : return 1000200400; // Ca40
1043 case 82 : return 1000210450; // Sc45
1044 case 83 : return 1000220480; // Ti48
1045 case 84 : return 1000230510; // V51
1046 case 85 : return 1000240520; // Cr52
1047 case 86 : return 1000250550; // Mn55
1048 case 87 : return 1000260560; // Fe56
1049 case 88 : return 1000270590; // Co59
1050 case 89 : return 1000280580; // Ni58
1051 case 90 : return 1000290630; // Cu63
1052 case 91 : return 1000300640; // Zn64
1053 case 92 : return 1000320740; // Ge74
1054 case 93 : return 1000340800; // Se80
1055 case 94 : return 1000360840; // Kr84
1056 case 95 : return 1000380880; // Sr88
1057 case 96 : return 1000400900; // Zr90
1058 case 97 : return 1000420980; // Mo98
1059 case 98 : return 1000461060; // Pd106
1060 case 99 : return 1000481140; // Cd114
1061 case 100 : return 1000501200; // Sn120
1062 case 101 : return 1000541320; // Xe132
1063 case 102 : return 1000561380; // Ba138
1064 case 103 : return 1000581400; // Ce140
1065 case 104 : return 1000621520; // Sm152
1066 case 105 : return 1000661640; // Dy164
1067 case 106 : return 1000701740; // Yb174
1068 case 107 : return 1000741840; // W184
1069 case 108 : return 1000781940; // Pt194
1070 case 109 : return 1000791970; // Au197
1071 case 110 : return 1000802020; // Hg202
1072 case 111 : return 1000822080; // Pb208
1073 case 112 : return 1000922380; // U238
1074
1075 // These are "private" Geant3 types that were defined in hdgeant
1076 case 33 : return 223; // omega(782)
1077 case 34 : return 333; // phi(1020)
1078 case 35 : return 331; // etaPrime(958)
1079 case 36 : return 0; // unused
1080 case 37 : return 0; // unused
1081 case 38 : return 0; // unused
1082 case 39 : return 0; // unused
1083 case 40 : return 0; // unused
1084 case 41 : return 0; // unused
1085 case 42 : return 213; // rho(770)+
1086 case 43 : return -213; // rho(770)-
1087 case 44 : return 113; // rho(770)0
1088
1089 case 182 : return 2224; // Delta++
1090 case 183 : return 443; // Jpsi
1091 case 184 : return 441; // Eta_c
1092 case 185 : return 10441; // Chi_c0
1093 case 186 : return 20443; // Chi_c1
1094 case 187 : return 445; // Chi_c2
1095 case 188 : return 100443; // Psi2s
1096 case 189 : return 421; // D0
1097 case 190 : return 411; // D+
1098 case 191 : return 10421; // Dstar0
1099 case 192 : return 10411; // Dstar+
1100 case 193 : return 4022; // Lambda_c+
1101 case 194 : return -421; // anti-D0
1102
1103 case 163 : return 9000111; // a0(980)
1104 case 164 : return 9010221; // f0(980)
1105 case 165 : return 313; // K*(892)0
1106 case 166 : return 323; // K*(892)+
1107 case 167 : return -323; // K*(892)-
1108 case 168 : return -313; // anti-K*(892)0
1109 case 169 : return 20323; // K1(1400)+
1110 case 170 : return -20323; // K1(1400)-
1111 case 171 : return 4122; // b1(1235)+
1112 case 172 : return 3224; // Sigma*(1385)+
1113 case 173 : return 3214; // Sigma*(1385)0
1114 case 174 : return 3114; // Sigma*(1385)-
1115
1116 default :
1117 G4cout(*G4cout_p) << "Warning in GlueXPrimaryGeneratorAction::"
1118 "ConvertGeant3ToPdg - lookup performed on unknown"
1119 " Geant3 particle type " << Geant3type << ","
1120 " returning 0 for the PDG particle code." << G4endlstd::endl;
1121 }
1122 return 0;
1123}
1124
1125// Convert particle types from PDG scheme to Geant3 types
1126
1127int GlueXPrimaryGeneratorAction::ConvertPdgToGeant3(int PDGtype)
1128{
1129 // Invert the table contained in ConvertGeant3ToPdg
1130
1131 switch (PDGtype) {
1132 case 0 : return 50; // optical photon
1133 case 22 : return 1; // photon
1134 case -11 : return 2; // e+
1135 case 11 : return 3; // e-
1136 case -12 : return 4; // anti-e-neutrino
1137 case 12 : return 4; // e-neutrino
1138 case -13 : return 5; // mu+
1139 case 13 : return 6; // mu-
1140 case -14 : return 4; // anti-mu-neutrino
1141 case 14 : return 4; // mu-neutrino
1142 case -16 : return 4; // anti-tau-neutrino
1143 case 16 : return 4; // tau-neutrino
1144 case 111 : return 7; // pi0
1145 case 211 : return 8; // pi+
1146 case -211 : return 9; // pi-
1147 case 130 : return 10; // K long
1148 case 321 : return 11; // K+
1149 case -321 : return 12; // K-
1150 case 2112 : return 13; // n
1151 case 2212 : return 14; // p
1152 case -2212 : return 15; // anti-proton
1153 case 310 : return 16; // K short
1154 case 221 : return 17; // eta
1155 case 3122 : return 18; // Lambda
1156 case 3222 : return 19; // Sigma+
1157 case 3212 : return 20; // Sigma0
1158 case 3112 : return 21; // Sigma-
1159 case 3322 : return 22; // Xi0
1160 case 3312 : return 23; // Xi-
1161 case 3334 : return 24; // Omega-
1162 case -2112 : return 25; // anti-neutron
1163 case -3122 : return 26; // anti-Lambda
1164 case -3222 : return 27; // Sigma-
1165 case -3212 : return 28; // Sigma0
1166 case -3112 : return 29; // Sigma+
1167 case -3322 : return 30; // Xi0
1168 case -3312 : return 31; // Xi+
1169 case -3334 : return 32; // Omega+
1170 case 1000010020 : return 45; // deuteron
1171 case 1000010030 : return 46; // triton
1172 case 1000020040 : return 47; // alpha
1173 case 1000020030 : return 49; // He3
1174 case 1000030060 : return 61; // Li6
1175 case 1000030070 : return 62; // Li7
1176 case 1000040070 : return 63; // Be7
1177 case 1000040090 : return 64; // Be9
1178 case 1000050100 : return 65; // B10
1179 case 1000050110 : return 66; // B11
1180 case 1000060120 : return 67; // C12
1181 case 1000070140 : return 68; // N14
1182 case 1000080160 : return 69; // O16
1183 case 1000090190 : return 70; // F19
1184 case 1000100200 : return 71; // Ne20
1185 case 1000110230 : return 72; // Na23
1186 case 1000120240 : return 73; // Mg24
1187 case 1000130270 : return 74; // Al27
1188 case 1000140280 : return 75; // Si28
1189 case 1000150310 : return 76; // P31
1190 case 1000160320 : return 77; // S32
1191 case 1000170350 : return 78; // Cl35
1192 case 1000180360 : return 79; // Ar36
1193 case 1000190390 : return 80; // K39
1194 case 1000200400 : return 81; // Ca40
1195 case 1000210450 : return 82; // Sc45
1196 case 1000220480 : return 83; // Ti48
1197 case 1000230510 : return 84; // V51
1198 case 1000240520 : return 85; // Cr52
1199 case 1000250550 : return 86; // Mn55
1200 case 1000260560 : return 87; // Fe56
1201 case 1000270590 : return 88; // Co59
1202 case 1000280580 : return 89; // Ni58
1203 case 1000290630 : return 90; // Cu63
1204 case 1000300640 : return 91; // Zn64
1205 case 1000320740 : return 92; // Ge74
1206 case 1000340800 : return 93; // Se80
1207 case 1000360840 : return 94; // Kr84
1208 case 1000380880 : return 95; // Sr88
1209 case 1000400900 : return 96; // Zr90
1210 case 1000420980 : return 97; // Mo98
1211 case 1000461060 : return 98; // Pd106
1212 case 1000481140 : return 99; // Cd114
1213 case 1000501200 : return 100; // Sn120
1214 case 1000541320 : return 101; // Xe132
1215 case 1000561380 : return 102; // Ba138
1216 case 1000581400 : return 103; // Ce140
1217 case 1000621520 : return 104; // Sm152
1218 case 1000661640 : return 105; // Dy164
1219 case 1000701740 : return 106; // Yb174
1220 case 1000741840 : return 107; // W184
1221 case 1000781940 : return 108; // Pt194
1222 case 1000791970 : return 109; // Au197
1223 case 1000802020 : return 110; // Hg202
1224 case 1000822080 : return 111; // Pb208
1225 case 1000922380 : return 112; // U238
1226
1227 // These are "private" Geant3 types that were defined in hdgeant
1228 case 223 : return 33; // omega(782)
1229 case 333 : return 34; // phi(1020)
1230 case 331 : return 35; // etaPrime(958)
1231 case 213 : return 42; // rho(770)+
1232 case -213 : return 43; // rho(770)-
1233 case 113 : return 44; // rho(770)0
1234
1235 case 2224 : return 182; // Delta++
1236 case 443 : return 183; // Jpsi
1237 case 441 : return 184; // Eta_c
1238 case 10441 : return 185; // Chi_c0
1239 case 20443 : return 186; // Chi_c1
1240 case 445 : return 187; // Chi_c2
1241 case 100443 : return 188; // Psi2s
1242 case 421 : return 189; // D0
1243 case 411 : return 190; // D+
1244 case 10421 : return 191; // Dstar0
1245 case 10411 : return 192; // Dstar+
1246 case 4022 : return 193; // Lambda_c+
1247 case -421 : return 194; // anti-D0
1248
1249 case 9000111 : return 163; // a0(980)
1250 case 9010221 : return 164; // f0(980)
1251 case 313 : return 165; // K*(892)0
1252 case 323 : return 166; // K*(892)+
1253 case -323 : return 167; // K*(892)-
1254 case -313 : return 168; // anti-K*(892)0
1255 case 20323 : return 169; // K1(1400)+
1256 case -20323 : return 170; // K1(1400)-
1257 case 4122 : return 171; // b1(1235)+
1258 case 3224 : return 172; // Sigma*(1385)+
1259 case 3214 : return 173; // Sigma*(1385)0
1260 case 3114 : return 174; // Sigma*(1385)-
1261
1262 default :
1263 if (PDGtype < 1000000000) {
1264 G4cout(*G4cout_p) << "Warning in GlueXPrimaryGeneratorAction::"
1265 "ConvertPdgToGeant3 - lookup performed on unknown"
1266 " PDG particle type " << PDGtype << ","
1267 " returning 0 for the Geant3 particle code." << G4endlstd::endl;
1268 }
1269 }
1270 return 0;
1271}
1272
1273double GlueXPrimaryGeneratorAction::GetMassPDG(int PDGtype)
1274{
1275 return fParticleTable->FindParticle(PDGtype)->GetPDGMass();
1276}
1277
1278double GlueXPrimaryGeneratorAction::GetMass(int Geant3Type)
1279{
1280 return GetMassPDG(ConvertGeant3ToPdg(Geant3Type));
1281}
1282
1283G4ParticleDefinition *GlueXPrimaryGeneratorAction::GetParticle(int PDGtype)
1284{
1285 G4ParticleDefinition *p = fParticleTable->FindParticle(PDGtype);
1286 if (p==0) {
1287 if (PDGtype > 1000000000) {
1288 p = fParticleTable->GetIonTable()->GetIon(PDGtype);
1289 }
1290 else {
1291 G4cout(*G4cout_p) << "unknown particle type " << PDGtype
1292 << ", substituting geantino in its place!"
1293 << G4endlstd::endl;
1294 p = fParticleTable->FindParticle("geantino");
1295 }
1296 }
1297 return p;
1298}
1299
1300G4ParticleDefinition *GlueXPrimaryGeneratorAction::GetParticle(const G4String &name)
1301{
1302 G4ParticleDefinition *p = fParticleTable->FindParticle(name);
1303 if (p==0) {
1304 G4cout(*G4cout_p) << "unknown particle type " << name
1305 << ", substituting geantino in its place!"
1306 << G4endlstd::endl;
1307 p = fParticleTable->FindParticle("geantino");
1308 }
1309 return p;
1310}
1311
1312void GlueXPrimaryGeneratorAction::configure_beam_vertex()
1313{
1314 GlueXUserOptions *user_opts = GlueXUserOptions::GetInstance();
1315 if (user_opts == 0) {
1316 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction::configure_beam_vertex"
1317 << " - GlueXUserOptions::GetInstance() returns null, "
1318 << "cannot continue." << G4endlstd::endl;
1319 exit(-1);
1320 }
1321 std::map<int, std::string> vertex_spec;
1322 if (user_opts->Find("VERTEX", vertex_spec)) {
1323 if (vertex_spec.find(1) == vertex_spec.end()) {
1324 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction::configure_beam_vertex"
1325 << " - VERTEX card found but no specification was given."
1326 << G4endlstd::endl;
1327 }
1328 }
1329 else {
1330 fBeamvertex_activated = -1;
1331 return;
1332 }
1333 const char *spec = vertex_spec[1].c_str();
1334 if (sscanf(spec, "beam_spot(ccdb) * %lf", &fBeamvertex.length)) {
1335 int runno = HddmOutput::getRunNo();
1336 extern jana::JApplication *japp;
1337 if (japp == 0) {
1338 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction::configure_beam_vertex"
1339 << " - jana global DApplication object not set, "
1340 << "cannot continue." << G4endlstd::endl;
1341 exit(-1);
1342 }
1343 jana::JCalibration *jcalib = japp->GetJCalibration(runno);
1344 std::map<string, float> beam_spot_params;
1345 jcalib->Get("/PHOTON_BEAM/beam_spot", beam_spot_params);
1346 fBeamvertex.x = beam_spot_params.at("x") * cm;
1347 fBeamvertex.y = beam_spot_params.at("y") * cm;
1348 fBeamvertex.z = beam_spot_params.at("z") * cm;
1349 fBeamvertex.var_xx = beam_spot_params.at("var_xx") * cm*cm;
1350 fBeamvertex.var_xy = beam_spot_params.at("var_xy") * cm*cm;
1351 fBeamvertex.var_yy = beam_spot_params.at("var_yy") * cm*cm;
1352 fBeamvertex.dxdz = beam_spot_params.at("dxdz");
1353 fBeamvertex.dydz = beam_spot_params.at("dydz");
1354 fBeamvertex.length *= cm;
1355 }
1356 else if (sscanf(spec, "beam_spot(%*[-+0-9.,]) * %lf", &fBeamvertex.length)) {
1357 sscanf(spec, "beam_spot(%lf,%lf,%lf,%lf,%lf,%lf,%lf,%lf)",
1358 &fBeamvertex.x, &fBeamvertex.y, &fBeamvertex.z,
1359 &fBeamvertex.var_xx, &fBeamvertex.var_xy,
1360 &fBeamvertex.var_yy, &fBeamvertex.dxdz,
1361 &fBeamvertex.dydz);
1362 fBeamvertex.x *= cm;
1363 fBeamvertex.y *= cm;
1364 fBeamvertex.z *= cm;
1365 fBeamvertex.var_xx *= cm*cm;
1366 fBeamvertex.var_xy *= cm*cm;
1367 fBeamvertex.var_yy *= cm*cm;
1368 fBeamvertex.length *= cm;
1369 }
1370 else {
1371 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction::configure_beam_vertex"
1372 << " - unrecognized VERTEX specification on control.in line"
1373 << G4endlstd::endl << spec << G4endlstd::endl;
1374 exit(-1);
1375 }
1376 double D = fBeamvertex.var_xx * fBeamvertex.var_yy -
1377 fBeamvertex.var_xy * fBeamvertex.var_xy;
1378 double A = (fBeamvertex.var_xx + fBeamvertex.var_yy)/2;
1379 double B = (fBeamvertex.var_xx - fBeamvertex.var_yy)/2;
1380 double C = fBeamvertex.var_xy;
1381 double evalue1 = (A + sqrt(B*B + C*C) + 1e-20) / (D + 1e-40);
1382 double evalue2 = (A - sqrt(B*B + C*C) + 1e-20) / (D + 1e-40);
1383 if (evalue1 < 0 || evalue2 < 0) {
1384 G4cerr(*G4cerr_p) << "Error in GlueXPrimaryGeneratorAction::configure_beam_vertex"
1385 << " - unphysical values given for the beam spot ellipse:"
1386 << " var_xx=" << fBeamvertex.var_xx
1387 << ", var_xy=" << fBeamvertex.var_xy
1388 << ", var_yy=" << fBeamvertex.var_yy
1389 << G4endlstd::endl;
1390 exit(1);
1391 }
1392 double alpha = 0;
1393 if (C == 0)
1394 alpha = (B < 0)? 0 : M_PI3.14159265358979323846/2;
1395 else if (B == 0)
1396 alpha = (C < 0)? -M_PI3.14159265358979323846/4 : M_PI3.14159265358979323846/4;
1397 else
1398 alpha = atan2(C,B) / 2;
1399 fBeamvertex.sigma[0] = sqrt(1 / evalue1);
1400 fBeamvertex.sigma[1] = sqrt(1 / evalue2);
1401 fBeamvertex.alpha = alpha;
1402 fBeamvertex_activated = 1;
1403
1404#ifndef QUIET_CONFIGURE_BEAM_VERTEX
1405 G4cout(*G4cout_p) << "Configured beam vertex parameters: (units of cm)" << G4endlstd::endl
1406 << " center x=" << fBeamvertex.x/cm
1407 << " y=" << fBeamvertex.y/cm
1408 << " z=" << fBeamvertex.z/cm << G4endlstd::endl
1409 << " sigma_x=" << sqrt(fBeamvertex.var_xx)/cm
1410 << " sigma_y=" << sqrt(fBeamvertex.var_yy)/cm
1411 << " sigma_xy=" << fBeamvertex.var_xy/(cm*cm) << G4endlstd::endl
1412 << " dxdz=" << fBeamvertex.dxdz
1413 << " dydz=" << fBeamvertex.dydz
1414 << " dz=" << fBeamvertex.length/cm << G4endlstd::endl;
1415#endif
1416}
1417
1418void GlueXPrimaryGeneratorAction::generate_beam_vertex(double v[3])
1419{
1420}
1421
1422void GlueXPrimaryGeneratorAction::clone_photon_beam_generator()
1423{
1424 fCobremsGeneration = new CobremsGeneration(fBeamEndpointEnergy/GeV,
1425 fBeamPeakEnergy/GeV);
1426 CobremsGeneration *gen0 = (*fInstance.begin())->fCobremsGeneration;
1427 double beamEmin = gen0->getPhotonEnergyMin()*GeV;
1428 double radColDist = gen0->getCollimatorDistance()*m;
1429 double colDiam = gen0->getCollimatorDiameter()*m;
1430 double beamEmit = gen0->getBeamEmittance()*(m*radian);
1431 double radThick = gen0->getTargetThickness()*m;
1432 double spotRMS = gen0->getCollimatorSpotrms()*m;
1433 GlueXPhotonBeamGenerator *gen1 = (*fInstance.begin())->fPhotonBeamGenerator;
1434 double spotX = gen1->getBeamOffset(0);
1435 double spotY = gen1->getBeamOffset(1);
1436 fCobremsGeneration->setPhotonEnergyMin(beamEmin/GeV);
1437 fCobremsGeneration->setCollimatorDistance(radColDist/m);
1438 fCobremsGeneration->setCollimatorDiameter(colDiam/m);
1439 fCobremsGeneration->setBeamEmittance(beamEmit/(m*radian));
1440 fCobremsGeneration->setTargetThickness(radThick/m);
1441 fCobremsGeneration->setCollimatorSpotrms(spotRMS/m);
1442 fPhotonBeamGenerator = new GlueXPhotonBeamGenerator(fCobremsGeneration);
1443 fPhotonBeamGenerator->setBeamOffset(spotX, spotY);
1444}