Bug Summary

File:programs/Simulation/gen_primex_eta_he4/gen_primex_eta_he4.cc
Location:line 219, column 12
Description:Value stored to 's' during its initialization is never read

Annotated Source Code

1/**************************************************************************
2* HallD software *
3* Copyright(C) 2019 GlueX and PrimEX-D Collaborations *
4* *
5* Author: The GlueX and PrimEX-D Collaborations *
6* Contributors: Igal Jaegle *
7* *
8* This software is provided "as is" without any warranty. *
9**************************************************************************/
10
11#include <iostream>
12#include <fstream>
13#include <complex>
14#include <string>
15#include <vector>
16#include <utility>
17#include <map>
18#include <cassert>
19#include <cstdlib>
20
21#include "particleType.h"
22
23#include "UTILITIES/BeamProperties.h"
24
25#include "AMPTOOLS_DATAIO/ROOTDataWriter.h"
26#include "AMPTOOLS_DATAIO/HDDMDataWriter.h"
27
28#include "AMPTOOLS_AMPS/Compton.h"
29
30#include "AMPTOOLS_MCGEN/GammaPToXP.h"
31
32#include "IUAmpTools/AmpToolsInterface.h"
33#include "IUAmpTools/ConfigFileParser.h"
34
35#include "TF1.h"
36#include "TH1F.h"
37#include "TH2F.h"
38#include "TFile.h"
39#include "TLorentzVector.h"
40#include "TLorentzRotation.h"
41#include <TGenPhaseSpace.h>
42#include "TRandom3.h"
43#include "TSystem.h"
44
45#include "HddmOut.h"
46#include "UTILITIES/MyReadConfig.h"
47
48using std::complex;
49using namespace std;
50
51#define eta_TYPE17 17
52#define pi0_TYPE7 7
53#define gamma_TYPE1 1
54#define Helium_TYPE47 47
55
56int main( int argc, char* argv[] ){
57
58 string beamconfigfile("");
59 TString genconfigfile("");// = "gen_config.dat";
60 string outname("");
61 string hddmname("");
62
63 ifstream in_coherent;
64 ifstream in_incoherent;
65
66 double beamLowE = 3.0;
67 double beamHighE = 12.0;
68 const double M_He4 = 3.727379378;
69 const double M_eta = 0.54730;
70 const double M_pi0 = 1.3957018e-01;
71 const double M_gamma = 0.0;
72 TLorentzVector Target_4Vec(0, 0, 0, M_He4);
73
74 int runNum = 9001;
75 int seed = 0;
76
77 int nEvents = 10000;
78
79 //double dummy;
80
81 //parse command line:
82 for (int i = 1; i < argc; i++) {
83
84 string arg(argv[i]);
85
86
87 if (arg == "-c") {
88 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
89 else beamconfigfile = argv[++i];
90 }
91 if (arg == "-e") {
92 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
93 else genconfigfile = argv[++i];
94 }
95 if (arg == "-o") {
96 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
97 else outname = argv[++i];
98 }
99 if (arg == "-hd") {
100 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
101 else hddmname = argv[++i];
102 }
103 if (arg == "-n"){
104 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
105 else nEvents = atoi( argv[++i] );
106 }
107 if (arg == "-a") {
108 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
109 else beamLowE = atof( argv[++i] );
110 }
111 if (arg == "-b") {
112 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
113 else beamHighE = atof( argv[++i] );
114 }
115 if (arg == "-r") {
116 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
117 else runNum = atoi( argv[++i] );
118 }
119 if (arg == "-s") {
120 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
121 else seed = atoi( argv[++i] );
122 }
123 if (arg == "-h") {
124 cout << endl << " Usage for: " << argv[0] << endl << endl;
125 cout << "\t -c <file>\t Beam config file" << endl;
126 cout << "\t -e <file>\t Generator config file" << endl;
127 cout << "\t -o <name>\t ASCII file output name" << endl;
128 cout << "\t -hd <name>\t HDDM file output name [optional]" << endl;
129 cout << "\t -n <value>\t Number of events to generate [optional]" << endl;
130 cout << "\t -a <value>\t Minimum photon energy to simulate events [optional]" << endl;
131 cout << "\t -b <value>\t Maximum photon energy to simulate events [optional]" << endl;
132 cout << "\t -r <value>\t Run number assigned to generated events [optional]" << endl;
133 cout << "\t -s <value>\t Random number seed initialization [optional]" << endl;
134 exit(1);
135 }
136 }
137
138 if (outname.size() == 0 && hddmname == "") {
139 cout << "No output specificed: run gen_compton_simple -h for help" << endl;
140 exit(1);
141 }
142
143 if (genconfigfile == "") {
144 cout << "No generator configuration file: run gen_primex_eta_he4 -h for help " << endl;
145 exit(1);
146 }
147 // random number initialization (set to 0 by default)
148 gRandom->SetSeed(seed);
149
150 // initialize HDDM output
151 HddmOut *hddmWriter = nullptr;
152 if (hddmname != "")
153 hddmWriter = new HddmOut(hddmname.c_str());
154
155 // initialize ASCII output
156 ofstream *asciiWriter = nullptr;
157 if (outname != "")
158 asciiWriter = new ofstream(outname.c_str());
159
160 // Assume a beam energy spectrum of 1/E(gamma)
161 TF1 ebeam_spectrum("beam_spectrum","1/x",beamLowE,beamHighE);
162
163 // Get beam properties from configuration file
164 TH1D * cobrem_vs_E = 0;
165 if (beamconfigfile != "") {
166 BeamProperties beamProp( beamconfigfile );
167 cobrem_vs_E = (TH1D*)beamProp.GetFlux();
168 }
169
170 // Get generator config file
171 MyReadConfig * ReadFile = new MyReadConfig();
172 ReadFile->ReadConfigFile(genconfigfile);
173 TString m_rfile = ReadFile->GetConfigName("rfile");
174 TString m_histo = ReadFile->GetConfigName("histo");
175 TString m_decay = ReadFile->GetConfigName("decay");
176 cout << "rfile " << m_rfile << endl;
177 cout << "histo " << m_histo << endl;
178 cout << "decay " << m_decay << endl;
179 // Load eta-meson differential cross-section based on Ilya Larin's calculation, see the *.F program in this directory
180 TFile * ifile = new TFile(m_rfile);
181 TH2F * h_dxs = (TH2F *) ifile->Get(m_histo);
182
183 double M_meson = 0;
184 if (m_decay.Contains("eta"))
185 M_meson = M_eta;
186 else if (m_decay.Contains("pi0"))
187 M_meson = M_pi0;
188
189 // Create decayGen
190 TGenPhaseSpace decayGen;
191
192 TFile* diagOut = new TFile( "gen_primex_eta_he4_diagnostic.root", "recreate" );
193 TH2F* h_Tkin_eta_vs_egam = new TH2F("Tkin_eta_vs_egam", ";E_{#gamma} [GeV];T^{kin}_{#eta} [GeV];Count [a.u.]", 1000, 0.0, 12.0, 1000, 0.0, 12.0);
194 TH2F* h_Tkin_photon_vs_egam = new TH2F("Tkin_photon_vs_egam", ";E_{#gamma} [GeV];T^{kin}_{#gamma} [GeV];Count [a.u.]", 1000, 0.0, 12.0, 1000, 0.0, 12.0);
195 TH2F* h_Tkin_recoilA_vs_egam = new TH2F("Tkin_recoilA_vs_egam", ";E_{#gamma} [GeV];T^{kin}_{A} [GeV];Count [a.u.]", 1000, 0.0, 12.0, 1000, 0.0, 1.0);
196 TH2F* h_theta_eta_vs_egam = new TH2F("theta_eta_vs_egam", ";E_{#gamma} [GeV];#theta_{#eta} [GeV];Count [a.u.]", 1000, 0.0, 12.0, 1000, 0., 10.);
197 TH2F* h_theta_photon_vs_egam = new TH2F("theta_photon_vs_egam", ";E_{#gamma} [GeV];#theta_{#gamma} [GeV];Count [a.u.]", 1000, 0.0, 12.0, 1000, 0., 180.);
198 TH2F* h_theta_recoilA_vs_egam = new TH2F("theta_recoilA_vs_egam", ";E_{#gamma} [GeV];#theta_{A} [GeV];Count [a.u.]", 1000, 0.0, 12.0, 1000, 0., 180.);
199
200 for (int i = 0; i < nEvents; ++i) {
201 if (i%1000 == 1)
202 cout << "event " << i <<endl;
203
204 // get beam energy
205 double ebeam = 0;
206 if (beamconfigfile == "" || cobrem_vs_E == 0)
207 ebeam = ebeam_spectrum.GetRandom();
208 else if (beamconfigfile != "")
209 ebeam = cobrem_vs_E->GetRandom();
210
211 // Incident photon-beam 4Vec
212 TLorentzVector InGamma_4Vec(0, 0, ebeam, ebeam);
213
214 // Initial state 4Vec
215 TLorentzVector IS_4Vec = InGamma_4Vec + Target_4Vec;
216
217 // Mass in the centre-of-mass frame
218 double sqrt_s = IS_4Vec.M();
219 double s = pow(sqrt_s, 2);
Value stored to 's' during its initialization is never read
220
221 // Histo. creation that will store the calculated diff. xs. vs. LAB polar angle
222 int ebeam_bin = h_dxs->GetXaxis()->FindBin(ebeam);
223 TH1F * h_ThetaLAB = (TH1F *) h_dxs->ProjectionY("h_ThetaLAB", ebeam_bin, ebeam_bin);
224 if (h_ThetaLAB->GetEntries() == 0) continue;
225
226 /*
227 // XS initialization
228 double xs_tot = 0, xs_Primakoff = 0, xs_interference = 0, xs_strong_coherent = 0, xs_incoherent = 0;
229
230 // read or calculate differential cross-section
231 // open coh. diff. xs
232 in_coherent.open(TString::Format("ds_eta_he4_%0.3f.dat", ebeam));
233
234 // open incoh. diff. xs
235 in_incoherent.open(TString::Format("inc_eta_he4_%0.3f.dat", ebeam));
236
237 // if already calculated, read
238 int j = 0;
239 if (in_coherent.good() && in_incoherent.good()) {
240
241 // Read Primakoff, interference, and strong coherent diff. xs terms
242 j = 0;
243 while (in_coherent.good()) {
244 xs_tot = 0; xs_Primakoff = 0; xs_interference = 0; xs_strong_coherent = 0; xs_incoherent = 0;
245 in_coherent >> dummy >> dummy >> xs_Primakoff >> xs_interference >> xs_strong_coherent >> xs_incoherent;
246 xs_tot = xs_Primakoff + xs_interference + xs_strong_coherent + xs_incoherent;
247 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_tot);
248 j ++;
249 }
250 in_coherent.close();
251
252 // Read incoh. term
253 j = 0;
254 while (in_incoherent.good()) {
255 xs_incoherent = 0;
256 in_incoherent >> dummy >> xs_incoherent >> dummy >> dummy;
257 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_incoherent);
258 j ++;
259 }
260 in_incoherent.close();
261
262 } else { // If not calculated, calculate and then read
263
264 // Close files if not already closed
265 in_coherent.close();
266 in_incoherent.close();
267
268 // Calculate and produce diff. xs dat files for 100 bin from 0 to 10 degree
269 // Calculate Coulomb term of the coherent diff. xs
270 system(TString::Format("ff_coulomb %0.03f ff_coulom_%0.3f.dat", ebeam, ebeam));
271
272 // Calculate strong term of the coherent diff. xs
273 system(TString::Format("ff_strong %0.03f ff_strong_%0.3f.dat", ebeam, ebeam));
274
275 // Calculate Primakoff, interference, and strong coherent diff. xs
276 system(TString::Format("ds_eta_he4 %0.03f ff_coulom_%0.3f.dat ff_strong_%0.3f.dat ds_eta_he4_%0.3f.dat", ebeam, ebeam, ebeam, ebeam));
277
278 // Calculate incoherent
279 system(TString::Format("inc_eta_he4 %0.03f inc_eta_he4_%0.3f.dat", ebeam, ebeam));
280
281 // Read Primakoff, interference, and strong coherent diff. xs terms
282 in_coherent.open(TString::Format("ds_eta_he4_%0.3f.dat", ebeam));
283 j = 0;
284 while (in_coherent.good()) {
285 xs_tot = 0; xs_Primakoff = 0; xs_interference = 0; xs_strong_coherent = 0; xs_incoherent = 0;
286 in_coherent >> dummy >> dummy >> xs_Primakoff >> xs_interference >> xs_strong_coherent >> xs_incoherent;
287 xs_tot = xs_Primakoff + xs_interference + xs_strong_coherent + xs_incoherent;
288 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_tot);
289 j ++;
290 }
291 in_coherent.close();
292
293 // Read incoh. term
294 j = 0;
295 while (in_incoherent.good()) {
296 xs_incoherent = 0;
297 in_incoherent >> dummy >> xs_incoherent >> dummy >> dummy;
298 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_incoherent);
299 j ++;
300 }
301 in_incoherent.close();
302 }
303 */
304 // Generate eta-meson theta in LAB
305 double ThetaLAB = h_ThetaLAB->GetRandom();
306 ThetaLAB *= TMath::DegToRad();
307
308 // Generate eta-meson phi in LAB
309 double PhiLAB = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
310
311 // Calculate eta-meson energy in COM
312 //double E_COM_eta = (s - pow(M_He4, 2) + pow(M_eta, 2)) / (2.0 * sqrt_s);
313
314 // Calculate eta-meson momentum in COM
315 //double P_COM_eta = sqrt(pow(E_COM_eta, 2) - pow(M_eta, 2));
316
317 // Calculate eta-meson energy in LAB
318 double E_LAB_eta = 0;
319 double p0 = IS_4Vec.P();
320 double E0 = IS_4Vec.E();
321 double m2 = M_meson;
322 double m3 = M_He4;
323 double a = pow(2.0 * p0 * cos(ThetaLAB), 2) - pow(2.0 * E0, 2);
324 double b = 4.0 * pow(E0, 3) - pow(2.0 * p0, 2) * E0 + pow(2.0 * m2, 2) * E0 - pow(2.0 * m3, 2) * E0;
325 double c = 2.0 * pow(p0 * E0, 2) - 2.0 * pow(m2 * E0, 2) + 2.0 * pow(m2 * p0, 2) + 2.0 * pow(m3 * E0, 2) - 2.0 * pow(m3 * p0, 2) +
326 2.0 * pow(m3 * m2, 2) - pow(E0, 4) - pow(p0, 4) - pow(m2, 4) - pow(m3, 4) - pow(2.0 * m2 * p0 * cos(ThetaLAB), 2);
327 double E_LAB_eta0 = (-b + sqrt(pow(b, 2) - 4.0 * a * c)) / (2.0 * a);
328 double E_LAB_eta1 = (-b - sqrt(pow(b, 2) - 4.0 * a * c)) / (2.0 * a);
329 if (E_LAB_eta0 < 0) E_LAB_eta = E_LAB_eta1;
330 if (E_LAB_eta1 < 0) E_LAB_eta = E_LAB_eta0;
331 if (E_LAB_eta0 > E_LAB_eta1) E_LAB_eta = E_LAB_eta0;
332 if (E_LAB_eta1 > E_LAB_eta0) E_LAB_eta = E_LAB_eta1;
333
334 // Calculate eta-meson momentun im LAB
335 double P_LAB_eta = sqrt(pow(E_LAB_eta, 2) - pow(M_eta, 2));
336
337 // Calculate the momentum for each direction
338 double Px_LAB_eta = P_LAB_eta * sin(ThetaLAB) * cos(PhiLAB);
339 double Py_LAB_eta = P_LAB_eta * sin(ThetaLAB) * sin(PhiLAB);
340 double Pz_LAB_eta = P_LAB_eta * cos(ThetaLAB);
341
342 // Store the results in TLorentzVector for the eta-meson
343 TLorentzVector eta_LAB_4Vec(Px_LAB_eta, Py_LAB_eta, Pz_LAB_eta, E_LAB_eta);
344 //TLorentzVector eta_COM_4Vec = eta_LAB_4Vec;
345 //eta_COM_4Vec.Boost(-IS_4Vec.BoostVector());
346
347 // Make the eta-meson decay into two photons
348
349 TLorentzVector photon_4Vec[6];
350 TLorentzVector pi0_4Vec[3];
351 int ng_max = 0;
352 if (m_decay == "pi0->2g") {
353 ng_max = 2;
354 double masses[] = {M_gamma, M_gamma};
355 if (decayGen.SetDecay(eta_LAB_4Vec, 2, masses)) {
356 decayGen.Generate();
357 photon_4Vec[0] = * decayGen.GetDecay(0);
358 photon_4Vec[1] = * decayGen.GetDecay(1);
359 }
360 } else if (m_decay == "eta->2g") {
361 ng_max = 2;
362 double masses[] = {M_gamma, M_gamma};
363 if (decayGen.SetDecay(eta_LAB_4Vec, 2, masses)) {
364 decayGen.Generate();
365 photon_4Vec[0] = * decayGen.GetDecay(0);
366 photon_4Vec[1] = * decayGen.GetDecay(1);
367 }
368 } else if (m_decay == "eta->6g") {
369 ng_max = 6;
370 double masses[] = {M_pi0, M_pi0, M_pi0};
371 if (decayGen.SetDecay(eta_LAB_4Vec, 3, masses)) {
372 decayGen.Generate();
373 pi0_4Vec[0] = * decayGen.GetDecay(0);
374 pi0_4Vec[1] = * decayGen.GetDecay(1);
375 pi0_4Vec[2] = * decayGen.GetDecay(2);
376 }
377 for (int j = 0; j < 3; j ++) {
378 double mass[] = {M_gamma, M_gamma};
379 if (decayGen.SetDecay(pi0_4Vec[j], 2, mass)) {
380 decayGen.Generate();
381 photon_4Vec[0 + 2 * j] = * decayGen.GetDecay(0);
382 photon_4Vec[1 + 2 * j] = * decayGen.GetDecay(1);
383 }
384 }
385 }
386
387
388 // Deduce by energy and mass conservation the recoil nucleus 4Vec
389 TLorentzVector He4_LAB_4Vec = IS_4Vec - eta_LAB_4Vec;
390
391 h_Tkin_recoilA_vs_egam->Fill(ebeam, He4_LAB_4Vec.E() - He4_LAB_4Vec.M());
392 h_theta_recoilA_vs_egam->Fill(ebeam, He4_LAB_4Vec.Theta() * TMath::RadToDeg());
393 h_Tkin_eta_vs_egam->Fill(ebeam, eta_LAB_4Vec.E() - eta_LAB_4Vec.M());
394 h_theta_eta_vs_egam->Fill(ebeam, eta_LAB_4Vec.Theta() * TMath::RadToDeg());
395 for (int j = 0; j < ng_max; j ++) {
396 h_Tkin_photon_vs_egam->Fill(ebeam, photon_4Vec[j].E());
397 h_theta_photon_vs_egam->Fill(ebeam, photon_4Vec[j].Theta() * TMath::RadToDeg());
398 }
399
400 if (hddmWriter) {
401 // ======= HDDM output =========
402 tmpEvt_t tmpEvt;
403 tmpEvt.beam = InGamma_4Vec;
404 tmpEvt.target = Target_4Vec;
405 if (m_decay == "pi0->2g") {
406 tmpEvt.q1 = photon_4Vec[0];
407 tmpEvt.q2 = photon_4Vec[1];
408 tmpEvt.q3 = He4_LAB_4Vec;
409 tmpEvt.nGen = 3;
410 } else if (m_decay == "eta->2g") {
411 tmpEvt.q1 = photon_4Vec[0];
412 tmpEvt.q2 = photon_4Vec[1];
413 tmpEvt.q3 = He4_LAB_4Vec;
414 tmpEvt.nGen = 3;
415 } else if (m_decay == "eta->6g") {
416 tmpEvt.q1 = photon_4Vec[0];
417 tmpEvt.q2 = photon_4Vec[1];
418 tmpEvt.q3 = photon_4Vec[2];
419 tmpEvt.q4 = photon_4Vec[3];
420 tmpEvt.q5 = photon_4Vec[4];
421 tmpEvt.q6 = photon_4Vec[5];
422 tmpEvt.q7 = He4_LAB_4Vec;
423 tmpEvt.nGen = 7;
424 } else if (ng_max == 0) {
425 tmpEvt.q1 = eta_LAB_4Vec;
426 tmpEvt.q2 = He4_LAB_4Vec;
427 tmpEvt.nGen = 2;
428 }
429 tmpEvt.weight = 1.;
430 hddmWriter->write(tmpEvt,runNum,i);
431 }
432 if (asciiWriter) {
433 // ======= ASCII output =========
434 (*asciiWriter)<<runNum<<" "<<i<<" 3"<<endl;
435 // photons from the eta
436 (*asciiWriter)<<"0 "<<gamma_TYPE1<<" "<<M_gamma<<endl;
437 (*asciiWriter)<<" "<<0<<" "<<photon_4Vec[0].Px()<<" "<<photon_4Vec[0].Py()<<" "<<photon_4Vec[0].Pz()<<" "<<photon_4Vec[0].E()<<endl;
438 (*asciiWriter)<<"1 "<<gamma_TYPE1<<" "<<M_gamma<<endl;
439 (*asciiWriter)<<" "<<0<<" "<<photon_4Vec[1].Px()<<" "<<photon_4Vec[1].Py()<<" "<<photon_4Vec[1].Pz()<<" "<<photon_4Vec[1].E()<<endl;
440 // Nucleus recoil
441 (*asciiWriter)<<"2 "<<Helium_TYPE47<<" "<<M_He4<<endl;
442 (*asciiWriter)<<" "<<1<<" "<<He4_LAB_4Vec.Px()<<" "<<He4_LAB_4Vec.Py()<<" "<<He4_LAB_4Vec.Pz()<<" "<<He4_LAB_4Vec.E()<<endl;
443 }
444
445 // deletion
446 delete h_ThetaLAB;
447 }
448 h_Tkin_eta_vs_egam->Write();
449 h_Tkin_photon_vs_egam->Write();
450 h_Tkin_recoilA_vs_egam->Write();
451 h_theta_eta_vs_egam->Write();
452 h_theta_photon_vs_egam->Write();
453 h_theta_recoilA_vs_egam->Write();
454 //h_dxs->Write();
455 //cobrem_vs_Erange->Write();
456 diagOut->Close();
457
458 if (hddmWriter) delete hddmWriter;
459 if (asciiWriter) delete asciiWriter;
460
461 return 0;
462}
463
464