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 | #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 | |
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 | double beamHighE = 12.0; |
71 | const double M_He4 = 3.727379378; |
72 | const double M_Be9 = 8.39479; |
73 | const double M_p = 0.93827208816; |
74 | const double M_n = 0.93956542052; |
75 | const double M_eta = 0.54730; |
76 | const double M_pi0 = 1.3957018e-01; |
77 | const double M_gamma = 0.0; |
78 | |
79 | int runNum = 9001; |
80 | int seed = 0; |
81 | |
82 | int nEvents = 10000; |
83 | |
84 | |
85 | |
86 | |
87 | for (int i = 1; i < argc; i++) { |
88 | |
89 | string arg(argv[i]); |
90 | |
91 | |
92 | if (arg == "-c") { |
93 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
94 | else beamconfigfile = argv[++i]; |
95 | } |
96 | if (arg == "-e") { |
97 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
98 | else genconfigfile = argv[++i]; |
99 | } |
100 | if (arg == "-o") { |
101 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
102 | else outname = argv[++i]; |
103 | } |
104 | if (arg == "-hd") { |
105 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
106 | else hddmname = argv[++i]; |
107 | } |
108 | if (arg == "-n"){ |
109 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
110 | else nEvents = atoi( argv[++i] ); |
111 | } |
112 | if (arg == "-a") { |
113 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
114 | else beamLowE = atof( argv[++i] ); |
115 | } |
116 | if (arg == "-b") { |
117 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
118 | else beamHighE = atof( argv[++i] ); |
119 | } |
120 | if (arg == "-r") { |
121 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
122 | else runNum = atoi( argv[++i] ); |
123 | } |
124 | if (arg == "-s") { |
125 | if ((i+1 == argc) || (argv[i+1][0] == '-')) arg = "-h"; |
126 | else seed = atoi( argv[++i] ); |
127 | } |
128 | if (arg == "-h") { |
129 | cout << endl << " Usage for: " << argv[0] << endl << endl; |
130 | cout << "\t -c <file>\t Beam config file" << endl; |
131 | cout << "\t -e <file>\t Generator config file" << endl; |
132 | cout << "\t -o <name>\t ASCII file output name" << endl; |
133 | cout << "\t -hd <name>\t HDDM file output name [optional]" << endl; |
134 | cout << "\t -n <value>\t Number of events to generate [optional]" << endl; |
135 | cout << "\t -a <value>\t Minimum photon energy to simulate events [optional]" << endl; |
136 | cout << "\t -b <value>\t Maximum photon energy to simulate events [optional]" << endl; |
137 | cout << "\t -r <value>\t Run number assigned to generated events [optional]" << endl; |
138 | cout << "\t -s <value>\t Random number seed initialization [optional]" << endl; |
139 | exit(1); |
140 | } |
141 | } |
142 | |
143 | if (outname.size() == 0 && hddmname == "") { |
144 | cout << "No output specificed: run gen_compton_simple -h for help" << endl; |
145 | exit(1); |
146 | } |
147 | |
148 | if (genconfigfile == "") { |
149 | cout << "No generator configuration file: run gen_primex_eta_he4 -h for help " << endl; |
150 | exit(1); |
151 | } |
152 | |
153 | gRandom->SetSeed(seed); |
154 | |
155 | |
156 | HddmOut *hddmWriter = nullptr; |
157 | if (hddmname != "") |
158 | hddmWriter = new HddmOut(hddmname.c_str()); |
159 | |
160 | |
161 | ofstream *asciiWriter = nullptr; |
162 | if (outname != "") |
163 | asciiWriter = new ofstream(outname.c_str()); |
164 | |
165 | |
166 | TF1 ebeam_spectrum("beam_spectrum","1/x",beamLowE,beamHighE); |
167 | |
168 | |
169 | TH1D * cobrem_vs_E = 0; |
170 | if (beamconfigfile != "") { |
171 | BeamProperties beamProp( beamconfigfile ); |
172 | cobrem_vs_E = (TH1D*)beamProp.GetFlux(); |
173 | } |
174 | |
175 | |
176 | MyReadConfig * ReadFile = new MyReadConfig(); |
177 | ReadFile->ReadConfigFile(genconfigfile); |
178 | TString m_rfile = ReadFile->GetConfigName("rfile"); |
179 | TString m_histo = ReadFile->GetConfigName("histo"); |
180 | TString m_decay = ReadFile->GetConfigName("decay"); |
181 | TString m_target = ReadFile->GetConfigName("target"); |
182 | TString m_Fermi_file = ReadFile->GetConfigName("Fermi_file"); |
183 | cout << "rfile " << m_rfile << endl; |
184 | cout << "histo " << m_histo << endl; |
185 | cout << "decay " << m_decay << endl; |
186 | cout << "target " << m_target << endl; |
187 | cout << "Fermi_file " << m_Fermi_file << endl; |
188 | |
189 | double M_target = M_He4; |
190 | if (m_target == "He4") M_target = M_He4; |
191 | if (m_target == "Be9") M_target = M_Be9; |
192 | if (m_target == "Proton") M_target = M_p; |
193 | if (m_target == "Neutron") M_target = M_n; |
194 | TLorentzVector Target_4Vec(0, 0, 0, M_target); |
195 | |
196 | |
197 | TFile * ifile = new TFile(m_rfile); |
198 | TH2F * h_dxs = (TH2F *) ifile->Get(m_histo); |
199 | |
200 | double M_meson = 0; |
201 | if (m_decay.Contains("eta")) |
202 | M_meson = M_eta; |
203 | else if (m_decay.Contains("pi0")) |
204 | M_meson = M_pi0; |
205 | |
206 | |
207 | TGenPhaseSpace decayGen; |
208 | TGenPhaseSpace decayGenTMP; |
209 | |
210 | TFile* diagOut = new TFile( "gen_primex_eta_he4_diagnostic.root", "recreate" ); |
211 | 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); |
212 | 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); |
213 | 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); |
214 | 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.); |
215 | 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.); |
216 | 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.); |
217 | |
218 | for (int i = 0; i < nEvents; ++i) { |
219 | if (i%1000 == 1) |
220 | cout << "event " << i <<endl; |
221 | |
222 | |
223 | double ebeam = 0; |
224 | if (beamconfigfile == "" || cobrem_vs_E == 0) |
225 | ebeam = ebeam_spectrum.GetRandom(); |
226 | else if (beamconfigfile != "") |
227 | ebeam = cobrem_vs_E->GetRandom(); |
228 | |
229 | |
230 | TLorentzVector InGamma_4Vec(0, 0, ebeam, ebeam); |
231 | |
232 | |
233 | TLorentzVector IS_4Vec = InGamma_4Vec + Target_4Vec; |
234 | |
235 | |
236 | double sqrt_s = IS_4Vec.M(); |
237 | double s = pow(sqrt_s, 2); |
| Value stored to 's' during its initialization is never read |
238 | |
239 | |
240 | int ebeam_bin = h_dxs->GetXaxis()->FindBin(ebeam); |
241 | TH1F * h_ThetaLAB = (TH1F *) h_dxs->ProjectionY("h_ThetaLAB", ebeam_bin, ebeam_bin); |
242 | if (h_ThetaLAB->GetEntries() == 0) continue; |
243 | |
244 | |
245 | |
246 | |
247 | |
248 | |
249 | |
250 | |
251 | |
252 | |
253 | |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | |
264 | |
265 | |
266 | |
267 | |
268 | |
269 | |
270 | |
271 | |
272 | |
273 | |
274 | |
275 | |
276 | |
277 | |
278 | |
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 | double ThetaLAB = h_ThetaLAB->GetRandom(); |
324 | ThetaLAB *= TMath::DegToRad(); |
325 | |
326 | |
327 | double PhiLAB = gRandom->Uniform(-TMath::Pi(), TMath::Pi()); |
328 | |
329 | |
330 | |
331 | |
332 | |
333 | |
334 | |
335 | |
336 | double E_LAB_eta = 0; |
337 | double p0 = IS_4Vec.P(); |
338 | double E0 = IS_4Vec.E(); |
339 | double m2 = M_meson; |
340 | double m3 = M_target; |
341 | double a = pow(2.0 * p0 * cos(ThetaLAB), 2) - pow(2.0 * E0, 2); |
342 | 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; |
343 | 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) + |
344 | 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); |
345 | double E_LAB_eta0 = (-b + sqrt(pow(b, 2) - 4.0 * a * c)) / (2.0 * a); |
346 | double E_LAB_eta1 = (-b - sqrt(pow(b, 2) - 4.0 * a * c)) / (2.0 * a); |
347 | if (E_LAB_eta0 < 0) E_LAB_eta = E_LAB_eta1; |
348 | if (E_LAB_eta1 < 0) E_LAB_eta = E_LAB_eta0; |
349 | if (E_LAB_eta0 > E_LAB_eta1) E_LAB_eta = E_LAB_eta0; |
350 | if (E_LAB_eta1 > E_LAB_eta0) E_LAB_eta = E_LAB_eta1; |
351 | |
352 | |
353 | double P_LAB_eta = sqrt(pow(E_LAB_eta, 2) - pow(M_eta, 2)); |
354 | |
355 | |
356 | double Px_LAB_eta = P_LAB_eta * sin(ThetaLAB) * cos(PhiLAB); |
357 | double Py_LAB_eta = P_LAB_eta * sin(ThetaLAB) * sin(PhiLAB); |
358 | double Pz_LAB_eta = P_LAB_eta * cos(ThetaLAB); |
359 | |
360 | |
361 | TLorentzVector eta_LAB_4Vec(Px_LAB_eta, Py_LAB_eta, Pz_LAB_eta, E_LAB_eta); |
362 | |
363 | |
364 | |
365 | |
366 | |
367 | TLorentzVector photon_4Vec[6]; |
368 | TLorentzVector pi0_4Vec[3]; |
369 | int ng_max = 0; |
370 | if (m_decay == "pi0->2g") { |
371 | ng_max = 2; |
372 | double masses[] = {M_gamma, M_gamma}; |
373 | if (decayGen.SetDecay(eta_LAB_4Vec, 2, masses)) { |
374 | decayGen.Generate(); |
375 | photon_4Vec[0] = * decayGen.GetDecay(0); |
376 | photon_4Vec[1] = * decayGen.GetDecay(1); |
377 | } |
378 | } else if (m_decay == "eta->2g") { |
379 | ng_max = 2; |
380 | double masses[] = {M_gamma, M_gamma}; |
381 | if (decayGen.SetDecay(eta_LAB_4Vec, 2, masses)) { |
382 | decayGen.Generate(); |
383 | photon_4Vec[0] = * decayGen.GetDecay(0); |
384 | photon_4Vec[1] = * decayGen.GetDecay(1); |
385 | } |
386 | } else if (m_decay == "eta->6g") { |
387 | ng_max = 6; |
388 | double masses[] = {M_pi0, M_pi0, M_pi0}; |
389 | if (decayGen.SetDecay(eta_LAB_4Vec, 3, masses)) { |
390 | decayGen.Generate(); |
391 | pi0_4Vec[0] = * decayGen.GetDecay(0); |
392 | pi0_4Vec[1] = * decayGen.GetDecay(1); |
393 | pi0_4Vec[2] = * decayGen.GetDecay(2); |
394 | } |
395 | for (int j = 0; j < 3; j ++) { |
396 | double mass[] = {M_gamma, M_gamma}; |
397 | if (decayGenTMP.SetDecay(pi0_4Vec[j], 2, mass)) { |
398 | decayGenTMP.Generate(); |
399 | photon_4Vec[0 + 2 * j] = * decayGenTMP.GetDecay(0); |
400 | photon_4Vec[1 + 2 * j] = * decayGenTMP.GetDecay(1); |
401 | } |
402 | } |
403 | } |
404 | |
405 | |
406 | |
407 | TLorentzVector He4_LAB_4Vec = IS_4Vec - eta_LAB_4Vec; |
408 | |
409 | h_Tkin_recoilA_vs_egam->Fill(ebeam, He4_LAB_4Vec.E() - He4_LAB_4Vec.M()); |
410 | h_theta_recoilA_vs_egam->Fill(ebeam, He4_LAB_4Vec.Theta() * TMath::RadToDeg()); |
411 | h_Tkin_eta_vs_egam->Fill(ebeam, eta_LAB_4Vec.E() - eta_LAB_4Vec.M()); |
412 | h_theta_eta_vs_egam->Fill(ebeam, eta_LAB_4Vec.Theta() * TMath::RadToDeg()); |
413 | for (int j = 0; j < ng_max; j ++) { |
414 | h_Tkin_photon_vs_egam->Fill(ebeam, photon_4Vec[j].E()); |
415 | h_theta_photon_vs_egam->Fill(ebeam, photon_4Vec[j].Theta() * TMath::RadToDeg()); |
416 | } |
417 | |
418 | if (hddmWriter) { |
419 | |
420 | tmpEvt_t tmpEvt; |
421 | tmpEvt.str_target = m_target; |
422 | tmpEvt.beam = InGamma_4Vec; |
423 | tmpEvt.target = Target_4Vec; |
424 | if (m_decay == "pi0->2g") { |
425 | tmpEvt.q1 = photon_4Vec[0]; |
426 | tmpEvt.q2 = photon_4Vec[1]; |
427 | tmpEvt.q3 = He4_LAB_4Vec; |
428 | tmpEvt.nGen = 3; |
429 | } else if (m_decay == "eta->2g") { |
430 | tmpEvt.q1 = photon_4Vec[0]; |
431 | tmpEvt.q2 = photon_4Vec[1]; |
432 | tmpEvt.q3 = He4_LAB_4Vec; |
433 | tmpEvt.nGen = 3; |
434 | } else if (m_decay == "eta->6g") { |
435 | tmpEvt.q1 = photon_4Vec[0]; |
436 | tmpEvt.q2 = photon_4Vec[1]; |
437 | tmpEvt.q3 = photon_4Vec[2]; |
438 | tmpEvt.q4 = photon_4Vec[3]; |
439 | tmpEvt.q5 = photon_4Vec[4]; |
440 | tmpEvt.q6 = photon_4Vec[5]; |
441 | tmpEvt.q7 = He4_LAB_4Vec; |
442 | tmpEvt.nGen = 7; |
443 | } else if (ng_max == 0) { |
444 | tmpEvt.q1 = eta_LAB_4Vec; |
445 | tmpEvt.q2 = He4_LAB_4Vec; |
446 | tmpEvt.nGen = 2; |
447 | } |
448 | tmpEvt.weight = 1.; |
449 | hddmWriter->write(tmpEvt,runNum,i); |
450 | } |
451 | if (asciiWriter) { |
452 | |
453 | (*asciiWriter)<<runNum<<" "<<i<<" 3"<<endl; |
454 | |
455 | (*asciiWriter)<<"0 "<<gamma_TYPE1<<" "<<M_gamma<<endl; |
456 | (*asciiWriter)<<" "<<0<<" "<<photon_4Vec[0].Px()<<" "<<photon_4Vec[0].Py()<<" "<<photon_4Vec[0].Pz()<<" "<<photon_4Vec[0].E()<<endl; |
457 | (*asciiWriter)<<"1 "<<gamma_TYPE1<<" "<<M_gamma<<endl; |
458 | (*asciiWriter)<<" "<<0<<" "<<photon_4Vec[1].Px()<<" "<<photon_4Vec[1].Py()<<" "<<photon_4Vec[1].Pz()<<" "<<photon_4Vec[1].E()<<endl; |
459 | |
460 | if (m_target == "He4") (*asciiWriter)<<"2 "<<Helium_TYPE47<<" "<<M_He4<<endl; |
461 | if (m_target == "Be4") (*asciiWriter)<<"2 "<<Be9_TYPE64<<" "<<M_Be9<<endl; |
462 | if (m_target == "Proton") (*asciiWriter)<<"2 "<<Proton_TYPE14<<" "<<M_p<<endl; |
463 | if (m_target == "Neutron") (*asciiWriter)<<"2 "<<Neutron_TYPE13<<" "<<M_n<<endl; |
464 | (*asciiWriter)<<" "<<1<<" "<<He4_LAB_4Vec.Px()<<" "<<He4_LAB_4Vec.Py()<<" "<<He4_LAB_4Vec.Pz()<<" "<<He4_LAB_4Vec.E()<<endl; |
465 | } |
466 | |
467 | |
468 | delete h_ThetaLAB; |
469 | } |
470 | h_Tkin_eta_vs_egam->Write(); |
471 | h_Tkin_photon_vs_egam->Write(); |
472 | h_Tkin_recoilA_vs_egam->Write(); |
473 | h_theta_eta_vs_egam->Write(); |
474 | h_theta_photon_vs_egam->Write(); |
475 | h_theta_recoilA_vs_egam->Write(); |
476 | |
477 | |
478 | diagOut->Close(); |
479 | |
480 | if (hddmWriter) delete hddmWriter; |
481 | if (asciiWriter) delete asciiWriter; |
482 | |
483 | return 0; |
484 | } |
485 | |
486 | |