Bug Summary

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