1 | |
2 | #include "HDDM/hddm_s.h" |
3 | #include "particleType.h" |
4 | |
5 | #include <TMath.h> |
6 | #include <TRandom3.h> |
7 | #include <TGenPhaseSpace.h> |
8 | #include <TFile.h> |
9 | #include <TH1D.h> |
10 | #include <TH2D.h> |
11 | #include <TGraph.h> |
12 | #include <cstdlib> |
13 | #include <sys/stat.h> |
14 | #include <iostream> |
15 | #include <fstream> |
16 | #include <iomanip> |
17 | #include <vector> |
18 | #include <string> |
19 | #include <algorithm> |
20 | #include <cctype> |
21 | #include <locale> |
22 | using namespace std; |
23 | |
24 | #include "UTILITIES/BeamProperties.h" |
25 | #ifdef HAVE_EVTGEN1 |
26 | #include "EVTGEN_MODELS/RegisterGlueXModels.h" |
27 | |
28 | #include "EvtGen/EvtGen.hh" |
29 | |
30 | #include "EvtGenBase/EvtParticle.hh" |
31 | #include "EvtGenBase/EvtParticleFactory.hh" |
32 | #include "EvtGenBase/EvtPatches.hh" |
33 | #include "EvtGenBase/EvtPDL.hh" |
34 | #include "EvtGenBase/EvtRandom.hh" |
35 | #include "EvtGenBase/EvtReport.hh" |
36 | #include "EvtGenBase/EvtHepMCEvent.hh" |
37 | #include "EvtGenBase/EvtSimpleRandomEngine.hh" |
38 | #include "EvtGenBase/EvtMTRandomEngine.hh" |
39 | #include "EvtGenBase/EvtAbsRadCorr.hh" |
40 | #include "EvtGenBase/EvtDecayBase.hh" |
41 | |
42 | #ifdef EVTGEN_EXTERNAL1 |
43 | #include "EvtGenExternal/EvtExternalGenList.hh" |
44 | #endif //EVTGEN_EXTERNAL |
45 | #endif //HAVE_EVTGEN |
46 | |
47 | typedef struct { |
48 | bool decayed = false; |
49 | int parent_id; |
50 | vector<Particle_t> ids; |
51 | vector<TLorentzVector> p4vs; |
52 | } secondary_decay_t; |
53 | |
54 | |
55 | static inline void ltrim(std::string &s) { |
56 | s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](int ch) { |
57 | return !std::isspace(ch); |
58 | })); |
59 | } |
60 | |
61 | |
62 | static inline void rtrim(std::string &s) { |
63 | s.erase(std::find_if(s.rbegin(), s.rend(), [](int ch) { |
64 | return !std::isspace(ch); |
65 | }).base(), s.end()); |
66 | } |
67 | |
68 | |
69 | static inline void trim(std::string &s) { |
70 | ltrim(s); |
71 | rtrim(s); |
72 | } |
73 | |
74 | |
75 | const double m_p=0.93827; |
76 | const double m_p_sq=m_p*m_p; |
77 | double m_eta=0.54775; |
78 | double m_eta_sq=m_eta*m_eta; |
79 | |
80 | double width=0.; |
81 | |
82 | double g_rho_eta_gamma=0.81; |
83 | double g_omega_eta_gamma=0.29; |
84 | double g_eta_gamma_gamma=0.0429; |
85 | double g_phi_eta_gamma=0.38; |
86 | |
87 | int Nevents=10000; |
88 | int runNo=10000; |
89 | bool debug=false; |
90 | |
91 | |
92 | TH1D *thrown_t; |
93 | TH1D *thrown_dalitzZ; |
94 | TH1D *thrown_Egamma; |
95 | TH2D *thrown_dalitzXY; |
96 | TH2D *thrown_theta_vs_p; |
97 | TH2D *thrown_theta_vs_p_eta; |
98 | TH1D *cobrems_vs_E; |
99 | |
100 | char input_file_name[250]="eta548.in"; |
101 | char output_file_name[250]="eta_gen.hddm"; |
102 | |
103 | |
104 | |
105 | bool gen_uniform_t=false; |
106 | float tflat_min=100.; |
107 | float tflat_max=100.; |
108 | |
109 | void Usage(void){ |
110 | printf("genEtaRegge: generator for eta production based on Regge trajectory formalism.\n"); |
111 | printf(" Usage: genEtaRegge <options>\n"); |
112 | printf(" Options: -N<number of events> (number of events to generate)\n"); |
113 | printf(" -O<output.hddm> (default: eta_gen.hddm)\n"); |
114 | printf(" -I<input.in> (default: eta548.in)\n"); |
115 | printf(" -R<run number> (default: 10000)\n"); |
116 | printf(" -h (Print this message and exit.)\n"); |
117 | printf("Coupling constants, photon beam energy range, and eta decay products are\n"); |
118 | printf("specified in the <input.in> file.\n"); |
119 | |
120 | exit(0); |
121 | } |
122 | |
123 | |
124 | |
125 | |
126 | |
127 | void ParseCommandLineArguments(int narg, char* argv[]) |
128 | { |
129 | int seed=0; |
130 | if (narg==1){ |
131 | Usage(); |
132 | } |
133 | for(int i=1; i<narg; i++){ |
134 | char *ptr = argv[i]; |
135 | if(ptr[0] == '-'){ |
136 | switch(ptr[1]){ |
137 | case 'h': Usage(); break; |
138 | case 'I': |
139 | sscanf(&ptr[2],"%s",input_file_name); |
140 | break; |
141 | case 'O': |
142 | sscanf(&ptr[2],"%s",output_file_name); |
143 | break; |
144 | case 'N': |
145 | sscanf(&ptr[2],"%d",&Nevents); |
146 | break; |
147 | case 'R': |
148 | sscanf(&ptr[2],"%d",&runNo); |
149 | break; |
150 | case 'S': |
151 | sscanf(&ptr[2],"%d",&seed); |
152 | break; |
153 | case 'd': |
154 | debug=true; |
155 | break; |
156 | default: |
157 | break; |
158 | } |
159 | } |
160 | } |
161 | } |
162 | |
163 | |
164 | double CrossSection(double s,double t,double p_gamma,double p_eta,double theta){ |
165 | |
166 | double c_rho_p_p=0.92/137.; |
167 | double c_omega_p_p=6.44/137.; |
168 | double c_gamma_p_p=1./(137.*137.); |
169 | double c_phi_p_p=0.72/137.; |
170 | |
171 | double kappa_rho=6.1; |
172 | double kappa_gamma=1.79; |
173 | |
174 | |
175 | |
176 | double sintheta=sin(theta); |
177 | double costheta=cos(theta); |
178 | |
179 | |
180 | double q0=p_gamma-sqrt(p_eta*p_eta+m_eta_sq); |
181 | double q3=p_gamma-p_eta*costheta; |
182 | double q0_minus_q3=q0-q3; |
183 | double q0_minus_q3_sq=q0_minus_q3*q0_minus_q3; |
184 | double q0_sq=q0*q0; |
185 | double q1sq_plus_q2sq=p_eta*p_eta*sintheta*sintheta; |
186 | |
187 | double pt3=-p_gamma; |
188 | double pt0=sqrt(m_p_sq+pt3*pt3); |
189 | double pt0_minus_pt3=pt0-pt3; |
190 | double pt0_plus_pt3=pt0+pt3; |
191 | |
192 | double p_gamma_sq=p_gamma*p_gamma; |
193 | |
194 | double pt_dot_q=pt0*q0-pt3*q3; |
195 | |
196 | |
197 | double s0=1.0; |
198 | |
199 | |
200 | double a_rho=0.55+0.8*t; |
201 | double a_rho_prime=0.8; |
202 | double regge_rho=pow(s/s0,a_rho-1.)*M_PI3.14159265358979323846*a_rho_prime/(sin(M_PI3.14159265358979323846*a_rho)*TMath::Gamma(a_rho)); |
203 | |
204 | |
205 | double a_omega=0.44+0.9*t; |
206 | double a_omega_prime=0.9; |
207 | double regge_omega=pow(s/s0,a_omega-1.)*M_PI3.14159265358979323846*a_omega_prime/(sin(M_PI3.14159265358979323846*a_omega)*TMath::Gamma(a_omega)); |
208 | |
209 | |
210 | double a_phi=0.23+0.7*t; |
211 | double a_phi_prime=0.7; |
212 | double regge_phi=pow(s/s0,a_phi-1.)*M_PI3.14159265358979323846*a_phi_prime/(sin(M_PI3.14159265358979323846*a_phi)*TMath::Gamma(a_phi)); |
213 | |
214 | |
215 | |
216 | |
217 | double amp_factor_kappa0 |
218 | =8.*p_gamma_sq*(q1sq_plus_q2sq*(s+q0_minus_q3*pt0_minus_pt3) |
219 | +q0_minus_q3_sq*pt_dot_q); |
220 | double amp_factor_kappa1=32.*p_gamma_sq*m_p*(q0_minus_q3_sq*t |
221 | +2.*q0_sq*q1sq_plus_q2sq); |
222 | double amp_factor_kappa2 |
223 | =32.*p_gamma_sq*(q1sq_plus_q2sq*(2.*q0_sq*(pt_dot_q-2.*m_p_sq) |
224 | -t*pt0_plus_pt3*pt0_plus_pt3 |
225 | +4.*pt_dot_q*q0*pt0_plus_pt3) |
226 | +q0_minus_q3_sq*(t*(pt_dot_q-2.*m_p_sq) |
227 | +2.*pt_dot_q*pt_dot_q) |
228 | ); |
229 | |
230 | |
231 | double M_rho_sq=16.*M_PI3.14159265358979323846*M_PI3.14159265358979323846*(g_rho_eta_gamma*g_rho_eta_gamma/m_eta_sq)*c_rho_p_p*regge_rho*regge_rho |
232 | *(amp_factor_kappa0-kappa_rho/(4.*m_p)*amp_factor_kappa1 |
233 | +kappa_rho*kappa_rho/(16.*m_p_sq)*amp_factor_kappa2); |
234 | double M_rho=-sqrt(M_rho_sq); |
235 | |
236 | |
237 | double M_omega_sq=16.*M_PI3.14159265358979323846*M_PI3.14159265358979323846*(g_omega_eta_gamma*g_omega_eta_gamma/m_eta_sq)*c_omega_p_p*regge_omega*regge_omega |
238 | *amp_factor_kappa0; |
239 | double M_omega=-sqrt(M_omega_sq); |
240 | |
241 | |
242 | double M_phi_sq=16.*M_PI3.14159265358979323846*M_PI3.14159265358979323846*(g_phi_eta_gamma*g_phi_eta_gamma/m_eta_sq)*c_phi_p_p*regge_phi*regge_phi |
243 | *amp_factor_kappa0; |
244 | double M_phi=+sqrt(M_phi_sq); |
245 | |
246 | |
247 | double M_primakoff_sq=16.*M_PI3.14159265358979323846*M_PI3.14159265358979323846*(g_eta_gamma_gamma*g_eta_gamma_gamma/m_eta_sq)*c_gamma_p_p/(t*t) |
248 | *(amp_factor_kappa0-kappa_gamma/(4.*m_p)*amp_factor_kappa1+kappa_gamma*kappa_gamma/(16.*m_p_sq)*amp_factor_kappa2); |
249 | double M_primakoff=sqrt(M_primakoff_sq); |
250 | |
251 | |
252 | |
253 | |
254 | |
255 | |
256 | |
257 | |
258 | |
259 | |
260 | |
261 | |
262 | |
263 | |
264 | double pi_a_omega=M_PI3.14159265358979323846*a_omega; |
265 | double pi_a_rho=M_PI3.14159265358979323846*a_rho; |
266 | double pi_a_phi=M_PI3.14159265358979323846*a_phi; |
267 | double M_sq =M_omega_sq+M_rho_sq+M_primakoff_sq+M_phi_sq |
268 | +2.*M_omega*M_phi*cos(pi_a_omega-pi_a_phi) |
269 | +2.*M_omega*M_rho*cos(pi_a_omega-pi_a_rho) |
270 | +2.*M_omega*M_primakoff*cos(pi_a_omega) |
271 | +2.*M_rho*M_primakoff*cos(pi_a_rho) |
272 | +2.*M_rho*M_phi*cos(pi_a_rho-pi_a_phi) |
273 | +2.*M_phi*M_primakoff*cos(pi_a_phi) |
274 | ; |
275 | |
276 | double hbarc_sq=389.; |
277 | double dsigma_dt=hbarc_sq*M_sq/(4.*64.*M_PI3.14159265358979323846*s*p_gamma_sq); |
278 | |
279 | |
280 | return(dsigma_dt); |
281 | |
282 | |
283 | } |
284 | |
285 | |
286 | void WriteEvent(unsigned int eventNumber,TLorentzVector &beam, float vert[3], |
287 | vector<Particle_t> &particle_types, |
288 | vector<TLorentzVector> &particle_vectors, |
289 | vector<bool> &particle_decayed, |
290 | vector< secondary_decay_t > &secondary_vertices, |
291 | s_iostream_t *file){ |
292 | s_PhysicsEvents_t* pes; |
293 | s_Reactions_t* rs; |
294 | s_Target_t* ta; |
295 | s_Beam_t* be; |
296 | s_Vertices_t* vs; |
297 | s_Origin_t* origin; |
298 | s_Products_t* ps; |
299 | s_HDDM_t *thisOutputEvent = make_s_HDDM(); |
300 | thisOutputEvent->physicsEvents = pes = make_s_PhysicsEvents(1); |
301 | pes->mult = 1; |
302 | pes->in[0].runNo = runNo; |
303 | pes->in[0].eventNo = eventNumber; |
304 | pes->in[0].reactions = rs = make_s_Reactions(1); |
305 | rs->mult = 1; |
306 | |
307 | rs->in[0].beam = be = make_s_Beam(); |
308 | be->type = Gamma; |
309 | be->properties = make_s_Properties(); |
310 | be->properties->charge = ParticleCharge(be->type); |
311 | be->properties->mass = ParticleMass(be->type); |
312 | be->momentum = make_s_Momentum(); |
313 | be->momentum->px = 0.; |
314 | be->momentum->py = 0.; |
315 | be->momentum->pz = beam.Pz(); |
316 | be->momentum->E = beam.E(); |
317 | |
318 | rs->in[0].target = ta = make_s_Target(); |
319 | ta->type = Proton; |
320 | ta->properties = make_s_Properties(); |
321 | ta->properties->charge = ParticleCharge(ta->type); |
322 | ta->properties->mass = ParticleMass(ta->type); |
323 | ta->momentum = make_s_Momentum(); |
324 | ta->momentum->px = 0.; |
325 | ta->momentum->py = 0.; |
326 | ta->momentum->pz = 0.; |
327 | ta->momentum->E = ParticleMass(ta->type); |
328 | |
329 | int num_vertices = 1 + secondary_vertices.size(); |
330 | rs->in[0].vertices = vs = make_s_Vertices(num_vertices); |
331 | vs->mult = num_vertices; |
332 | vs->in[0].origin = origin = make_s_Origin(); |
333 | vs->in[0].products = ps = make_s_Products(particle_vectors.size()); |
334 | ps->mult = 0; |
335 | origin->t = 0.0; |
336 | origin->vx = vert[0]; |
337 | origin->vy = vert[1]; |
338 | origin->vz = vert[2]; |
339 | |
340 | int part_ind = 1; |
341 | for (unsigned int i=0;i<particle_vectors.size();i++,ps->mult++){ |
342 | Particle_t my_particle=particle_types[i]; |
343 | if(particle_decayed[i]) |
344 | ps->in[ps->mult].type = Unknown; |
345 | else |
346 | ps->in[ps->mult].type = my_particle; |
347 | ps->in[ps->mult].pdgtype = PDGtype(my_particle); |
348 | ps->in[ps->mult].id = part_ind++; |
349 | ps->in[ps->mult].parentid = 0; |
350 | ps->in[ps->mult].mech = 0; |
351 | ps->in[ps->mult].momentum = make_s_Momentum(); |
352 | ps->in[ps->mult].momentum->px = particle_vectors[i].Px(); |
353 | ps->in[ps->mult].momentum->py = particle_vectors[i].Py(); |
354 | ps->in[ps->mult].momentum->pz = particle_vectors[i].Pz(); |
355 | ps->in[ps->mult].momentum->E = particle_vectors[i].E(); |
356 | } |
357 | |
358 | |
359 | if(secondary_vertices.size() > 0) { |
360 | int vertex_ind = 1; |
361 | for(auto& the_vertex : secondary_vertices) { |
362 | |
363 | vs->in[vertex_ind].origin = origin = make_s_Origin(); |
364 | vs->in[vertex_ind].products = ps = make_s_Products(the_vertex.ids.size()); |
365 | ps->mult = 0; |
366 | origin->t = 0.0; |
367 | origin->vx = vert[0]; |
368 | origin->vy = vert[1]; |
369 | origin->vz = vert[2]; |
370 | |
371 | |
372 | for (unsigned int i=0; i<the_vertex.ids.size(); i++,ps->mult++){ |
373 | Particle_t my_particle = the_vertex.ids[i]; |
374 | ps->in[ps->mult].decayVertex = vertex_ind; |
375 | ps->in[ps->mult].type = my_particle; |
376 | ps->in[ps->mult].pdgtype = PDGtype(my_particle); |
377 | ps->in[ps->mult].id = part_ind++; |
378 | ps->in[ps->mult].parentid = the_vertex.parent_id; |
379 | ps->in[ps->mult].mech = 0; |
380 | ps->in[ps->mult].momentum = make_s_Momentum(); |
381 | ps->in[ps->mult].momentum->px = the_vertex.p4vs[i].Px(); |
382 | ps->in[ps->mult].momentum->py = the_vertex.p4vs[i].Py(); |
383 | ps->in[ps->mult].momentum->pz = the_vertex.p4vs[i].Pz(); |
384 | ps->in[ps->mult].momentum->E = the_vertex.p4vs[i].E(); |
385 | } |
386 | |
387 | vertex_ind++; |
388 | } |
389 | } |
390 | flush_s_HDDM(thisOutputEvent,file); |
391 | } |
392 | |
393 | |
394 | void CreateHistograms(string beamConfigFile,int num_decay_particles){ |
395 | |
396 | if(gen_uniform_t) thrown_t=new TH1D("thrown_t","Thrown -t distribution",1000,0.,tflat_max); |
397 | else thrown_t=new TH1D("thrown_t","Thrown -t distribution",1000,0.,3); |
398 | thrown_t->SetXTitle("-t [GeV^{2}]"); |
399 | thrown_Egamma=new TH1D("thrown_Egamma","Thrown E_{#gamma} distribution", |
400 | 1000,0,12.); |
401 | thrown_Egamma->SetTitle("E_{#gamma} [GeV]"); |
402 | |
403 | thrown_theta_vs_p=new TH2D("thrown_theta_vs_p","Proton #theta_{LAB} vs. p", |
404 | 200,0,2.,180,0.,90.); |
405 | thrown_theta_vs_p->SetXTitle("p [GeV/c]"); |
406 | thrown_theta_vs_p->SetYTitle("#theta [degrees]"); |
407 | |
408 | thrown_theta_vs_p_eta=new TH2D("thrown_theta_vs_p_eta","#eta #theta_{LAB} vs. p", |
409 | 120,0,12.,180,0.,180.); |
410 | thrown_theta_vs_p_eta->SetXTitle("p [GeV/c]"); |
411 | thrown_theta_vs_p_eta->SetYTitle("#theta [degrees]"); |
412 | |
413 | if(num_decay_particles==3) { |
414 | thrown_dalitzZ=new TH1D("thrown_dalitzZ","thrown dalitz Z",110,-0.05,1.05); |
415 | thrown_dalitzXY=new TH2D("thrown_dalitzXY","Dalitz distribution Y vs X",100,-1.,1.,100,-1.,1); |
416 | } |
417 | |
418 | BeamProperties beamProp(beamConfigFile); |
419 | cobrems_vs_E = (TH1D*)beamProp.GetFlux(); |
420 | } |
421 | |
422 | |
423 | |
424 | void GraphCrossSection(double &xsec_max){ |
425 | |
426 | double Egamma=cobrems_vs_E->GetBinLowEdge(1); |
427 | |
428 | |
429 | double s=m_p*(m_p+2.*Egamma); |
430 | double Ecm=sqrt(s); |
431 | |
432 | |
433 | double p_gamma=(s-m_p_sq)/(2.*Ecm); |
434 | double E_eta=(s+m_eta_sq-m_p_sq)/(2.*Ecm); |
435 | double p_eta=sqrt(E_eta*E_eta-m_eta_sq); |
436 | |
437 | |
438 | double p_diff=p_gamma-p_eta; |
439 | double t0=m_eta_sq*m_eta_sq/(4.*s)-p_diff*p_diff; |
440 | |
441 | double sum=0.; |
442 | double t_old=t0; |
443 | double t_array[10000]; |
444 | double xsec_array[10000]; |
445 | xsec_max=0.; |
446 | for (unsigned int k=0;k<10000;k++){ |
447 | double theta_cm=M_PI3.14159265358979323846*double(k)/10000.; |
448 | double sin_theta_over_2=sin(0.5*theta_cm); |
449 | double t=t0-4.*p_gamma*p_eta*sin_theta_over_2*sin_theta_over_2; |
450 | double xsec=CrossSection(s,t,p_gamma,p_eta,theta_cm); |
451 | if (xsec>xsec_max) xsec_max=xsec; |
452 | |
453 | t_array[k]=-t; |
454 | xsec_array[k]=xsec; |
455 | |
456 | sum-=xsec*(t-t_old); |
457 | t_old=t; |
458 | } |
459 | TGraph *Gxsec=new TGraph(10000,t_array,xsec_array); |
460 | TString xsec_title = "#eta Cross Section at E_{#gamma}="+to_string(Egamma)+" GeV;-t [GeV^{2}];d#sigma/dt [#mub/GeV^{2}]"; |
461 | Gxsec->SetTitle(xsec_title); |
462 | Gxsec->Write("Cross section"); |
463 | |
464 | cout << "Total cross section at " << Egamma << " GeV = "<< sum |
465 | << " micro-barns"<<endl; |
466 | } |
467 | |
468 | |
469 | |
470 | |
471 | |
472 | int main(int narg, char *argv[]) |
473 | { |
474 | ParseCommandLineArguments(narg, argv); |
475 | |
476 | |
477 | |
478 | string rootfilename="eta_gen.root"; |
479 | TFile *rootfile=new TFile(rootfilename.c_str(),"RECREATE", |
480 | "Produced by genEta"); |
481 | |
482 | |
483 | s_iostream_t *file = init_s_HDDM(output_file_name); |
484 | |
485 | |
486 | |
487 | TRandom3 *myrand=new TRandom3(0); |
488 | |
489 | |
490 | TLorentzVector target(0.,0.,0.,m_p); |
491 | |
492 | |
493 | |
494 | |
495 | |
496 | |
497 | ifstream infile(input_file_name); |
498 | if (!infile.is_open()){ |
499 | cerr << "Input file missing! Exiting..." <<endl; |
500 | exit(-1); |
501 | } |
502 | |
503 | |
504 | string comment_line; |
505 | getline(infile,comment_line); |
506 | string beamConfigFile; |
507 | infile >> beamConfigFile; |
508 | infile.ignore(); |
509 | |
510 | cout << "Photon beam configuration file " << beamConfigFile.data() << endl; |
511 | |
512 | |
513 | string comment_line2; |
514 | getline(infile,comment_line); |
515 | double m_eta_R=0.; |
516 | infile >> m_eta_R; |
517 | infile >> width; |
518 | infile.ignore(); |
519 | |
520 | m_eta=m_eta_R; |
521 | m_eta_sq=m_eta*m_eta; |
522 | |
523 | cout << "Mass, width of decaying particle [GeV] = "<< m_eta <<"," << width << endl; |
524 | |
525 | |
526 | getline(infile,comment_line); |
527 | infile >> g_eta_gamma_gamma; |
528 | infile >> g_rho_eta_gamma; |
529 | infile >> g_omega_eta_gamma; |
530 | infile >> g_phi_eta_gamma; |
531 | infile.ignore(); |
532 | |
533 | cout << "Coupling constants:" <<endl; |
534 | cout << " g_eta_gamma_gamma = " << g_eta_gamma_gamma <<endl; |
535 | cout << " g_rho_eta_gamma = " << g_rho_eta_gamma <<endl; |
536 | cout << " g_omega_eta_gamma = " << g_omega_eta_gamma << endl; |
537 | cout << " g_phi_eta_gamma = " << g_phi_eta_gamma <<endl; |
538 | |
539 | |
540 | int num_decay_particles=0; |
541 | getline(infile,comment_line); |
542 | infile >> num_decay_particles; |
543 | infile.ignore(); |
544 | |
545 | cout << "number of decay particles = " << num_decay_particles << endl; |
546 | |
547 | |
548 | |
549 | bool use_evtgen = false; |
550 | #ifdef HAVE_EVTGEN1 |
551 | |
552 | EvtParticle* parent(0); |
553 | EvtGen *myGenerator = nullptr; |
554 | EvtId EtaId; |
555 | |
556 | |
557 | |
558 | if(num_decay_particles <= 0) { |
559 | num_decay_particles = 0; |
560 | use_evtgen = true; |
561 | |
562 | cout << "Using EvtGen to decay particles ..." << endl; |
563 | |
564 | |
565 | getline(infile,comment_line); |
566 | string particle_type; |
567 | getline(infile,particle_type); |
568 | trim(particle_type); |
569 | cout << "Generating particle: " << particle_type << endl; |
570 | |
571 | |
572 | const char* evtgen_home_env_ptr = std::getenv("EVTGENDIR"); |
573 | string EVTGEN_HOME = (evtgen_home_env_ptr==nullptr) ? "." : evtgen_home_env_ptr; |
574 | |
575 | |
576 | EvtRandomEngine* eng = 0; |
577 | #ifdef EVTGEN_CPP11 |
578 | |
579 | eng = new EvtMTRandomEngine(); |
580 | #else |
581 | eng = new EvtSimpleRandomEngine(); |
582 | #endif |
583 | |
584 | EvtRandom::setRandomEngine(eng); |
585 | |
586 | EvtAbsRadCorr* radCorrEngine = 0; |
587 | std::list<EvtDecayBase*> extraModels; |
588 | |
589 | #ifdef EVTGEN_EXTERNAL1 |
590 | bool convertPythiaCodes(false); |
591 | bool useEvtGenRandom(true); |
592 | EvtExternalGenList genList(convertPythiaCodes, "", "gamma", useEvtGenRandom); |
593 | radCorrEngine = genList.getPhotosModel(); |
594 | extraModels = genList.getListOfModels(); |
595 | #endif |
596 | |
597 | |
598 | const char* evtgen_decay_file_ptr = std::getenv("EVTGEN_DECAY_FILE"); |
599 | string evtGenDecayFile = (evtgen_decay_file_ptr==nullptr) ? EVTGEN_HOME + "/DECAY.DEC" : evtgen_decay_file_ptr; |
600 | const char* evtgen_particle_defs_ptr = std::getenv("EVTGEN_PARTICLE_DEFINITIONS"); |
601 | string evtGenParticleDefs = (evtgen_particle_defs_ptr==nullptr) ? EVTGEN_HOME + "/evt.pdl" : evtgen_particle_defs_ptr; |
602 | myGenerator = new EvtGen(evtGenDecayFile.c_str(), evtGenParticleDefs.c_str(), eng, |
603 | radCorrEngine, &extraModels); |
604 | |
605 | |
606 | GlueX_EvtGen::RegisterGlueXModels(); |
607 | |
608 | |
609 | struct stat buffer; |
610 | if(stat("userDecay.dec", &buffer) == 0) |
611 | myGenerator->readUDecay("userDecay.dec"); |
612 | |
613 | |
614 | EtaId = EvtPDL::getId(std::string(particle_type)); |
615 | |
616 | } |
617 | #endif // HAVE_EVTGEN |
618 | |
619 | |
620 | vector<Particle_t>particle_types; |
621 | |
622 | vector<double> decay_masses(num_decay_particles); |
623 | double *res_decay_masses=NULL__null; |
624 | vector<Particle_t>res_particle_types; |
625 | |
626 | |
627 | getline(infile,comment_line); |
628 | cout << comment_line << endl; |
629 | int reson_index=-1; |
630 | if(!use_evtgen) { |
631 | cout << "Particle id's of decay particles ="; |
632 | for (int k=0;k<num_decay_particles;k++){ |
633 | int ipart; |
634 | infile >> ipart; |
635 | cout << " " << ipart; |
636 | particle_types.push_back((Particle_t)ipart); |
637 | if (ipart>0){ |
638 | decay_masses[k]=ParticleMass((Particle_t)ipart); |
639 | } |
640 | else { |
641 | reson_index=k; |
642 | } |
643 | } |
644 | cout << endl; |
645 | } |
646 | unsigned int num_res_decay_particles=0; |
647 | double reson_mass=0.,reson_width=0.; |
648 | int reson_L=0; |
649 | if (reson_index>=0){ |
650 | infile.ignore(); |
651 | getline(infile,comment_line); |
652 | cout << comment_line << endl; |
653 | infile >> reson_mass; |
654 | decay_masses[reson_index]=reson_mass; |
655 | cout << "Resonance mass = " << reson_mass << " [GeV]" << endl; |
656 | infile >> reson_width; |
657 | cout << "Resonance width = " << reson_width << " [GeV]" << endl; |
658 | infile >> reson_L; |
659 | cout << "Resonance orbital angular momentum L = " << reson_L << endl; |
660 | infile >> num_res_decay_particles; |
661 | if (num_res_decay_particles>1) res_decay_masses=new double[num_res_decay_particles]; |
662 | else{ |
663 | cout << "Invalid number of decay particles! " << endl; |
664 | exit(0); |
665 | } |
666 | cout << " Decay particles: "; |
667 | for (unsigned int i=0;i<num_res_decay_particles;i++){ |
668 | int ipart; |
669 | infile >> ipart; |
670 | cout << " " << ipart; |
671 | res_particle_types.push_back((Particle_t)ipart); |
672 | res_decay_masses[i]=ParticleMass((Particle_t)ipart); |
673 | } |
674 | } |
675 | |
676 | |
677 | while( !infile.eof() ) { |
678 | string line = ""; |
679 | string tflat_string = ""; |
680 | getline(infile,line); |
681 | if(line.length() < 11) continue; |
682 | |
683 | if(line.substr(0,9) == "tflat_min") { |
684 | string str_tval=""; |
685 | for(size_t loc_i=9; loc_i<line.length(); loc_i++) { |
686 | string this_char; this_char += line[loc_i]; |
687 | if(isdigit(line[loc_i]) || this_char=="." ) { |
688 | str_tval+=line[loc_i]; |
689 | } |
690 | if(str_tval.length()>0 && this_char==" ") break; |
691 | } |
692 | tflat_min = -1*fabs(atof(str_tval.c_str())); |
693 | } |
694 | |
695 | if(line.substr(0,9) == "tflat_max") { |
696 | string str_tval=""; |
697 | for(size_t loc_i=9; loc_i<line.length(); loc_i++) { |
698 | string this_char; this_char += line[loc_i]; |
699 | if(isdigit(line[loc_i]) || this_char=="." ) { |
700 | str_tval+=line[loc_i]; |
701 | } |
702 | if(str_tval.length()>0 && this_char==" ") break; |
703 | } |
704 | tflat_max = -1*fabs(atof(str_tval.c_str())); |
705 | } |
706 | } |
707 | if(tflat_min<0.&&tflat_max<0.&&tflat_max<tflat_min) gen_uniform_t = true; |
708 | if(gen_uniform_t) cout << "GENERATING DATA WITH UNIFORM T-DIST FROM " << tflat_min << " TO " << tflat_max << endl; |
709 | |
710 | |
711 | infile.close(); |
712 | |
713 | |
714 | CreateHistograms(beamConfigFile,num_decay_particles); |
715 | |
716 | |
717 | double xsec_max=0.; |
718 | GraphCrossSection(xsec_max); |
719 | |
720 | |
721 | |
722 | |
723 | for (int i=1;i<=Nevents;i++){ |
724 | double Egamma=0.; |
725 | double xsec=0.,xsec_test=0.; |
726 | |
727 | |
728 | double theta_cm=0.; |
729 | |
730 | |
731 | double p_eta=0.; |
732 | |
733 | |
734 | double t=0.; |
735 | |
736 | |
737 | float vert[4]={0.,0.,0.,0.}; |
738 | |
739 | |
740 | do{ |
741 | |
742 | Egamma = cobrems_vs_E->GetRandom(); |
743 | |
744 | |
745 | double s=m_p*(m_p+2.*Egamma); |
746 | double Ecm=sqrt(s); |
747 | |
748 | |
749 | double p_gamma=(s-m_p_sq)/(2.*Ecm); |
750 | |
751 | |
752 | |
753 | if(!use_evtgen) { |
754 | double mass_check=0.; |
755 | do { |
756 | if (reson_index>-1 && reson_width>0.){ |
757 | if (num_res_decay_particles==2){ |
758 | double BW=0.,BWtest=0.; |
759 | double m1sq=res_decay_masses[0]*res_decay_masses[0]; |
760 | double m2sq=res_decay_masses[1]*res_decay_masses[1]; |
761 | double m0sq=reson_mass*reson_mass; |
762 | double m1sq_minus_m2sq=m1sq-m2sq; |
763 | double q0sq=(m0sq*m0sq-2.*m0sq*(m1sq+m2sq) |
764 | +m1sq_minus_m2sq*m1sq_minus_m2sq)/(4.*m0sq); |
765 | double BlattWeisskopf=1.; |
766 | double d=5.; |
767 | double dsq=d*d; |
768 | double Gamma0sq=reson_width*reson_width; |
769 | double BWmax=0.; |
770 | double BWmin=0.; |
771 | double m_min=res_decay_masses[0]+res_decay_masses[1]; |
772 | |
773 | |
774 | if (reson_L==0){ |
775 | BWmax=1./(m0sq*Gamma0sq); |
776 | } |
777 | else{ |
778 | BWmax=pow(q0sq,2*reson_L)/(m0sq*Gamma0sq); |
779 | } |
780 | double m_max=m_p*(sqrt(1.+2.*Egamma/m_p)-1.); |
781 | for (int im=0;im<num_decay_particles;im++){ |
782 | if (im==reson_index) continue; |
783 | m_max-=decay_masses[im]; |
784 | } |
785 | double m=0.; |
786 | do { |
787 | m=m_min+myrand->Uniform(m_max-m_min); |
788 | double msq=m*m; |
789 | double qsq=(msq*msq-2.*msq*(m1sq+m2sq) |
790 | +m1sq_minus_m2sq*m1sq_minus_m2sq)/(4.*msq); |
791 | double z=dsq*qsq; |
792 | double z0=dsq*q0sq; |
793 | double m0sq_minus_msq=m0sq-msq; |
794 | if (reson_L==0){ |
795 | BW=1./(m0sq_minus_msq*m0sq_minus_msq+m0sq*qsq/q0sq*Gamma0sq); |
796 | } |
797 | else{ |
798 | if (reson_L==1){ |
799 | BlattWeisskopf=(1.+z0)/(1.+z); |
800 | } |
801 | BW=pow(qsq,2*reson_L) |
802 | /(m0sq_minus_msq*m0sq_minus_msq |
803 | +m0sq*pow(qsq/q0sq,2*reson_L+1)*Gamma0sq*BlattWeisskopf); |
804 | } |
805 | BWtest=BWmin+myrand->Uniform(BWmax-BWmin); |
806 | } while (BWtest>BW); |
807 | decay_masses[reson_index]=m; |
808 | } |
809 | } |
810 | |
811 | if (width>0.001){ |
812 | |
813 | |
814 | |
815 | double m_max_=m_p*(sqrt(1.+2.*Egamma/m_p)-1.); |
816 | double m_min_=decay_masses[0]; |
817 | double m1sq_=decay_masses[0]*decay_masses[0]; |
818 | double m2sq_=0.; |
819 | switch(num_decay_particles){ |
820 | case 2: |
821 | { |
822 | m_min_+=decay_masses[1]; |
823 | m2sq_=decay_masses[1]*decay_masses[1]; |
824 | break; |
825 | } |
826 | case 3: |
827 | |
828 | |
829 | |
830 | |
831 | |
832 | |
833 | { |
834 | double m2_=0.5*(m_eta_R-decay_masses[0]+decay_masses[1]+decay_masses[2]); |
835 | m_min_+=decay_masses[1]+decay_masses[2]; |
836 | m2sq_=m2_*m2_; |
837 | break; |
838 | } |
839 | default: |
840 | break; |
841 | } |
842 | double m0sq_=m_eta_R*m_eta_R; |
843 | double BW_=0.,BWtest_=0.; |
844 | double Gamma0sq_=width*width; |
845 | double m1sq_minus_m2sq_=m1sq_-m2sq_; |
846 | double q0sq_=(m0sq_*m0sq_-2.*m0sq_*(m1sq_+m2sq_) |
847 | +m1sq_minus_m2sq_*m1sq_minus_m2sq_)/(4.*m0sq_); |
848 | double BWmax_=1./(Gamma0sq_*m0sq_); |
849 | double BWmin_=0.; |
850 | double m_=0.; |
851 | do{ |
852 | m_=m_min_+myrand->Uniform(m_max_-m_min_); |
853 | double msq_=m_*m_; |
854 | double qsq_=(msq_*msq_-2.*msq_*(m1sq_+m2sq_) |
855 | +m1sq_minus_m2sq_*m1sq_minus_m2sq_)/(4.*msq_); |
856 | double m0sq_minus_msq_=m0sq_-msq_; |
857 | BW_=1./(m0sq_minus_msq_*m0sq_minus_msq_ |
858 | +m0sq_*m0sq_*Gamma0sq_*qsq_/q0sq_); |
859 | BWtest_=BWmin_+myrand->Uniform(BWmax_-BWmin_); |
860 | } |
861 | while (BWtest_>BW_); |
862 | m_eta=m_; |
863 | m_eta_sq=m_*m_; |
864 | } |
865 | |
866 | |
867 | mass_check=decay_masses[0]; |
868 | for (int im=1;im<num_decay_particles;im++){ |
869 | mass_check+=decay_masses[im]; |
870 | } |
871 | } while (mass_check>m_eta); |
872 | } else { |
873 | |
874 | |
875 | if(width > 0.000001) { |
876 | bool is_good_mass = false; |
877 | do { |
878 | double m_ = myrand->BreitWigner(m_eta_R,width); |
879 | cout << m_ << " " << m_eta_R << " " << width << endl; |
880 | |
881 | is_good_mass = (m_ > m_eta_R-5.*width) && (m_ < m_eta_R+5.*width); |
882 | } while (!is_good_mass); |
883 | } else { |
884 | m_eta = m_eta_R; |
885 | } |
886 | } |
887 | |
888 | |
889 | double E_eta=(s+m_eta_sq-m_p_sq)/(2.*Ecm); |
890 | p_eta=sqrt(E_eta*E_eta-m_eta_sq); |
891 | |
892 | |
893 | double p_diff=p_gamma-p_eta; |
894 | double t0=m_eta_sq*m_eta_sq/(4.*s)-p_diff*p_diff; |
895 | double sin_theta_over_2=0.; |
896 | t=t0; |
| Value stored to 't' is never read |
897 | |
898 | |
899 | |
900 | double cos_theta_cm=-1.0+myrand->Uniform(2.); |
901 | theta_cm=acos(cos_theta_cm); |
902 | |
903 | sin_theta_over_2=sin(0.5*theta_cm); |
904 | t=t0-4.*p_gamma*p_eta*sin_theta_over_2*sin_theta_over_2; |
905 | xsec=CrossSection(s,t,p_gamma,p_eta,theta_cm); |
906 | |
907 | |
908 | if(gen_uniform_t&&t<tflat_min&&t>tflat_max) { |
909 | |
910 | double t_tmp = tflat_min; |
911 | |
912 | double theta_cm_tmp = 2.*asin(0.5*sqrt( (t0-t_tmp)/(p_gamma*p_eta) ) ); |
913 | xsec=CrossSection(s,t_tmp,p_gamma,p_eta,theta_cm_tmp); |
914 | |
915 | t=myrand->Uniform(tflat_max, min( float(t0),tflat_min) ); |
916 | theta_cm=2.*asin(0.5*sqrt( (t0-t)/(p_gamma*p_eta) ) ); |
917 | if( std::isnan(theta_cm)==true ) xsec=-1.; |
918 | } |
919 | |
920 | |
921 | xsec_test=myrand->Uniform(xsec_max); |
922 | } |
923 | while (xsec_test>xsec); |
924 | |
925 | |
926 | double phi_cm=myrand->Uniform(2.*M_PI3.14159265358979323846); |
927 | |
928 | |
929 | TLorentzVector beam(0.,0.,Egamma,Egamma); |
930 | thrown_Egamma->Fill(Egamma); |
931 | |
932 | |
933 | TVector3 v_cm=(1./(Egamma+m_p))*beam.Vect(); |
934 | |
935 | double pt=p_eta*sin(theta_cm); |
936 | TLorentzVector eta4(pt*cos(phi_cm),pt*sin(phi_cm),p_eta*cos(theta_cm), |
937 | sqrt(p_eta*p_eta+m_eta_sq)); |
938 | |
939 | |
940 | eta4.Boost(v_cm); |
941 | |
942 | |
943 | TLorentzVector proton4=beam+target-eta4; |
944 | |
945 | |
946 | thrown_theta_vs_p->Fill(proton4.P(),180./M_PI3.14159265358979323846*proton4.Theta()); |
947 | thrown_theta_vs_p_eta->Fill(eta4.P(),180./M_PI3.14159265358979323846*eta4.Theta()); |
948 | |
949 | |
950 | thrown_t->Fill(-t); |
951 | |
952 | |
953 | vector<TLorentzVector>output_particle_vectors; |
954 | output_particle_vectors.push_back(proton4); |
955 | |
956 | vector<Particle_t>output_particle_types; |
957 | output_particle_types.push_back(Proton); |
958 | |
959 | vector<bool>output_particle_decays; |
960 | output_particle_decays.push_back(false); |
961 | |
962 | vector<secondary_decay_t>secondary_vertices; |
963 | #ifdef HAVE_EVTGEN1 |
964 | if(use_evtgen) { |
965 | |
966 | EvtVector4R pInit(eta4.E(), eta4.X(), eta4.Y(), eta4.Z()); |
967 | parent = EvtParticleFactory::particleFactory(EtaId, pInit); |
968 | |
969 | if(num_decay_particles == 0) { |
970 | |
971 | |
972 | TLorentzVector vec4v( eta4.X(), eta4.Y(), eta4.Z(), eta4.E()); |
973 | |
974 | output_particle_vectors.push_back(vec4v); |
975 | output_particle_types.push_back(PDGtoPType(parent->getPDGId())); |
976 | } else { |
977 | |
978 | myGenerator->generateDecay(parent); |
979 | |
980 | |
981 | for(unsigned int i=0; i<parent->getNDaug(); i++) { |
982 | TLorentzVector vec4v( parent->getDaug(i)->getP4Lab().get(1), |
983 | parent->getDaug(i)->getP4Lab().get(2), |
984 | parent->getDaug(i)->getP4Lab().get(3), |
985 | parent->getDaug(i)->getP4Lab().get(0) ); |
986 | output_particle_vectors.push_back(vec4v); |
987 | output_particle_types.push_back(PDGtoPType(parent->getDaug(i)->getPDGId())); |
988 | |
989 | |
990 | |
991 | if(parent->getDaug(i)->getNDaug()>0) { |
992 | output_particle_decays.push_back(true); |
993 | |
994 | secondary_decay_t secondary_vertex; |
995 | secondary_vertex.parent_id = i; |
996 | for(unsigned int j=0; j<parent->getDaug(i)->getNDaug(); j++) { |
997 | TLorentzVector vec4v( parent->getDaug(i)->getDaug(j)->getP4Lab().get(1), |
998 | parent->getDaug(i)->getDaug(j)->getP4Lab().get(2), |
999 | parent->getDaug(i)->getDaug(j)->getP4Lab().get(3), |
1000 | parent->getDaug(i)->getDaug(j)->getP4Lab().get(0) ); |
1001 | secondary_vertex.p4vs.push_back(vec4v); |
1002 | secondary_vertex.ids.push_back(PDGtoPType(parent->getDaug(i)->getDaug(j)->getPDGId())); |
1003 | } |
1004 | secondary_vertices.push_back(secondary_vertex); |
1005 | } else { |
1006 | output_particle_decays.push_back(false); |
1007 | } |
1008 | } |
1009 | |
1010 | parent->deleteTree(); |
1011 | } |
1012 | } else { |
1013 | #endif //HAVE_EVTGEN |
1014 | |
1015 | TGenPhaseSpace phase_space; |
1016 | phase_space.SetDecay(eta4,num_decay_particles,decay_masses.data()); |
1017 | double weight=0.,rand_weight=1.; |
1018 | do{ |
1019 | weight=phase_space.Generate(); |
1020 | rand_weight=myrand->Uniform(1.); |
1021 | } |
1022 | while (rand_weight>weight); |
1023 | |
1024 | |
1025 | if (num_decay_particles==3){ |
1026 | TLorentzVector one=*phase_space.GetDecay(0); |
1027 | TLorentzVector two=*phase_space.GetDecay(1); |
1028 | TLorentzVector three=*phase_space.GetDecay(2); |
1029 | TLorentzVector one_two=one+two; |
1030 | TLorentzVector one_three=one+three; |
1031 | |
1032 | TLorentzVector eta=one_two+three; |
1033 | TVector3 boost=-eta.BoostVector(); |
1034 | |
1035 | double eta_mass=eta.M(); |
1036 | eta.Boost(boost); |
1037 | |
1038 | one.Boost(boost); |
1039 | two.Boost(boost); |
1040 | three.Boost(boost); |
1041 | |
1042 | double m1=one.M(),m2=two.M(),m3=three.M(); |
1043 | double E1=one.E(),E2=two.E(),E3=three.E(); |
1044 | double T1=E1-m1; |
1045 | double T2=E2-m2; |
1046 | double T3=E3-m3; |
1047 | double Q_eta=eta_mass-m1-m2-m3; |
1048 | double X=sqrt(3.)*(T2-T3)/Q_eta; |
1049 | double Y=3.*T1/Q_eta-1.; |
1050 | thrown_dalitzXY->Fill(X,Y); |
1051 | |
1052 | double z_dalitz=X*X+Y*Y; |
1053 | |
1054 | thrown_dalitzZ->Fill(z_dalitz); |
1055 | } |
1056 | |
1057 | for (int j=0;j<num_decay_particles;j++){ |
1058 | if (particle_types[j]!=Unknown) { |
1059 | output_particle_vectors.push_back(*phase_space.GetDecay(j)); |
1060 | output_particle_types.push_back(particle_types[j]); |
1061 | } else { |
1062 | TGenPhaseSpace phase_space2; |
1063 | phase_space2.SetDecay(*phase_space.GetDecay(j),num_res_decay_particles, |
1064 | res_decay_masses); |
1065 | weight=0.,rand_weight=1.; |
1066 | do{ |
1067 | weight=phase_space2.Generate(); |
1068 | rand_weight=myrand->Uniform(1.); |
1069 | } while (rand_weight>weight); |
1070 | for (unsigned int im=0;im<num_res_decay_particles;im++){ |
1071 | output_particle_types.push_back(res_particle_types[im]); |
1072 | output_particle_vectors.push_back(*phase_space2.GetDecay(im)); |
1073 | } |
1074 | } |
1075 | } |
1076 | #ifdef HAVE_EVTGEN1 |
1077 | } |
1078 | #endif //HAVE_EVTGEN |
1079 | |
1080 | |
1081 | WriteEvent(i,beam,vert,output_particle_types,output_particle_vectors,output_particle_decays,secondary_vertices,file); |
1082 | |
1083 | if (((10*i)%Nevents)==0) cout << 100.*double(i)/double(Nevents) << "\% done" << endl; |
1084 | } |
1085 | |
1086 | |
1087 | |
1088 | rootfile->Write(); |
1089 | rootfile->Close(); |
1090 | |
1091 | |
1092 | close_s_HDDM(file); |
1093 | cout<<endl<<"Closed HDDM file"<<endl; |
1094 | cout<<" "<<Nevents<<" event written to "<<output_file_name<<endl; |
1095 | |
1096 | |
1097 | |
1098 | if (res_decay_masses!=NULL__null) delete []res_decay_masses; |
1099 | |
1100 | return 0; |
1101 | } |