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