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 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 | |
60 | int main( int argc, char* argv[] ){ |
61 | |
62 | string beamconfigfile(""); |
63 | TString genconfigfile(""); |
64 | string outname(""); |
65 | string hddmname(""); |
66 | |
67 | ifstream in_coherent; |
68 | ifstream in_incoherent; |
69 | |
70 | double beamLowE = 3.0; |
71 | |
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 | |
88 | |
89 | |
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 | |
98 | |
99 | |
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 | |
166 | gRandom->SetSeed(seed); |
167 | |
168 | |
169 | HddmOut *hddmWriter = nullptr; |
170 | if (hddmname != "") |
171 | hddmWriter = new HddmOut(hddmname.c_str()); |
172 | |
173 | |
174 | ofstream *asciiWriter = nullptr; |
175 | if (outname != "") |
176 | asciiWriter = new ofstream(outname.c_str()); |
177 | |
178 | |
179 | TF1 ebeam_spectrum("beam_spectrum","1/x",beamLowE,beamHighE); |
180 | |
181 | |
182 | TH1D * cobrem_vs_E = 0; |
183 | if (beamconfigfile != "") { |
184 | BeamProperties beamProp( beamconfigfile ); |
185 | cobrem_vs_E = (TH1D*)beamProp.GetFlux(); |
186 | } |
187 | |
188 | |
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"); |
| Value stored to 'm_binning' during its initialization is never read |
197 | |
198 | |
199 | |
200 | |
201 | |
202 | |
203 | |
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 | |
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 | |
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 | |
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 | |
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 | |
267 | TLorentzVector InGamma_4Vec(0, 0, ebeam, ebeam); |
268 | |
269 | |
270 | TLorentzVector IS_4Vec = InGamma_4Vec + Target_4Vec; |
271 | |
272 | |
273 | double sqrt_s = IS_4Vec.M(); |
274 | double s = pow(sqrt_s, 2); |
275 | |
276 | |
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 | |
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 | |
360 | double ThetaLAB = h_ThetaLAB->GetRandom(); |
361 | ThetaLAB *= TMath::DegToRad(); |
362 | |
363 | |
364 | double PhiLAB = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); |
365 | |
366 | |
367 | |
368 | |
369 | |
370 | |
371 | |
372 | |
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 | |
390 | double P_LAB_eta = sqrt(pow(E_LAB_eta, 2) - pow(M_meson, 2)); |
391 | |
392 | |
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 | |
398 | TLorentzVector eta_LAB_4Vec(Px_LAB_eta, Py_LAB_eta, Pz_LAB_eta, E_LAB_eta); |
399 | |
400 | |
401 | |
402 | |
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 | |
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 | |
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 | |
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 | |
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 | |
506 | (*asciiWriter)<<runNum<<" "<<i<<" 3"<<endl; |
507 | |
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 | |
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 | |
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 | |
530 | |
531 | diagOut->Close(); |
532 | |
533 | if (hddmWriter) delete hddmWriter; |
534 | if (asciiWriter) delete asciiWriter; |
535 | |
536 | return 0; |
537 | } |
538 | |
539 | |