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