Bug Summary

File:programs/Simulation/genEtaRegge/genEtaRegge.cc
Location:line 1010, column 7
Description:Value stored to 't' 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/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
48typedef struct {
49 bool decayed = false; // not used for now
50 int parent_id;
51 vector<Particle_t> ids;
52 vector<TLorentzVector> p4vs;
53} secondary_decay_t;
54
55// trim from start (in place)
56static 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// trim from end (in place)
63static 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// trim from both ends (in place)
70static inline void trim(std::string &s) {
71 ltrim(s);
72 rtrim(s);
73}
74
75//IA
76TString m_str_Nucleus = "";
77TString m_str_Participant = "";
78TString m_str_Spectator = "";
79TH1F * m_h_PFermi;
80double m_mass_nuclei = 0;
81double m_ParticipantMass = 0;
82double m_SpectatorMass = 0;
83
84// Masses
85double m_p=0.93827; // GeV
86double m_p_sq=m_p*m_p;
87double m_eta=0.54775; // GeV
88double m_eta_sq=m_eta*m_eta;
89// Width
90double width=0.;
91// Coupling constants
92double g_rho_eta_gamma=0.81;
93double g_omega_eta_gamma=0.29;
94double g_eta_gamma_gamma=0.0429;
95double g_phi_eta_gamma=0.38;
96
97int Nevents=10000;
98int runNo=10000;
99bool debug=false;
100
101// Diagnostic histograms
102TH1D *thrown_t;
103TH1D *thrown_dalitzZ;
104TH1D *thrown_Egamma;
105TH2D *thrown_dalitzXY;
106TH2D *thrown_theta_vs_p;
107TH2D *thrown_theta_vs_p_eta;
108TH1D *cobrems_vs_E;
109TH1D *thrown_FermiP;
110TH1D *thrown_f;
111
112char input_file_name[250]="eta548.in";
113char output_file_name[250]="eta_gen.hddm";
114
115// Non-default option to generate uniform t-distribution from tmin to tmax
116/// (fixed to cross section at tflat_min)
117bool gen_uniform_t=false;
118float tflat_min=100.; // Physical values are negative. Cross section at larger |t| is equal to cross section at this number (for fixed E_gamma)
119float tflat_max=100.; // Physical values are negative
120
121void 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// ParseCommandLineArguments
138//-----------
139void 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// Cross section dsigma/dt derived from Laget(2005)
176double CrossSection(double s,double t,double p_gamma,double p_eta,double theta){
177 // Coupling constants
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 //double kappa_phi=0.;
186
187 // Angular quantities
188 double sintheta=sin(theta);
189 double costheta=cos(theta);
190
191 // Kinematic quantities for exchange particle
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 // Kinematic quantities for target
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 // Kinematic quantities for beam
204 double p_gamma_sq=p_gamma*p_gamma;
205 // other kinematic quantities
206 double pt_dot_q=pt0*q0-pt3*q3;
207
208 // Mass scale for Regge propagators
209 double s0=1.0;
210
211 // Regge trajectory for rho
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 // Regge trajectory for omega
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 // Regge trajectory for phi(1020)
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 // printf("> 36: %f\n",p_gamma*p_eta*sqrt(mass_factor));
227
228 // amplitude factors for terms involving 0, 1 and 2 powers of kappa
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 // rho amplitude
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 // omega amplitude
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 // phi amplitude
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 // Primakoff amplitude
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 // M_primakoff=0.;
265 //M_primakoff_sq=0.;
266
267 //M_omega=0.;
268 // M_omega_sq=0.;
269
270 //M_rho=0.;
271 // M_rho_sq=0.;
272
273 //M_phi_sq=0.;
274 //M_phi=0.;
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.; // Convert to micro-barns
289 double dsigma_dt=hbarc_sq*M_sq/(4.*64.*M_PI3.14159265358979323846*s*p_gamma_sq);
290 // the extra factor for is for 2 photon spins x 2 proton spins
291
292 return(dsigma_dt);
293
294
295}
296
297// Put particle data into hddm format and output to file
298void 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 // Beam
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 // Target
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 //ta->type = Proton;
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 //ta->momentum->px = 0.;
343 //ta->momentum->py = 0.;
344 //ta->momentum->pz = 0.;
345 //ta->momentum->E = ParticleMass(ta->type);
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 // Primary vertex
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 // Final state particles
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; // zero out particle type info so that hdgeant won't decay the particle. maybe there is a better way?
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++; /* unique value for this particle within the event */
378 ps->in[ps->mult].parentid = 0; /* 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 = 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 // write out any secondary vertices (like pi0 decays)
387 // vector< secondary_decay_t > &secondary_vertices,
388 if(secondary_vertices.size() > 0) {
389 int vertex_ind = 1;
390 for(auto& the_vertex : secondary_vertices) {
391 // assume that all particles are generated at the same vertex
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 // add in the particles associated with this vertex
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++; /* unique value for this particle within the event */
407 ps->in[ps->mult].parentid = the_vertex.parent_id; /* All internally generated particles have no parent */
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++; // get ready for next iteration
417 }
418 }
419 flush_s_HDDM(thisOutputEvent,file);
420}
421
422// Create some diagnostic histograms
423void 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// Create a graph of the cross section dsigma/dt as a function of -t
454void GraphCrossSection(double &xsec_max){
455 // beam energy in lab
456 double Egamma=cobrems_vs_E->GetBinLowEdge(1); // get from CobremsGenerator histogram
457
458 // CM energy
459 double s=m_p*(m_p+2.*Egamma);
460 double Ecm=sqrt(s);
461
462 // Momenta of incoming photon and outgoing eta and proton in cm frame
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 // Momentum transfer t
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// main
501//-----------
502int main(int narg, char *argv[])
503{
504 ParseCommandLineArguments(narg, argv);
505
506
507 // open ROOT file
508 string rootfilename="eta_gen.root";
509 TFile *rootfile=new TFile(rootfilename.c_str(),"RECREATE",
510 "Produced by genEta");
511
512 // open HDDM file
513 s_iostream_t *file = init_s_HDDM(output_file_name);
514
515
516 // Initialize random number generator
517 TRandom3 *myrand=new TRandom3(0);// If seed is 0, the seed is automatically computed via a TUUID object, according to TRandom3 documentation
518
519 // Fixed target
520 TLorentzVector target(0.,0.,0.,m_p);
521
522 //----------------------------------------------------------------------------
523 // Get production (Egamma range) and decay parameters from input file
524 //----------------------------------------------------------------------------
525
526 // Start reading the input file
527 ifstream infile(input_file_name);
528 if (!infile.is_open()){
529 cerr << "Input file missing! Exiting..." <<endl;
530 exit(-1);
531 }
532
533 // IA, get generator config file
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;//003 001 3H
555 } else if (m_str_Participant == "Neutron") {
556 m_ParticipantMass = 0.93956;
557 m_SpectatorMass = 2.808391;//003 002 3He
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;//011 005 11B
564 } else if (m_str_Participant == "Neutron") {
565 m_ParticipantMass = 0.93956;
566 m_SpectatorMass = 10.254018;//011 006 11C
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 // Get beam properties configuration file
599 string comment_line;
600 getline(infile,comment_line);
601 string beamConfigFile;
602 infile >> beamConfigFile;
603 infile.ignore(); // ignore the '\n' at the end of this line
604
605 cout << "Photon beam configuration file " << beamConfigFile.data() << endl;
606
607 // Get decaying particle mass and width
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(); // ignore the '\n' at the end of this line
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 // Get coupling constants for photon vertex
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(); // ignore the '\n' at the end of this line
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 // Get number of decay particles
635 int num_decay_particles=0;
636 getline(infile,comment_line);
637 infile >> num_decay_particles;
638 infile.ignore(); // ignore the '\n' at the end of this line
639
640 cout << "number of decay particles = " << num_decay_particles << endl;
641
642
643
644 bool use_evtgen = false;
645#ifdef HAVE_EVTGEN1
646 // check to see if we should use EvtGen
647 EvtParticle* parent(0);
648 EvtGen *myGenerator = nullptr;
649 EvtId EtaId; // read this in from file below - this is actually a class, which gets filled with reasonable defaults
650
651 // if we don't explicitly define the decay particles in the config file,
652 // then use EvtGen to generate the decays
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 // get the produced particle type
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 // initialize EvtGen
667 const char* evtgen_home_env_ptr = std::getenv("EVTGENDIR");
668 string EVTGEN_HOME = (evtgen_home_env_ptr==nullptr) ? "." : evtgen_home_env_ptr; // default to the current directory
669
670 // Define the random number generator
671 EvtRandomEngine* eng = 0;
672#ifdef EVTGEN_CPP11
673 // Use the Mersenne-Twister generator (C++11 only)
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 //Initialize the generator - read in the decay table and particle properties
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 // Make special GlueX definitions
701 GlueX_EvtGen::RegisterGlueXModels();
702
703 // open optional user decay file, if it exists
704 struct stat buffer;
705 if(stat("userDecay.dec", &buffer) == 0)
706 myGenerator->readUDecay("userDecay.dec");
707
708 // Now that we have initialized EvtGen, we can access things like the particle DB
709 EtaId = EvtPDL::getId(std::string(particle_type));
710
711 }
712#endif // HAVE_EVTGEN
713
714 // Set up vectors of particle ids
715 vector<Particle_t>particle_types;
716 //double *decay_masses =new double[num_decay_particles];
717 vector<double> decay_masses(num_decay_particles);
718 double *res_decay_masses=NULL__null;
719 vector<Particle_t>res_particle_types;
720
721 // GEANT ids of decay particles
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(); // ignore the '\n' at the end of this line
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 // Search for lines in input file starting with "tflat_min" or "tflat_max", if found we reset globals
772 while( !infile.eof() ) {
773 string line = "";
774 string tflat_string = "";
775 getline(infile,line);
776 if(line.length() < 11) continue;
777 // Yes this code is ugly, but works. I miss python.
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 // Yes this code is ugly, but works. I miss python.
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 // Create some diagonistic histographs
809 CreateHistograms(beamConfigFile,num_decay_particles);
810
811 // Make a TGraph of the cross section at a fixed beam energy
812 double xsec_max=0.;
813 GraphCrossSection(xsec_max);
814
815 //----------------------------------------------------------------------------
816 // Event generation loop
817 //----------------------------------------------------------------------------
818 for (int i=1;i<=Nevents;i++){
819 double Egamma=0.;
820 double xsec=0.,xsec_test=0.;
821
822 // Polar angle in center of mass frame
823 double theta_cm=0.;
824
825 // Eta momentum in cm
826 double p_eta=0.;
827
828 // Transfer 4-momentum;
829 double t=0.;
830
831 // vertex position at target
832 float vert[4]={0.,0.,0.,0.};
833
834 // IA variables
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 // use the rejection method to produce eta's based on the cross section
847 do{
848 // First generate a beam photon using bremsstrahlung spectrum
849 Egamma = cobrems_vs_E->GetRandom();
850
851 // CM energy
852 double s=m_p*(m_p+2.*Egamma);
853 double Ecm=sqrt(s);
854
855 // IA, momenta of incoming photon and outgoing eta and proton in cm frame
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 // Momenta of incoming photon and outgoing eta and proton in cm frame
863 double p_gamma=(s-m_p_sq)/(2.*Ecm);
864
865 // Generate mass distribution for unstable particle in the final state
866 // with non-negligible width if specified in the input file
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.; // Meson radius in GeV^-1: 5 GeV^-1 -> 1 fm
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 // double BW_at_m_min=0.;
887 //double BW_at_m_max=0.;
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 // Take into account width of resonance, but apply a practical minimum
927 // for the width, overwise we are just wasting cpu cycles...
928 // Use a relativistic Breit-Wigner distribution for the shape.
929 double m_max_=m_p*(sqrt(1.+2.*Egamma/m_p)-1.);
930 double m_min_=decay_masses[0]; // will add the second mass below
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 // Define an effective mass in an ad hoc way: we assume that in the
942 // CM one particle goes in one direction and the two other particles
943 // go in the opposite direction such that p1=-p2-p3. The effective
944 // mass of the 2-3 system must be something between min=m2+m3
945 // and max=M-m1, where M is the mass of the resonance. For
946 // simplicity use the average of these two extremes.
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 // Check that the decay products are consistent with a particle of mass
980 // m_eta...
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 // if using EvtGen, we don't know what the decay products are a priori
988 // use a non-relativistic BW instead
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 // use a cutoff of +- 5 times the width so we don't populate the far tails too much
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 // Momentum transfer t
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 // Generate cos(theta) with a uniform distribution and compute the cross
1013 // section at this value
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 // 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.
1022 if(gen_uniform_t&&t<tflat_min&&t>tflat_max) {
1023 //Cross section at fixed t value (tflat_min)
1024 double t_tmp = tflat_min;
1025 // if(t0<t_tmp && t_tmp < 0. ) t_tmp=t0-0.00001; //If tflat_min is unphysically small, use (essentially) t0 instead
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 // Make t uniform, calculate theta_cm based off of it
1029 t=myrand->Uniform(tflat_max, min( float(t0),tflat_min) ); // If t_min_uniform provided is unphysical, then use physical t_min.
1030 theta_cm=2.*asin(0.5*sqrt( (t0-t)/(p_gamma*p_eta) ) );
1031 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.
1032 }
1033
1034 // Generate a test value for the cross section
1035 xsec_test=myrand->Uniform(xsec_max);
1036 }
1037 while (xsec_test>xsec);
1038
1039 // Generate phi using uniform distribution
1040 double phi_cm=myrand->Uniform(2.*M_PI3.14159265358979323846);
1041
1042 // beam 4-vector (ignoring px and py, which are extremely small)
1043 TLorentzVector beam(0.,0.,Egamma,Egamma);
1044 thrown_Egamma->Fill(Egamma);
1045
1046 // IA nucleon/nuclei target
1047 if (m_str_Nucleus != "") {
1048 target = TLorentzVector(p_Fermi_x, p_Fermi_y, p_Fermi_z, ParticipantEnergy);
1049 //cout <<"rewrite target p4 w/ fermi"<<endl;
1050 } else if (m_str_Nucleus == "" && m_str_Participant != "") {
1051 target = TLorentzVector(0, 0, 0, m_p);
1052 }
1053
1054 // Velocity of the cm frame with respect to the lab frame
1055 TVector3 v_cm=(1./(Egamma+m_p))*beam.Vect();
1056 // Four-moementum of the eta in the CM frame
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 //Boost the eta 4-momentum into the lab
1062 //eta4.Boost(v_cm);
1063 // IA modified boost
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 // Compute the 4-momentum for the recoil proton
1072 TLorentzVector proton4=beam+target-eta4;
1073
1074 //proton4.Print();
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 // Other diagnostic histograms
1079 thrown_t->Fill(-t);
1080
1081 // Gather the particles in the reaction and write out event in hddm format
1082 vector<TLorentzVector>output_particle_vectors;
1083 output_particle_vectors.push_back(proton4);
1084
1085 vector<Particle_t>output_particle_types;
1086 //output_particle_types.push_back(Proton);
1087 // IA modifications
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 // Set up the parent particle
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 // just generate the eta and let the decay be done by an external program like decay_evtgen
1107 // this plays better with MCWrapper, apparently...
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 // Generate the event
1114 myGenerator->generateDecay(parent);
1115
1116 // Write out resulting particles
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 // see if any of the particles decay and add info on them
1126 // should be mostly pi0's, but we should go recursive...
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 { // no evtgen
1149#endif //HAVE_EVTGEN
1150 // Generate 3-body decay of eta according to phase space
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 // Histograms of Dalitz distribution
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; // pi0 for charged channel
1181 double T2=E2-m2; // pi+ for charged channel
1182 double T3=E3-m3; // pi- for charged channel
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 //printf("z %f\n",z_dalitz);
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 // Write Event to HDDM file
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 // Write histograms and close root file
1224 rootfile->Write();
1225 rootfile->Close();
1226
1227 // Close HDDM file
1228 close_s_HDDM(file);
1229 cout<<endl<<"Closed HDDM file"<<endl;
1230 cout<<" "<<Nevents<<" event written to "<<output_file_name<<endl;
1231
1232 // Cleanup
1233 //delete []decay_masses;
1234 if (res_decay_masses!=NULL__null) delete []res_decay_masses;
1235
1236 return 0;
1237}