Bug Summary

File:programs/Simulation/genEtaRegge/genEtaRegge.cc
Location:line 1056, column 5
Description:Value stored to 'weight' is never read

Annotated Source Code

1// Main program for generating eta events.
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>
22using 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
47typedef struct {
48 bool decayed = false; // not used for now
49 int parent_id;
50 vector<Particle_t> ids;
51 vector<TLorentzVector> p4vs;
52} secondary_decay_t;
53
54// trim from start (in place)
55static 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// trim from end (in place)
62static 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// trim from both ends (in place)
69static inline void trim(std::string &s) {
70 ltrim(s);
71 rtrim(s);
72}
73
74// Masses
75const double m_p=0.93827; // GeV
76const double m_p_sq=m_p*m_p;
77double m_eta=0.54775; // GeV
78double m_eta_sq=m_eta*m_eta;
79// Width
80double width=0.;
81// Coupling constants
82double g_rho_eta_gamma=0.81;
83double g_omega_eta_gamma=0.29;
84double g_eta_gamma_gamma=0.0429;
85double g_phi_eta_gamma=0.38;
86
87int Nevents=10000;
88int runNo=10000;
89bool debug=false;
90
91// Diagnostic histograms
92TH1D *thrown_t;
93TH1D *thrown_dalitzZ;
94TH1D *thrown_Egamma;
95TH2D *thrown_dalitzXY;
96TH2D *thrown_theta_vs_p;
97TH2D *thrown_theta_vs_p_eta;
98TH1D *cobrems_vs_E;
99
100char input_file_name[250]="eta548.in";
101char output_file_name[250]="eta_gen.hddm";
102
103// Non-default option to generate uniform t-distribution from tmin to tmax
104/// (calculating cross section at fixed t_uniform_eval)
105bool gen_uniform_t=false;
106double t_uniform_eval=-1; // Only used if gen_uniform is true. Currently,
107float t_min_uniform=0; // takes min(t_min_uniform,t_0) so as to avoid unphysical t_values
108float t_max_uniform=-3; // takes max(t_max_uniform,t_max) so as to avoid unphysical t_values
109
110void 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// ParseCommandLineArguments
130//-----------
131void 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 // For command line options -tmin and -tmax
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); // Make sure value is negative, no matter what was supplied
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); // Make sure value is negative, no matter what was supplied
154 }
155 }
156
157 // For all other command line options (which are exactly one character long)
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// Cross section dsigma/dt derived from Laget(2005)
190double CrossSection(double s,double t,double p_gamma,double p_eta,double theta){
191 // Coupling constants
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 //double kappa_phi=0.;
200
201 // Angular quantities
202 double sintheta=sin(theta);
203 double costheta=cos(theta);
204
205 // Kinematic quantities for exchange particle
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 // Kinematic quantities for target
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 // Kinematic quantities for beam
218 double p_gamma_sq=p_gamma*p_gamma;
219 // other kinematic quantities
220 double pt_dot_q=pt0*q0-pt3*q3;
221
222 // Mass scale for Regge propagators
223 double s0=1.0;
224
225 // Regge trajectory for rho
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 // Regge trajectory for omega
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 // Regge trajectory for phi(1020)
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 // printf("> 36: %f\n",p_gamma*p_eta*sqrt(mass_factor));
241
242 // amplitude factors for terms involving 0, 1 and 2 powers of kappa
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 // rho amplitude
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 // omega amplitude
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 // phi amplitude
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 // Primakoff amplitude
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 // M_primakoff=0.;
279 //M_primakoff_sq=0.;
280
281 //M_omega=0.;
282 // M_omega_sq=0.;
283
284 //M_rho=0.;
285 // M_rho_sq=0.;
286
287 //M_phi_sq=0.;
288 //M_phi=0.;
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.; // Convert to micro-barns
303 double dsigma_dt=hbarc_sq*M_sq/(4.*64.*M_PI3.14159265358979323846*s*p_gamma_sq);
304 // the extra factor for is for 2 photon spins x 2 proton spins
305
306 return(dsigma_dt);
307
308
309}
310
311// Put particle data into hddm format and output to file
312void 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 // Beam
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 // Target
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 // Primary vertex
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 // Final state particles
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; // zero out particle type info so that hdgeant won't decay the particle. maybe there is a better way?
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++; /* unique value for this particle within the event */
375 ps->in[ps->mult].parentid = 0; /* All internally generated particles have no parent */
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 // write out any secondary vertices (like pi0 decays)
384 // vector< secondary_decay_t > &secondary_vertices,
385 if(secondary_vertices.size() > 0) {
386 int vertex_ind = 1;
387 for(auto& the_vertex : secondary_vertices) {
388 // assume that all particles are generated at the same vertex
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 // add in the particles associated with this vertex
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++; /* unique value for this particle within the event */
404 ps->in[ps->mult].parentid = the_vertex.parent_id; /* All internally generated particles have no parent */
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++; // get ready for next iteration
414 }
415 }
416 flush_s_HDDM(thisOutputEvent,file);
417}
418
419// Create some diagnostic histograms
420void 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// Create a graph of the cross section dsigma/dt as a function of -t
446void GraphCrossSection(double &xsec_max){
447 // beam energy in lab
448 double Egamma=cobrems_vs_E->GetBinLowEdge(1); // get from CobremsGenerator histogram
449
450 // CM energy
451 double s=m_p*(m_p+2.*Egamma);
452 double Ecm=sqrt(s);
453
454 // Momenta of incoming photon and outgoing eta and proton in cm frame
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 // Momentum transfer t
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// main
493//-----------
494int main(int narg, char *argv[])
495{
496 ParseCommandLineArguments(narg, argv);
497
498
499 // open ROOT file
500 string rootfilename="eta_gen.root";
501 TFile *rootfile=new TFile(rootfilename.c_str(),"RECREATE",
502 "Produced by genEta");
503
504 // open HDDM file
505 s_iostream_t *file = init_s_HDDM(output_file_name);
506
507
508 // Initialize random number generator
509 TRandom3 *myrand=new TRandom3(0);// If seed is 0, the seed is automatically computed via a TUUID object, according to TRandom3 documentation
510
511 // Fixed target
512 TLorentzVector target(0.,0.,0.,m_p);
513
514 //----------------------------------------------------------------------------
515 // Get production (Egamma range) and decay parameters from input file
516 //----------------------------------------------------------------------------
517
518 // Start reading the input file
519 ifstream infile(input_file_name);
520 if (!infile.is_open()){
521 cerr << "Input file missing! Exiting..." <<endl;
522 exit(-1);
523 }
524
525 // Get beam properties configuration file
526 string comment_line;
527 getline(infile,comment_line);
528 string beamConfigFile;
529 infile >> beamConfigFile;
530 infile.ignore(); // ignore the '\n' at the end of this line
531
532 cout << "Photon beam configuration file " << beamConfigFile.data() << endl;
533
534 // Get decaying particle mass and width
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(); // ignore the '\n' at the end of this line
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 // Get coupling constants for photon vertex
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(); // ignore the '\n' at the end of this line
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 // Get number of decay particles
562 int num_decay_particles=0;
563 getline(infile,comment_line);
564 infile >> num_decay_particles;
565 infile.ignore(); // ignore the '\n' at the end of this line
566
567 cout << "number of decay particles = " << num_decay_particles << endl;
568
569 if( abs(t_max_uniform)>21. ) { // Max t at GlueX endpoint energy is about 20 GeV^2. Reset value to protect against inefficient accept/reject.
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 // check to see if we should use EvtGen
579 EvtParticle* parent(0);
580 EvtGen *myGenerator = nullptr;
581 EvtId EtaId; // read this in from file below - this is actually a class, which gets filled with reasonable defaults
582
583 // if we don't explicitly define the decay particles in the config file,
584 // then use EvtGen to generate the decays
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 // get the produced particle type
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 // initialize EvtGen
599 const char* evtgen_home_env_ptr = std::getenv("EVTGENDIR");
600 string EVTGEN_HOME = (evtgen_home_env_ptr==nullptr) ? "." : evtgen_home_env_ptr; // default to the current directory
601
602 // Define the random number generator
603 EvtRandomEngine* eng = 0;
604#ifdef EVTGEN_CPP11
605 // Use the Mersenne-Twister generator (C++11 only)
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 //Initialize the generator - read in the decay table and particle properties
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 // Make special GlueX definitions
633 GlueX_EvtGen::RegisterGlueXModels();
634
635 // open optional user decay file, if it exists
636 struct stat buffer;
637 if(stat("userDecay.dec", &buffer) == 0)
638 myGenerator->readUDecay("userDecay.dec");
639
640 // Now that we have initialized EvtGen, we can access things like the particle DB
641 EtaId = EvtPDL::getId(std::string(particle_type));
642
643 }
644#endif // HAVE_EVTGEN
645
646 // Set up vectors of particle ids
647 vector<Particle_t>particle_types;
648 //double *decay_masses =new double[num_decay_particles];
649 vector<double> decay_masses(num_decay_particles);
650 double *res_decay_masses=NULL__null;
651 vector<Particle_t>res_particle_types;
652
653 // GEANT ids of decay particles
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(); // ignore the '\n' at the end of this line
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 // Create some diagonistic histographs
706 CreateHistograms(beamConfigFile);
707
708 // Make a TGraph of the cross section at a fixed beam energy
709 double xsec_max=0.;
710 GraphCrossSection(xsec_max);
711
712 //----------------------------------------------------------------------------
713 // Event generation loop
714 //----------------------------------------------------------------------------
715 for (int i=1;i<=Nevents;i++){
716 double Egamma=0.;
717 double xsec=0.,xsec_test=0.;
718
719 // Polar angle in center of mass frame
720 double theta_cm=0.;
721
722 // Eta momentum in cm
723 double p_eta=0.;
724
725 // Transfer 4-momentum;
726 double t=0.;
727
728 // vertex position at target
729 float vert[4]={0.,0.,0.,0.};
730
731 // use the rejection method to produce eta's based on the cross section
732 do{
733 // First generate a beam photon using bremsstrahlung spectrum
734 Egamma = cobrems_vs_E->GetRandom();
735
736 // CM energy
737 double s=m_p*(m_p+2.*Egamma);
738 double Ecm=sqrt(s);
739
740 // Momenta of incoming photon and outgoing eta and proton in cm frame
741 double p_gamma=(s-m_p_sq)/(2.*Ecm);
742
743 // Generate mass distribution for unstable particle in the final state
744 // with non-negligible width if specified in the input file
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.; // Meson radius in GeV^-1: 5 GeV^-1 -> 1 fm
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 // double BW_at_m_min=0.;
765 //double BW_at_m_max=0.;
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 // Take into account width of resonance, but apply a practical minimum
805 // for the width, overwise we are just wasting cpu cycles...
806 // Use a relativistic Breit-Wigner distribution for the shape.
807 double m_max_=m_p*(sqrt(1.+2.*Egamma/m_p)-1.);
808 double m_min_=decay_masses[0]; // will add the second mass below
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 // Define an effective mass in an ad hoc way: we assume that in the
820 // CM one particle goes in one direction and the two other particles
821 // go in the opposite direction such that p1=-p2-p3. The effective
822 // mass of the 2-3 system must be something between min=m2+m3
823 // and max=M-m1, where M is the mass of the resonance. For
824 // simplicity use the average of these two extremes.
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 // Check that the decay products are consistent with a particle of mass
858 // m_eta...
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 // if using EvtGen, we don't know what the decay products are a priori
866 // use a non-relativistic BW instead
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 // use a cutoff of +- 5 times the width so we don't populate the far tails too much
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 // Momentum transfer t
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 // Generate cos(theta) with a uniform distribution and compute the cross
891 // section at this value
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 // If generating a sample uniform in t, we need to fix t and re-calculate theta_cm based on it. Others do not depend on t.
900 if(gen_uniform_t) {
901 //Cross section at fixed t value (t_tmp)
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 // Make t uniform, calculate theta_cm based off of it
906 t=myrand->Uniform(t_max_uniform, min( float(t0),t_min_uniform) ); // If t_min_uniform provided is unphysical, then use physical t_min.
907 theta_cm=2.*asin(0.5*sqrt( (t0-t)/(p_gamma*p_eta) ) );
908 if( std::isnan(theta_cm)==true ) xsec=-1.; // Lazy person's way of skipping unphysical theta_cm. Breaking do/while to accept event will never be satisfied for this case.
909 }
910
911 // Generate a test value for the cross section
912 xsec_test=myrand->Uniform(xsec_max);
913 }
914 while (xsec_test>xsec);
915
916 // Generate phi using uniform distribution
917 double phi_cm=myrand->Uniform(2.*M_PI3.14159265358979323846);
918
919 // beam 4-vector (ignoring px and py, which are extremely small)
920 TLorentzVector beam(0.,0.,Egamma,Egamma);
921 thrown_Egamma->Fill(Egamma);
922
923 // Velocity of the cm frame with respect to the lab frame
924 TVector3 v_cm=(1./(Egamma+m_p))*beam.Vect();
925 // Four-moementum of the eta in the CM frame
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 //Boost the eta 4-momentum into the lab
931 eta4.Boost(v_cm);
932
933 // Compute the 4-momentum for the recoil proton
934 TLorentzVector proton4=beam+target-eta4;
935
936 //proton4.Print();
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 // Other diagnostic histograms
941 thrown_t->Fill(-t);
942
943 // Gather the particles in the reaction and write out event in hddm format
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 // Set up the parent particle
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 // just generate the eta and let the decay be done by an external program like decay_evtgen
962 // this plays better with MCWrapper, apparently...
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 // Generate the event
969 myGenerator->generateDecay(parent);
970
971 // Write out resulting particles
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 // see if any of the particles decay and add info on them
981 // should be mostly pi0's, but we should go recursive...
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 { // no evtgen
1004#endif //HAVE_EVTGEN
1005 // Generate 3-body decay of eta according to phase space
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 // Histograms of Dalitz distribution
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; // pi0 for charged channel
1036 double T2=E2-m2; // pi+ for charged channel
1037 double T3=E3-m3; // pi- for charged channel
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 //printf("z %f\n",z_dalitz);
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 // Write Event to HDDM file
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 // Write histograms and close root file
1079 rootfile->Write();
1080 rootfile->Close();
1081
1082 // Close HDDM file
1083 close_s_HDDM(file);
1084 cout<<endl<<"Closed HDDM file"<<endl;
1085 cout<<" "<<Nevents<<" event written to "<<output_file_name<<endl;
1086
1087 // Cleanup
1088 //delete []decay_masses;
1089 if (res_decay_masses!=NULL__null) delete []res_decay_masses;
1090
1091 return 0;
1092}