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