Bug Summary

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