From 6cb9108ab6ac79eef05572ef3f82d66467262806 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Fri, 18 Sep 2026 11:08:34 -0400 Subject: [PATCH 1/3] Add amplitude for eta->e+e-pi+pi- decay --- .../Simulation/genEtaRegge/genEtaRegge.cc | 104 ++++++++++++++++-- 1 file changed, 96 insertions(+), 8 deletions(-) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index 358e6aeba..87c563948 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -93,6 +93,9 @@ double g_rho_eta_gamma=0.81; double g_omega_eta_gamma=0.29; double g_eta_gamma_gamma=0.0429; double g_phi_eta_gamma=0.38; +double G=0.; // CV coupling constant + +double max_asq=0.; int Nevents=10000; int runNo=10000; @@ -108,6 +111,7 @@ TH2D *thrown_theta_vs_p_eta; TH1D *cobrems_vs_E; TH1D *thrown_FermiP; TH1D *thrown_f; +TH1D *thrown_qsq; char input_file_name[250]="eta548.in"; char output_file_name[250]="eta_gen.hddm"; @@ -147,6 +151,10 @@ void ParseCommandLineArguments(int narg, char* argv[]) if(ptr[0] == '-'){ switch(ptr[1]){ case 'h': Usage(); break; + case 'G': + sscanf(&ptr[2],"%lf",&G); + cout << "G=" << G << endl; + break; case 'I': sscanf(&ptr[2],"%s",input_file_name); break; @@ -445,6 +453,9 @@ void CreateHistograms(string beamConfigFile,int num_decay_particles){ thrown_dalitzZ=new TH1D("thrown_dalitzZ","thrown dalitz Z",110,-0.05,1.05); thrown_dalitzXY=new TH2D("thrown_dalitzXY","Dalitz distribution Y vs X",100,-1.,1.,100,-1.,1); } + if (num_decay_particles==4){ + thrown_qsq=new TH1D("thrown_qsq",";q^{2} [Gev^{2}]",800,0,0.08); + } BeamProperties beamProp(beamConfigFile); cobrems_vs_E = (TH1D*)beamProp.GetFlux(); @@ -496,7 +507,78 @@ void GraphCrossSection(double &xsec_max){ << " micro-barns"<1e-16){ + TVector3 qdir=qvec.Unit(); + TVector3 pipvec=piplus.Vect(); + kvec.RotateUz(qdir); + pipvec.RotateUz(qdir); + + double k_dot_pm=k.Dot(piminus); + double k_dot_pp=k.Dot(piplus); + double q_dot_pm=q4.Dot(piminus); + double q_dot_pp=q4.Dot(piplus); + + amp_sq=2*efac*efac/(q_sq*q_sq)*weight + *(M*M*m_eta_sq*qvec.Mag2() + *(q_sq*pipvec.Perp2()-pow(kvec.x()*pipvec.y()-kvec.y()*pipvec.x(),2)) + +pow(Escale*(q_sq+2*q_dot_pp),2) + *(q_dot_pm*q_dot_pm-k_dot_pm*k_dot_pm-q_sq*pion_mass_sq) + +pow(Escale*(q_sq+2*q_dot_pm),2) + *(q_dot_pp*q_dot_pp-k_dot_pp*k_dot_pp-q_sq*pion_mass_sq) + -2.*Escale*Escale*(q_sq+2*q_dot_pm)*(q_sq+2*q_dot_pp) + *(q_dot_pp*q_dot_pm-q_sq*piplus.Dot(piminus)-k_dot_pp*k_dot_pm) + -2.*M*Escale*((q_sq+2*q_dot_pm)*k_dot_pp-(q_sq+2*q_dot_pp)*k_dot_pm) + *m_eta*qvec.Mag()*(kvec.x()*pipvec.y()-kvec.y()*pipvec.x()) + ); + if (amp_sq>0.){ + //thrown_qsq->Fill(q_sq,amp_sq*weight); + if (amp_sq>max_asq) max_asq=amp_sq; + } + } + rand_amp_sq=myrand->Uniform(15.0); + } while (rand_amp_sq>amp_sq); + if (q_sq>0.){ + thrown_qsq->Fill(q_sq); + } +} + //----------- // main //----------- @@ -1058,7 +1140,7 @@ int main(int narg, char *argv[]) double pt=p_eta*sin(theta_cm); TLorentzVector eta4(pt*cos(phi_cm),pt*sin(phi_cm),p_eta*cos(theta_cm), sqrt(p_eta*p_eta+m_eta_sq)); - + //Boost the eta 4-momentum into the lab //eta4.Boost(v_cm); // IA modified boost @@ -1148,16 +1230,21 @@ int main(int narg, char *argv[]) } } else { // no evtgen #endif //HAVE_EVTGEN - // Generate 3-body decay of eta according to phase space + // Generate n-body decay of eta according to phase space TGenPhaseSpace phase_space; phase_space.SetDecay(eta4,num_decay_particles,decay_masses.data()); double weight=0.,rand_weight=1.; - do{ - weight=phase_space.Generate(); - rand_weight=myrand->Uniform(1.); + if (num_decay_particles<4){ + do{ + weight=phase_space.Generate(); + rand_weight=myrand->Uniform(1.); + } + while (rand_weight>weight); } - while (rand_weight>weight); - + else { + GenerateEpEmPipPim(myrand,phase_space); + } + // Histograms of Dalitz distribution if (num_decay_particles==3){ TLorentzVector one=*phase_space.GetDecay(0); @@ -1222,6 +1309,7 @@ int main(int narg, char *argv[]) if (((10*i)%Nevents)==0) cout << 100.*double(i)/double(Nevents) << "\% done" << endl; } + cout << "max " << max_asq << endl; // Write histograms and close root file rootfile->Write(); From aadc3d609d9d5d99778a6a1f6e58c33009e5be76 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Fri, 18 Sep 2026 14:15:04 -0400 Subject: [PATCH 2/3] Add some comments and clean up code a little --- .../Simulation/genEtaRegge/genEtaRegge.cc | 28 +++++++++++-------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index 87c563948..36d84fb8e 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -93,9 +93,9 @@ double g_rho_eta_gamma=0.81; double g_omega_eta_gamma=0.29; double g_eta_gamma_gamma=0.0429; double g_phi_eta_gamma=0.38; -double G=0.; // CV coupling constant +double G=0.; // CP-violation parameter for eta->e+e-pi+pi- -double max_asq=0.; +//double max_asq=0.; int Nevents=10000; int runNo=10000; @@ -507,6 +507,8 @@ void GraphCrossSection(double &xsec_max){ << " micro-barns"<e+e-pi+pi- decay using the model +// described in Gao hep-ph/0202002 void GenerateEpEmPipPim(TRandom3 *myrand,TGenPhaseSpace &phase_space){ double amp_sq=0., rand_amp_sq=0.; double q_sq=0.; @@ -517,6 +519,8 @@ void GenerateEpEmPipPim(TRandom3 *myrand,TGenPhaseSpace &phase_space){ double f0=1.1*fpi; double efac=sqrt(4.*M_PI/137.); double theta_mix=-20.*M_PI/180.; // eta-eta' mixing angle + // Scale factor for electromagnetic transition. Form is taken from eq. 16 + // and eq. 17 in Gao double Escale=efac*0.19*G/pow(m_eta,3); do { double weight=phase_space.Generate(); @@ -532,7 +536,8 @@ void GenerateEpEmPipPim(TRandom3 *myrand,TGenPhaseSpace &phase_space){ double s=pippim.M2(); TLorentzVector k=positron-electron; q_sq=q4.M2(); - + + // Magnetic transition term double M=efac/(8.*M_PI*M_PI*fpi*fpi) *(cos(theta_mix)/(sqrt(3.)*f8)-sqrt(2.)*sin(theta_mix)/(sqrt(3.)*f0)) *(1.-3.*mV2/(mV2-s)); @@ -554,23 +559,24 @@ void GenerateEpEmPipPim(TRandom3 *myrand,TGenPhaseSpace &phase_space){ double k_dot_pp=k.Dot(piplus); double q_dot_pm=q4.Dot(piminus); double q_dot_pp=q4.Dot(piplus); + double kx_py_minus_ky_px=kvec.x()*pipvec.y()-kvec.y()*pipvec.x(); amp_sq=2*efac*efac/(q_sq*q_sq)*weight *(M*M*m_eta_sq*qvec.Mag2() - *(q_sq*pipvec.Perp2()-pow(kvec.x()*pipvec.y()-kvec.y()*pipvec.x(),2)) + *(q_sq*pipvec.Perp2()-pow(kx_py_minus_ky_px,2)) +pow(Escale*(q_sq+2*q_dot_pp),2) *(q_dot_pm*q_dot_pm-k_dot_pm*k_dot_pm-q_sq*pion_mass_sq) +pow(Escale*(q_sq+2*q_dot_pm),2) *(q_dot_pp*q_dot_pp-k_dot_pp*k_dot_pp-q_sq*pion_mass_sq) -2.*Escale*Escale*(q_sq+2*q_dot_pm)*(q_sq+2*q_dot_pp) *(q_dot_pp*q_dot_pm-q_sq*piplus.Dot(piminus)-k_dot_pp*k_dot_pm) - -2.*M*Escale*((q_sq+2*q_dot_pm)*k_dot_pp-(q_sq+2*q_dot_pp)*k_dot_pm) - *m_eta*qvec.Mag()*(kvec.x()*pipvec.y()-kvec.y()*pipvec.x()) + -2.*M*Escale*((q_sq+2*q_dot_pm)*k_dot_pp-(q_sq+2*q_dot_pp)*k_dot_pm) + *m_eta*qvec.Mag()*kx_py_minus_ky_px ); - if (amp_sq>0.){ - //thrown_qsq->Fill(q_sq,amp_sq*weight); - if (amp_sq>max_asq) max_asq=amp_sq; - } + //if (amp_sq>0.){ + //thrown_qsq->Fill(q_sq,amp_sq*weight); + //if (amp_sq>max_asq) max_asq=amp_sq; + //} } rand_amp_sq=myrand->Uniform(15.0); } while (rand_amp_sq>amp_sq); @@ -1309,7 +1315,7 @@ int main(int narg, char *argv[]) if (((10*i)%Nevents)==0) cout << 100.*double(i)/double(Nevents) << "\% done" << endl; } - cout << "max " << max_asq << endl; + // cout << "max " << max_asq << endl; // Write histograms and close root file rootfile->Write(); From c5de2ed9e7d2852e7bee9e2255736a6e149036d0 Mon Sep 17 00:00:00 2001 From: staylorjlab Date: Mon, 21 Sep 2026 07:55:01 -0400 Subject: [PATCH 3/3] Sort particle types and masses --- src/programs/Simulation/genEtaRegge/genEtaRegge.cc | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc index 36d84fb8e..394d6a0b4 100644 --- a/src/programs/Simulation/genEtaRegge/genEtaRegge.cc +++ b/src/programs/Simulation/genEtaRegge/genEtaRegge.cc @@ -827,6 +827,14 @@ int main(int narg, char *argv[]) } cout << endl; } + + // Set specific particle order for eta->e+e-pi+pi- + if (num_decay_particles==4){ + // Order according to particle type and mass + sort(particle_types.begin(),particle_types.end()); + sort(decay_masses.begin(),decay_masses.end(),[&](double a,double b){return a