1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
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 | |
26 | |
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 | |
48 | using std::complex; |
49 | using 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 | |
59 | int main( int argc, char* argv[] ){ |
60 | |
61 | string beamconfigfile(""); |
62 | TString genconfigfile(""); |
63 | string outname(""); |
64 | string hddmname(""); |
65 | |
66 | ifstream in_coherent; |
67 | ifstream in_incoherent; |
68 | |
69 | double beamLowE = 3.0; |
70 | |
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 | |
87 | |
88 | |
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 | |
97 | |
98 | |
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 | |
165 | gRandom->SetSeed(seed); |
166 | |
167 | |
168 | HddmOut *hddmWriter = nullptr; |
169 | if (hddmname != "") |
170 | hddmWriter = new HddmOut(hddmname.c_str()); |
171 | |
172 | |
173 | ofstream *asciiWriter = nullptr; |
174 | if (outname != "") |
175 | asciiWriter = new ofstream(outname.c_str()); |
176 | |
177 | |
178 | TF1 ebeam_spectrum("beam_spectrum","1/x",beamLowE,beamHighE); |
179 | |
180 | |
181 | TH1D * cobrem_vs_E = 0; |
182 | if (beamconfigfile != "") { |
183 | BeamProperties beamProp( beamconfigfile ); |
184 | cobrem_vs_E = (TH1D*)beamProp.GetFlux(); |
185 | } |
186 | |
187 | |
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"); |
| Value stored to 'm_binning' during its initialization is never read |
196 | |
197 | |
198 | |
199 | |
200 | |
201 | |
202 | |
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 | |
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 | |
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 | |
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 | |
266 | TLorentzVector InGamma_4Vec(0, 0, ebeam, ebeam); |
267 | |
268 | |
269 | TLorentzVector IS_4Vec = InGamma_4Vec + Target_4Vec; |
270 | |
271 | |
272 | double sqrt_s = IS_4Vec.M(); |
273 | double s = pow(sqrt_s, 2); |
274 | |
275 | |
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 | |
282 | |
283 | |
284 | |
285 | |
286 | |
287 | |
288 | |
289 | |
290 | |
291 | |
292 | |
293 | |
294 | |
295 | |
296 | |
297 | |
298 | |
299 | |
300 | |
301 | |
302 | |
303 | |
304 | |
305 | |
306 | |
307 | |
308 | |
309 | |
310 | |
311 | |
312 | |
313 | |
314 | |
315 | |
316 | |
317 | |
318 | |
319 | |
320 | |
321 | |
322 | |
323 | |
324 | |
325 | |
326 | |
327 | |
328 | |
329 | |
330 | |
331 | |
332 | |
333 | |
334 | |
335 | |
336 | |
337 | |
338 | |
339 | |
340 | |
341 | |
342 | |
343 | |
344 | |
345 | |
346 | |
347 | |
348 | |
349 | |
350 | |
351 | |
352 | |
353 | |
354 | |
355 | |
356 | |
357 | |
358 | |
359 | double ThetaLAB = h_ThetaLAB->GetRandom(); |
360 | ThetaLAB *= TMath::DegToRad(); |
361 | |
362 | |
363 | double PhiLAB = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); |
364 | |
365 | |
366 | |
367 | |
368 | |
369 | |
370 | |
371 | |
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 | |
389 | double P_LAB_eta = sqrt(pow(E_LAB_eta, 2) - pow(M_meson, 2)); |
390 | |
391 | |
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 | |
397 | TLorentzVector eta_LAB_4Vec(Px_LAB_eta, Py_LAB_eta, Pz_LAB_eta, E_LAB_eta); |
398 | |
399 | |
400 | |
401 | |
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 | |
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 | |
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 | |
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 | |
504 | (*asciiWriter)<<runNum<<" "<<i<<" 3"<<endl; |
505 | |
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 | |
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 | |
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 | |
528 | |
529 | diagOut->Close(); |
530 | |
531 | if (hddmWriter) delete hddmWriter; |
532 | if (asciiWriter) delete asciiWriter; |
533 | |
534 | return 0; |
535 | } |
536 | |
537 | |