Bug Summary

File:programs/Simulation/gen_primex_eta_he4/gen_primex_eta_he4.cc
Location:line 274, 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 etaprime_TYPE35 35
53#define pi0_TYPE7 7
54#define gamma_TYPE1 1
55#define Helium_TYPE47 47
56#define Be9_TYPE64 64
57#define Proton_TYPE14 14
58#define Neutron_TYPE13 13
59
60int main( int argc, char* argv[] ){
61
62 string beamconfigfile("");
63 TString genconfigfile("");// = "gen_config.dat";
64 string outname("");
65 string hddmname("");
66
67 ifstream in_coherent;
68 ifstream in_incoherent;
69
70 double beamLowE = 3.0;
71 //double beamHighE = 12.0;
72 double beamHighE = 22.011;
73 const double M_He4 = 3.727379378;
74 const double M_Be9 = 8.39479;
75 const double M_p = 0.93827208816;
76 const double M_n = 0.93956542052;
77 const double M_eta = 0.54730;
78 const double M_etapr = 0.95778;
79 const double M_pi0 = 0.13497685;
80 const double M_gamma = 0.0;
81
82 int runNum = 9001;
83 int seed = 0;
84
85 int nEvents = 10000;
86
87 //int bin_egam = 9400;
88 //double egam_min = 2.5995;
89 //double egam_max = 12.0005;
90 int bin_egam = 1000;
91 double egam_min = 12.011;
92 double egam_max = 22.011;
93 int bin_theta = 650;
94 double theta_min = 0.0;
95 double theta_max = 6.5;
96
97 //double dummy;
98
99 //parse command line:
100 for (int i = 1; i < argc; i++) {
101
102 string arg(argv[i]);
103
104
105 if (arg == "-c") {
106 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
107 else beamconfigfile = argv[++i];
108 }
109 if (arg == "-e") {
110 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
111 else genconfigfile = argv[++i];
112 }
113 if (arg == "-o") {
114 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
115 else outname = argv[++i];
116 }
117 if (arg == "-hd") {
118 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
119 else hddmname = argv[++i];
120 }
121 if (arg == "-n"){
122 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
123 else nEvents = atoi( argv[++i] );
124 }
125 if (arg == "-a") {
126 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
127 else beamLowE = atof( argv[++i] );
128 }
129 if (arg == "-b") {
130 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
131 else beamHighE = atof( argv[++i] );
132 }
133 if (arg == "-r") {
134 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
135 else runNum = atoi( argv[++i] );
136 }
137 if (arg == "-s") {
138 if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h";
139 else seed = atoi( argv[++i] );
140 }
141 if (arg == "-h") {
142 cout << endl << " Usage for: " << argv[0] << endl << endl;
143 cout << "\t -c <file>\t Beam config file" << endl;
144 cout << "\t -e <file>\t Generator config file" << endl;
145 cout << "\t -o <name>\t ASCII file output name" << endl;
146 cout << "\t -hd <name>\t HDDM file output name [optional]" << endl;
147 cout << "\t -n <value>\t Number of events to generate [optional]" << endl;
148 cout << "\t -a <value>\t Minimum photon energy to simulate events [optional]" << endl;
149 cout << "\t -b <value>\t Maximum photon energy to simulate events [optional]" << endl;
150 cout << "\t -r <value>\t Run number assigned to generated events [optional]" << endl;
151 cout << "\t -s <value>\t Random number seed initialization [optional]" << endl;
152 exit(1);
153 }
154 }
155
156 if (outname.size() == 0 && hddmname == "") {
157 cout << "No output specificed: run gen_compton_simple -h for help" << endl;
158 exit(1);
159 }
160
161 if (genconfigfile == "") {
162 cout << "No generator configuration file: run gen_primex_eta_he4 -h for help " << endl;
163 exit(1);
164 }
165 // random number initialization (set to 0 by default)
166 gRandom->SetSeed(seed);
167
168 // initialize HDDM output
169 HddmOut *hddmWriter = nullptr;
170 if (hddmname != "")
171 hddmWriter = new HddmOut(hddmname.c_str());
172
173 // initialize ASCII output
174 ofstream *asciiWriter = nullptr;
175 if (outname != "")
176 asciiWriter = new ofstream(outname.c_str());
177
178 // Assume a beam energy spectrum of 1/E(gamma)
179 TF1 ebeam_spectrum("beam_spectrum","1/x",beamLowE,beamHighE);
180
181 // Get beam properties from configuration file
182 TH1D * cobrem_vs_E = 0;
183 if (beamconfigfile != "") {
184 BeamProperties beamProp( beamconfigfile );
185 cobrem_vs_E = (TH1D*)beamProp.GetFlux();
186 }
187
188 // Get generator config file
189 MyReadConfig * ReadFile = new MyReadConfig();
190 ReadFile->ReadConfigFile(genconfigfile);
191 TString m_rfile = ReadFile->GetConfigName("rfile");
192 TString m_histo = ReadFile->GetConfigName("histo");
193 TString m_decay = ReadFile->GetConfigName("decay");
194 TString m_target = ReadFile->GetConfigName("target");
195 TString m_Fermi_file = ReadFile->GetConfigName("Fermi_file");
196 Double_t * m_binning = ReadFile->GetConfig6Par("binning");
197 /*if (m_decay != 0) {
198 bin_egam = (int) m_binning[0];
199 egam_min = m_binning[1];
200 egam_max = m_binning[2];
201 bin_theta = (int) m_binning[3];
202 theta_min = m_binning[4];
203 theta_max = m_binning[5];
204 }*/
205 cout << "rfile " << m_rfile << endl;
206 cout << "histo " << m_histo << endl;
207 cout << "decay " << m_decay << endl;
208 cout << "target " << m_target << endl;
209 //cout << "Fermi_file " << m_Fermi_file << endl;
210
211 double M_target = M_He4;
212 if (m_target == "He4" || m_target == "Helium") M_target = M_He4;
213 if (m_target == "Be9" || m_target == "Beryllium-9") M_target = M_Be9;
214 if (m_target == "Proton") M_target = M_p;
215 if (m_target == "Neutron") M_target = M_n;
216 TLorentzVector Target_4Vec(0, 0, 0, M_target);
217
218 // Load eta-meson differential cross-section based on Ilya Larin's calculation, see the *.F program in this directory
219 TFile * ifile = new TFile(m_rfile);
220 TH2F * h_dxs = (TH2F *) ifile->Get(m_histo);
221 bin_egam = h_dxs->GetNbinsX();
222 egam_min = h_dxs->GetXaxis()->GetXmin();
223 egam_max = h_dxs->GetXaxis()->GetXmax();
224 bin_theta = h_dxs->GetNbinsY();
225 theta_min = h_dxs->GetYaxis()->GetXmin();
226 theta_max = h_dxs->GetYaxis()->GetXmax();
227
228 double M_meson = 0;
229 if (m_decay == "eta")
230 M_meson = M_eta;
231 else if (m_decay == "pi0")
232 M_meson = M_pi0;
233 else if (m_decay == "eta'")
234 M_meson = M_etapr;
235 else if (m_decay == "eta->2g")
236 M_meson = M_eta;
237 else if (m_decay == "eta->6g")
238 M_meson = M_eta;
239 else if (m_decay == "pi0->2g")
240 M_meson = M_pi0;
241 // Create decayGen
242 TGenPhaseSpace decayGen;
243 TGenPhaseSpace decayGenTMP1;
244 TGenPhaseSpace decayGenTMP2;
245 TGenPhaseSpace decayGenTMP3;
246
247 TFile* diagOut = new TFile( "gen_primex_eta_he4_diagnostic.root", "recreate" );
248 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);
249 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);
250 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);
251 TH2F* h_theta_eta_vs_egam = new TH2F("theta_eta_vs_egam", ";E_{#gamma} [GeV];#theta_{#eta} [GeV];Count [a.u.]", bin_egam, egam_min, egam_max, bin_theta, theta_min, theta_max);
252 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.);
253 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.);
254
255 for (int i = 0; i < nEvents; ++i) {
256 if (i%1000 == 1)
257 cout << "event " << i <<endl;
258
259 // get beam energy
260 double ebeam = 0;
261 if (beamconfigfile == "" || cobrem_vs_E == 0)
262 ebeam = ebeam_spectrum.GetRandom();
263 else if (beamconfigfile != "")
264 ebeam = cobrem_vs_E->GetRandom();
265
266 // Incident photon-beam 4Vec
267 TLorentzVector InGamma_4Vec(0, 0, ebeam, ebeam);
268
269 // Initial state 4Vec
270 TLorentzVector IS_4Vec = InGamma_4Vec + Target_4Vec;
271
272 // Mass in the centre-of-mass frame
273 double sqrt_s = IS_4Vec.M();
274 double s = pow(sqrt_s, 2);
Value stored to 's' during its initialization is never read
275
276 // Histo. creation that will store the calculated diff. xs. vs. LAB polar angle
277 int ebeam_bin = h_dxs->GetXaxis()->FindBin(ebeam);
278 TH1F * h_ThetaLAB = (TH1F *) h_dxs->ProjectionY("h_ThetaLAB", ebeam_bin, ebeam_bin);
279 if (h_ThetaLAB->GetEntries() == 0) continue;
280
281 /*
282 // XS initialization
283 double xs_tot = 0, xs_Primakoff = 0, xs_interference = 0, xs_strong_coherent = 0, xs_incoherent = 0;
284
285 // read or calculate differential cross-section
286 // open coh. diff. xs
287 in_coherent.open(TString::Format("ds_eta_he4_%0.3f.dat", ebeam));
288
289 // open incoh. diff. xs
290 in_incoherent.open(TString::Format("inc_eta_he4_%0.3f.dat", ebeam));
291
292 // if already calculated, read
293 int j = 0;
294 if (in_coherent.good() && in_incoherent.good()) {
295
296 // Read Primakoff, interference, and strong coherent diff. xs terms
297 j = 0;
298 while (in_coherent.good()) {
299 xs_tot = 0; xs_Primakoff = 0; xs_interference = 0; xs_strong_coherent = 0; xs_incoherent = 0;
300 in_coherent >> dummy >> dummy >> xs_Primakoff >> xs_interference >> xs_strong_coherent >> xs_incoherent;
301 xs_tot = xs_Primakoff + xs_interference + xs_strong_coherent + xs_incoherent;
302 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_tot);
303 j ++;
304 }
305 in_coherent.close();
306
307 // Read incoh. term
308 j = 0;
309 while (in_incoherent.good()) {
310 xs_incoherent = 0;
311 in_incoherent >> dummy >> xs_incoherent >> dummy >> dummy;
312 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_incoherent);
313 j ++;
314 }
315 in_incoherent.close();
316
317 } else { // If not calculated, calculate and then read
318
319 // Close files if not already closed
320 in_coherent.close();
321 in_incoherent.close();
322
323 // Calculate and produce diff. xs dat files for 100 bin from 0 to 10 degree
324 // Calculate Coulomb term of the coherent diff. xs
325 system(TString::Format("ff_coulomb %0.03f ff_coulom_%0.3f.dat", ebeam, ebeam));
326
327 // Calculate strong term of the coherent diff. xs
328 system(TString::Format("ff_strong %0.03f ff_strong_%0.3f.dat", ebeam, ebeam));
329
330 // Calculate Primakoff, interference, and strong coherent diff. xs
331 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));
332
333 // Calculate incoherent
334 system(TString::Format("inc_eta_he4 %0.03f inc_eta_he4_%0.3f.dat", ebeam, ebeam));
335
336 // Read Primakoff, interference, and strong coherent diff. xs terms
337 in_coherent.open(TString::Format("ds_eta_he4_%0.3f.dat", ebeam));
338 j = 0;
339 while (in_coherent.good()) {
340 xs_tot = 0; xs_Primakoff = 0; xs_interference = 0; xs_strong_coherent = 0; xs_incoherent = 0;
341 in_coherent >> dummy >> dummy >> xs_Primakoff >> xs_interference >> xs_strong_coherent >> xs_incoherent;
342 xs_tot = xs_Primakoff + xs_interference + xs_strong_coherent + xs_incoherent;
343 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_tot);
344 j ++;
345 }
346 in_coherent.close();
347
348 // Read incoh. term
349 j = 0;
350 while (in_incoherent.good()) {
351 xs_incoherent = 0;
352 in_incoherent >> dummy >> xs_incoherent >> dummy >> dummy;
353 h_ThetaLAB->Fill(h_ThetaLAB->GetBinCenter(j + 1), xs_incoherent);
354 j ++;
355 }
356 in_incoherent.close();
357 }
358 */
359 // Generate eta-meson theta in LAB
360 double ThetaLAB = h_ThetaLAB->GetRandom();
361 ThetaLAB *= TMath::DegToRad();
362
363 // Generate eta-meson phi in LAB
364 double PhiLAB = gRandom->Uniform(-TMath::Pi(), TMath::Pi());
365
366 // Calculate eta-meson energy in COM
367 //double E_COM_eta = (s - pow(M_He4, 2) + pow(M_eta, 2)) / (2.0 * sqrt_s);
368
369 // Calculate eta-meson momentum in COM
370 //double P_COM_eta = sqrt(pow(E_COM_eta, 2) - pow(M_eta, 2));
371
372 // Calculate eta-meson energy in LAB
373 double E_LAB_eta = 0;
374 double p0 = IS_4Vec.P();
375 double E0 = IS_4Vec.E();
376 double m2 = M_meson;
377 double m3 = M_target;
378 double a = pow(2.0 * p0 * cos(ThetaLAB), 2) - pow(2.0 * E0, 2);
379 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;
380 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) +
381 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);
382 double E_LAB_eta0 = (-b + sqrt(pow(b, 2) - 4.0 * a * c)) / (2.0 * a);
383 double E_LAB_eta1 = (-b - sqrt(pow(b, 2) - 4.0 * a * c)) / (2.0 * a);
384 if (E_LAB_eta0 < 0) E_LAB_eta = E_LAB_eta1;
385 if (E_LAB_eta1 < 0) E_LAB_eta = E_LAB_eta0;
386 if (E_LAB_eta0 > E_LAB_eta1) E_LAB_eta = E_LAB_eta0;
387 if (E_LAB_eta1 > E_LAB_eta0) E_LAB_eta = E_LAB_eta1;
388
389 // Calculate eta-meson momentun im LAB
390 double P_LAB_eta = sqrt(pow(E_LAB_eta, 2) - pow(M_meson, 2));
391
392 // Calculate the momentum for each direction
393 double Px_LAB_eta = P_LAB_eta * sin(ThetaLAB) * cos(PhiLAB);
394 double Py_LAB_eta = P_LAB_eta * sin(ThetaLAB) * sin(PhiLAB);
395 double Pz_LAB_eta = P_LAB_eta * cos(ThetaLAB);
396
397 // Store the results in TLorentzVector for the eta-meson
398 TLorentzVector eta_LAB_4Vec(Px_LAB_eta, Py_LAB_eta, Pz_LAB_eta, E_LAB_eta);
399 //TLorentzVector eta_COM_4Vec = eta_LAB_4Vec;
400 //eta_COM_4Vec.Boost(-IS_4Vec.BoostVector());
401
402 // Make the eta-meson decay into two photons
403
404 TLorentzVector photon_4Vec[6];
405 TLorentzVector pi0_4Vec[3];
406 for (int i = 0; i < 3; i ++) pi0_4Vec[i] = TLorentzVector(0, 0, 0, 0);
407 for (int i = 0; i < 6; i ++) photon_4Vec[i] = TLorentzVector(0, 0, 0, 0);
408 int ng_max = 0;
409 if (m_decay == "pi0->2g") {
410 ng_max = 2;
411 double masses[] = {M_gamma, M_gamma};
412 if (decayGen.SetDecay(eta_LAB_4Vec, 2, masses)) {
413 decayGen.Generate();
414 photon_4Vec[0] = * decayGen.GetDecay(0);
415 photon_4Vec[1] = * decayGen.GetDecay(1);
416 }
417 } else if (m_decay == "eta->2g") {
418 ng_max = 2;
419 double masses[] = {M_gamma, M_gamma};
420 if (decayGen.SetDecay(eta_LAB_4Vec, 2, masses)) {
421 decayGen.Generate();
422 photon_4Vec[0] = * decayGen.GetDecay(0);
423 photon_4Vec[1] = * decayGen.GetDecay(1);
424 }
425 } else if (m_decay == "eta->6g") {
426 ng_max = 6;
427 double masses[] = {M_pi0, M_pi0, M_pi0};
428 if (decayGen.SetDecay(eta_LAB_4Vec, 3, masses)) {
429 decayGen.Generate();
430 pi0_4Vec[0] = * decayGen.GetDecay(0);
431 pi0_4Vec[1] = * decayGen.GetDecay(1);
432 pi0_4Vec[2] = * decayGen.GetDecay(2);
433 //for (int j = 0; j < 3; j ++) {
434 double mass[] = {M_gamma, M_gamma};
435 if (decayGenTMP1.SetDecay(pi0_4Vec[0], 2, mass)) {
436 decayGenTMP1.Generate();
437 photon_4Vec[0] = * decayGenTMP1.GetDecay(0);
438 photon_4Vec[1] = * decayGenTMP1.GetDecay(1);
439 }
440 if (decayGenTMP2.SetDecay(pi0_4Vec[1], 2, mass)) {
441 decayGenTMP2.Generate();
442 photon_4Vec[2] = * decayGenTMP2.GetDecay(0);
443 photon_4Vec[3] = * decayGenTMP2.GetDecay(1);
444 }
445 if (decayGenTMP3.SetDecay(pi0_4Vec[2], 2, mass)) {
446 decayGenTMP3.Generate();
447 photon_4Vec[4] = * decayGenTMP3.GetDecay(0);
448 photon_4Vec[5] = * decayGenTMP3.GetDecay(1);
449 }
450 }
451 }
452
453 // Deduce by energy and mass conservation the recoil nucleus 4Vec
454 TLorentzVector He4_LAB_4Vec = IS_4Vec - eta_LAB_4Vec;
455 if (m_target == "Neutron")
456 He4_LAB_4Vec = TLorentzVector(He4_LAB_4Vec.Px(), He4_LAB_4Vec.Py(), He4_LAB_4Vec.Pz(), sqrt(pow(M_n, 2) + pow(He4_LAB_4Vec.P(), 2)));
457 if (m_target == "Proton")
458 He4_LAB_4Vec = TLorentzVector(He4_LAB_4Vec.Px(), He4_LAB_4Vec.Py(), He4_LAB_4Vec.Pz(), sqrt(pow(M_p, 2) + pow(He4_LAB_4Vec.P(), 2)));
459
460 h_Tkin_recoilA_vs_egam->Fill(ebeam, He4_LAB_4Vec.E() - He4_LAB_4Vec.M());
461 h_theta_recoilA_vs_egam->Fill(ebeam, He4_LAB_4Vec.Theta() * TMath::RadToDeg());
462 h_Tkin_eta_vs_egam->Fill(ebeam, eta_LAB_4Vec.E() - eta_LAB_4Vec.M());
463 h_theta_eta_vs_egam->Fill(ebeam, eta_LAB_4Vec.Theta() * TMath::RadToDeg());
464 for (int j = 0; j < ng_max; j ++) {
465 h_Tkin_photon_vs_egam->Fill(ebeam, photon_4Vec[j].E());
466 h_theta_photon_vs_egam->Fill(ebeam, photon_4Vec[j].Theta() * TMath::RadToDeg());
467 }
468
469 if (hddmWriter) {
470 // ======= HDDM output =========
471 tmpEvt_t tmpEvt;
472 tmpEvt.str_target = m_target;
473 tmpEvt.beam = InGamma_4Vec;
474 tmpEvt.target = Target_4Vec;
475 if (m_decay == "pi0->2g") {
476 tmpEvt.q1 = photon_4Vec[0];
477 tmpEvt.q2 = photon_4Vec[1];
478 tmpEvt.q3 = He4_LAB_4Vec;
479 tmpEvt.nGen = 3;
480 } else if (m_decay == "eta->2g") {
481 tmpEvt.q1 = photon_4Vec[0];
482 tmpEvt.q2 = photon_4Vec[1];
483 tmpEvt.q3 = He4_LAB_4Vec;
484 tmpEvt.nGen = 3;
485 } else if (m_decay == "eta->6g") {
486 tmpEvt.q1 = photon_4Vec[0];
487 tmpEvt.q2 = photon_4Vec[1];
488 tmpEvt.q3 = photon_4Vec[2];
489 tmpEvt.q4 = photon_4Vec[3];
490 tmpEvt.q5 = photon_4Vec[4];
491 tmpEvt.q6 = photon_4Vec[5];
492 tmpEvt.q7 = He4_LAB_4Vec;
493 tmpEvt.nGen = 7;
494 } else if (ng_max == 0) {
495 //cout <<"decay " << m_decay << endl;
496 tmpEvt.str_meson = m_decay;
497 tmpEvt.q1 = eta_LAB_4Vec;
498 tmpEvt.q2 = He4_LAB_4Vec;
499 tmpEvt.nGen = 2;
500 }
501 tmpEvt.weight = 1.;
502 hddmWriter->write(tmpEvt,runNum,i);
503 }
504 if (asciiWriter) {
505 // ======= ASCII output =========
506 (*asciiWriter)<<runNum<<" "<<i<<" 3"<<endl;
507 // photons from the eta
508 (*asciiWriter)<<"0 "<<gamma_TYPE1<<" "<<M_gamma<<endl;
509 (*asciiWriter)<<" "<<0<<" "<<photon_4Vec[0].Px()<<" "<<photon_4Vec[0].Py()<<" "<<photon_4Vec[0].Pz()<<" "<<photon_4Vec[0].E()<<endl;
510 (*asciiWriter)<<"1 "<<gamma_TYPE1<<" "<<M_gamma<<endl;
511 (*asciiWriter)<<" "<<0<<" "<<photon_4Vec[1].Px()<<" "<<photon_4Vec[1].Py()<<" "<<photon_4Vec[1].Pz()<<" "<<photon_4Vec[1].E()<<endl;
512 // Nucleus recoil
513 if (m_target == "He4") (*asciiWriter)<<"2 "<<Helium_TYPE47<<" "<<M_He4<<endl;
514 if (m_target == "Be4") (*asciiWriter)<<"2 "<<Be9_TYPE64<<" "<<M_Be9<<endl;
515 if (m_target == "Proton") (*asciiWriter)<<"2 "<<Proton_TYPE14<<" "<<M_p<<endl;
516 if (m_target == "Neutron") (*asciiWriter)<<"2 "<<Neutron_TYPE13<<" "<<M_n<<endl;
517 (*asciiWriter)<<" "<<1<<" "<<He4_LAB_4Vec.Px()<<" "<<He4_LAB_4Vec.Py()<<" "<<He4_LAB_4Vec.Pz()<<" "<<He4_LAB_4Vec.E()<<endl;
518 }
519
520 // deletion
521 delete h_ThetaLAB;
522 }
523 h_Tkin_eta_vs_egam->Write();
524 h_Tkin_photon_vs_egam->Write();
525 h_Tkin_recoilA_vs_egam->Write();
526 h_theta_eta_vs_egam->Write();
527 h_theta_photon_vs_egam->Write();
528 h_theta_recoilA_vs_egam->Write();
529 //h_dxs->Write();
530 //cobrem_vs_Erange->Write();
531 diagOut->Close();
532
533 if (hddmWriter) delete hddmWriter;
534 if (asciiWriter) delete asciiWriter;
535
536 return 0;
537}
538
539