programs/Simulation/genScalarRegge/genScalarRegge.cc

Bug Summary

File:	programs/Simulation/genScalarRegge/genScalarRegge.cc
Location:	line 2853, column 7
Description:	Value stored to 't' is never read

Annotated Source Code

1	// Main program for generating scalar 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 <iostream>
13	#include <fstream>
14	#include <iomanip>
15	#include <vector>
16	#include <string>
17	using namespace std;
18
19	#include "UTILITIES/BeamProperties.h"
20
21	// Masses
22	const double m_p=0.93827; // GeV
23	const double m_p_sq=m_p*m_p;
24	// Width
25	double width=0.;
26	// Coupling constant
27	double g0_sq=110.5; // GeV^-2
28	// Regge cut parameters
29	double regge_cuts[5]; // dc c_P_omega c_f2_omega c_P_rho c_f2_rho
30
31	double zmin=50.0,zmax=80.0; // cm, target extent
32	int Nevents=10000;
33	int runNo=30300;
34	bool debug=false;
35
36	// Diagnostic histograms
37	TH1D *thrown_t;
38	TH1D *thrown_mass;
39	TH1D *thrown_dalitzZ;
40	TH1D *thrown_Egamma;
41	TH2D *thrown_dalitzXY;
42	TH2D *thrown_theta_vs_p;
43	TH2D *thrown_mass_vs_E;
44	TH1D *cobrems_vs_E;
45
46	char input_file_name[50]="scalar.in";
47	char output_file_name[50]="scalar_gen.hddm";
48
49	void Usage(void){
50	printf("genScalarRegge: generator for eta production based on Regge trajectory formalism.\n");
51	printf(" Usage: genScalarRegge <options>\n");
52	printf(" Options: -N<number of events> (number of events to generate)\n");
53	printf(" -O<output.hddm> (default: scalar_gen.hddm)\n");
54	printf(" -I<input.in> (default: scalar.in)\n");
55	printf(" -R<run number> (default: 30300)\n");
56	printf(" -h (Print this message and exit.)\n");
57	printf("Photon beam energy range, Regge cut parameters, and decay products are\n");
58	printf("specified in the <input.in> file.\n");
59
60	exit(0);
61	}
62
63	//-----------
64	// ParseCommandLineArguments
65	//-----------
66	void ParseCommandLineArguments(int narg, char* argv[])
67	{
68	int seed=0;
69	if (narg==1){
70	Usage();
71	}
72	for(int i=1; i<narg; i++){
73	char *ptr = argv[i];
74
75	if(ptr[0] == '-'){
76	switch(ptr[1]){
77	case 'h': Usage(); break;
78	case 'I':
79	sscanf(&ptr[2],"%s",input_file_name);
80	break;
81	case 'O':
82	sscanf(&ptr[2],"%s",output_file_name);
83	break;
84	case 'N':
85	sscanf(&ptr[2],"%d",&Nevents);
86	break;
87	case 'R':
88	sscanf(&ptr[2],"%d",&runNo);
89	break;
90	case 'S':
91	sscanf(&ptr[2],"%d",&seed);
92	break;
93	case 'd':
94	debug=true;
95	break;
96	default:
97	break;
98	}
99	}
100	}
101	}
102
103	// Non-resonant pi-pi or pi-eta background following Donnachie and Kalashnikova,
104	// arXiv:0806.3698v1
105	double BackgroundCrossSection(TLorentzVector &q /* beam */,
106	vector<Particle_t>&particle_types,
107	vector<TLorentzVector>&particles){
108
109	int two_particles=particle_types[0]+particle_types[1];
110	TLorentzVector p1(0,0,0.,ParticleMass(Proton));
111	TLorentzVector p2=particles[2];
112	TLorentzVector p=p1-p2;
113	double t=p.M2();
114	double s=(q+p1).M2();
115	TLorentzVector v1=particles[0]-q;
116	TLorentzVector v2=particles[1]-q;
117
118	double p1_dot_p2=p1.Dot(p2);
119	double q_dot_p=q.Dot(p);
120	double q_dot_v1=q.Dot(v1);
121	double p_dot_v1=p.Dot(v1);
122	double q_dot_v2=q.Dot(v2);
123	double p_dot_v2=p.Dot(v2);
124	double v1sq=v1.M2();
125	double v2sq=v2.M2();
126	double v1sq_minus_v2sq=v1sq-v2sq;
127	double v1_dot_v2=v1.Dot(v2);
128	double psq=p.M2();
129	double b1=q_dot_pv1sq-q_dot_v1p_dot_v1;
130	double b2=q_dot_pv2sq-q_dot_v2p_dot_v2;
131	TLorentzVector c1=p_dot_v1q-q_dot_pv1;
132	TLorentzVector c2=p_dot_v2q-q_dot_pv2;
133	TLorentzVector d1=q_dot_v1v1-v1sqq;
134	TLorentzVector d2=q_dot_v2v2-v2sqq;
135	TLorentzVector N1=b1p1+p1.Dot(c1)v1+p1.Dot(d1)*p;
136	TLorentzVector N2=b2p1+p1.Dot(c2)v1+p1.Dot(d2)*p;
137	double N1_N1=N1.Dot(N1);
138	double N2_N2=N2.Dot(N2);
139	double N1_N2=N1.Dot(N2);
140	double M1_M2=p_dot_v1p_dot_v1q_dot_v2*q_dot_v2
141	- 2.v1_dot_v2q_dot_p(p_dot_v1q_dot_v2+p_dot_v2*q_dot_v1)
142	+ v1_dot_v2v1_dot_v2q_dot_pq_dot_p + p_dot_v2p_dot_v2q_dot_v1q_dot_v1
143	+ v2sqp_dot_v1q_dot_v1q_dot_p + v1sqp_dot_v2q_dot_v2q_dot_p
144	+ psqv1_dot_v2q_dot_v1q_dot_v2 - psqv2sqq_dot_v1q_dot_v1
145	- psqv1sqq_dot_v2*q_dot_v2;
146	double M1_M1=2.p_dot_v1p_dot_v1q_dot_v1q_dot_v1
147	- 2.v1sqp_dot_v1q_dot_v1q_dot_p + v1sqv1sqq_dot_p*q_dot_p
148	- psqv1sqq_dot_v1*q_dot_v1;
149	double M2_M2=2.p_dot_v2p_dot_v2q_dot_v2q_dot_v2
150	- 2.v2sqp_dot_v2q_dot_v2q_dot_p + v2sqv2sqq_dot_p*q_dot_p
151	- psqv2sqq_dot_v2*q_dot_v2;
152
153	// Rho propagators for top exchange
154	double m_rho=0.7685;
155	double Gamma_rho=0.1462;
156	double m_rhosq_minus_v1sq=m_rho*m_rho-v1sq;
157	double Pi_rho_1_sq=1./(m_rhosq_minus_v1sq*m_rhosq_minus_v1sq
158	+m_rhom_rhoGamma_rho*Gamma_rho);
159	double m_rhosq_minus_v2sq=m_rho*m_rho-v2sq;
160	double Pi_rho_2_sq=1./(m_rhosq_minus_v2sq*m_rhosq_minus_v2sq
161	+m_rhom_rhoGamma_rho*Gamma_rho);
162
163	double s0=1.;
164	// Regge trajectory for rho
165	double a_rho=0.55+0.8*t;
166	double a_rho_prime=0.8;
167	double cos_rho=cos(M_PI3.14159265358979323846*a_rho);
168	double sin_rho=sin(M_PI3.14159265358979323846*a_rho);
169	double regge_rho=pow(s/s0,a_rho-1.)M_PI3.14159265358979323846a_rho_prime/(sin_rho*TMath::Gamma(a_rho)); // excluding phase factor
170	double regge_rho_sq=regge_rho*regge_rho;
171
172	// omega propagators for top exchange
173	double m_omega=0.78265;
174	double Gamma_omega=0.00849;
175	double m_omegasq_minus_v1sq=m_omega*m_omega-v1sq;
176	double Pi_omega_1_sq=1./(m_omegasq_minus_v1sq*m_omegasq_minus_v1sq
177	+m_omegam_omegaGamma_omega*Gamma_omega);
178	double m_omegasq_minus_v2sq=m_omega*m_omega-v2sq;
179	double Pi_omega_2_sq=1./(m_omegasq_minus_v2sq*m_omegasq_minus_v2sq
180	+m_omegam_omegaGamma_omega*Gamma_omega);
181
182	// Regge trajectory for omega
183	double a_omega=0.44+0.9*t;
184	double a_omega_prime=0.9;
185	double cos_omega=cos(M_PI3.14159265358979323846*a_omega);
186	double sin_omega=sin(M_PI3.14159265358979323846*a_omega);
187	double regge_omega=pow(s/s0,a_omega-1.)M_PI3.14159265358979323846a_omega_prime/(sin_omega*TMath::Gamma(a_omega)); // excluding phase factor
188	double regge_omega_sq=regge_omega*regge_omega;
189
190	// rho-omega interference
191	double regge_rho_omega=regge_rho*regge_omega;
192	double cos_rho_omega=cos(M_PI3.14159265358979323846*(a_rho-a_omega));
193	double sin_rho_omega=sin(M_PI3.14159265358979323846*(a_rho-a_omega));
194
195	// Coupling constants
196	double g_omega_V=15.;
197	double gsq_omega_V=g_omega_V*g_omega_V;
198	double g_rho_V=3.4;
199	double g_rho_T=11.0; // GeV^-1
200	double gsq_rho_T=g_rho_T*g_rho_T;
201	double g_rho_V_and_T=g_rho_V+2.m_pg_rho_T;
202	double gsq_rho_V_and_T=g_rho_V_and_T*g_rho_V_and_T;
203
204	// terms involving complex conjugates of Regge propagators and rho/omega propagtors
205	double a1_a1=0.,b1_b1=0.;
206	double a2_a2=0.,b2_b2=0.;
207	double a1_a2=0.,b1_b2=0.;
208	double b1_a1=0,b2_a2=0.,b1_a2_b2_a1=0.;
209	double Csq=1.,zetasq=1./2.;
210	// Compute square of amplitude
211	if (two_particles==(7+7)){ // pi0 pi0
212	double a1_a1_rho_V_and_T_sq=gsq_rho_V_and_Tregge_rho_sqPi_omega_1_sq;
213	double a1_a1_omega_V_sq=gsq_omega_Vregge_omega_sqPi_rho_1_sq;
214	double a1_a1_omega_V_rho_V_and_T
215	=g_rho_V_and_Tg_omega_Vregge_rho_omega
216	Pi_rho_1_sqPi_omega_1_sq*( cos_rho_omega
217	(m_rhosq_minus_v1sqm_omegasq_minus_v1sq
218	+m_rhom_omegaGamma_rho*Gamma_omega)
219	+ sin_rho_omega // check sign
220	(m_rhoGamma_rho*m_omegasq_minus_v1sq
221	-m_omegaGamma_omegam_rhosq_minus_v1sq));
222	double a2_a2_rho_V_and_T_sq=gsq_rho_V_and_Tregge_rho_sqPi_omega_2_sq;
223	double a2_a2_omega_V_sq=gsq_omega_Vregge_omega_sqPi_rho_2_sq;
224	double a2_a2_omega_V_rho_V_and_T=g_rho_V_and_Tg_omega_Vregge_rho_omega
225	Pi_rho_2_sqPi_omega_2_sq*( cos_rho_omega
226	(m_rhosq_minus_v2sqm_omegasq_minus_v2sq
227	+m_rhom_omegaGamma_rho*Gamma_omega)
228	+ sin_rho_omega // check sign
229	(m_rhoGamma_rho*m_omegasq_minus_v2sq
230	-m_omegaGamma_omegam_rhosq_minus_v2sq));
231	double a1_a2_rho_V_and_T_sq
232	=gsq_rho_V_and_Tregge_rhoregge_rhoPi_omega_1_sqPi_omega_2_sq
233	(1.-cos_rho)(m_omegasq_minus_v1sq*m_omegasq_minus_v2sq
234	+ m_omegaGamma_omegam_omega*Gamma_omega);
235	double a1_a2_omega_V_sq
236	=gsq_omega_Vregge_omegaregge_omegaPi_rho_1_sqPi_rho_2_sq
237	(1.-cos_omega)(m_rhosq_minus_v1sq*m_rhosq_minus_v2sq
238	+m_rhom_rhoGamma_rho*Gamma_rho);
239	double a1_a2_omega_V_rho_V_and_T_term1=g_rho_V_and_T*g_omega_V
240	Pi_omega_1_sqPi_rho_2_sqregge_rhoregge_omega
241	(cos_rho_omega(m_omegasq_minus_v1sq*m_rhosq_minus_v2sq
242	+m_rhoGamma_rhom_omega*Gamma_omega)
243	+sin_rho_omega(m_omegaGamma_omega*m_rhosq_minus_v2sq
244	-m_rhoGamma_rhom_omegasq_minus_v1sq));
245	double a1_a2_omega_V_rho_V_and_T_term2=g_rho_V_and_T*g_omega_V
246	Pi_omega_1_sqPi_omega_2_sqregge_rhoregge_omega
247	(cos_rho_omega(m_omegasq_minus_v1sq*m_omegasq_minus_v2sq
248	+m_omegaGamma_omegam_omega*Gamma_omega)
249	-sin_rho_omegam_omegaGamma_omega*v1sq_minus_v2sq);
250	double b1_a1_rho_T_rho_V_and_T
251	=-4.g_rho_Tg_rho_V_and_Tregge_rho_sqPi_omega_1_sq;
252	double b1_a1_omega_V_rho_V_and_T
253	=-g_rho_Tg_omega_VPi_rho_1_sqPi_omega_1_sqregge_omega*regge_rho
254	(cos_rho_omega(m_rhosq_minus_v1sq*m_omegasq_minus_v1sq
255	+m_rhoGamma_rhom_omega*Gamma_omega)
256	+sin_rho_omega(m_omegaGamma_omega*m_rhosq_minus_v1sq
257	-m_rhoGamma_rhom_omegasq_minus_v1sq));
258	double b2_a2_rho_T_rho_V_and_T
259	=-4.g_rho_Tg_rho_V_and_Tregge_rho_sqPi_omega_2_sq;
260	double b2_a2_omega_V_rho_V_and_T
261	=-g_rho_Tg_omega_VPi_rho_2_sqPi_omega_2_sqregge_rho*regge_omega
262	(cos_rho_omega(m_rhosq_minus_v2sq*m_omegasq_minus_v2sq
263	+m_rhoGamma_rhom_omega*Gamma_omega)
264	+sin_rho_omega(m_omegaGamma_omega*m_rhosq_minus_v2sq
265	-m_rhoGamma_rhom_omegasq_minus_v2sq));
266	double b1_a2_rho_T_rho_V_and_T
267	=-2.regge_rhoregge_rhog_rho_Tg_rho_V_and_T
268	Pi_omega_1_sqPi_omega_2_sq
269	(1.-cos_rho)(m_omegasq_minus_v1sq*m_omegasq_minus_v2sq
270	+m_omegaGamma_omegam_omega*Gamma_omega);
271	double b1_a2_omega_V_rho_T=-g_rho_Tg_omega_VPi_omega_1_sq*Pi_rho_2_sq
272	regge_rhoregge_omega
273	(cos_rho_omega(m_omegasq_minus_v1sq*m_rhosq_minus_v2sq
274	+m_rhom_omegaGamma_omega*Gamma_rho)
275	+sin_rho_omega(m_omegaGamma_omega*m_rhosq_minus_v2sq
276	-m_rhoGamma_rhom_omegasq_minus_v1sq));
277	double b2_a1_rho_T_rho_V_and_T=b1_a2_rho_T_rho_V_and_T;
278	double b2_a1_omega_V_rho_T=-g_rho_Tg_omega_VPi_rho_1_sq*Pi_omega_2_sq
279	regge_rhoregge_omega
280	(cos_rho_omega(m_rhosq_minus_v1sq*m_omegasq_minus_v2sq
281	+m_rhom_omegaGamma_rho*Gamma_omega)
282	+sin_rho_omega(m_omegaGamma_omega*m_rhosq_minus_v1sq
283	-m_rhoGamma_rhom_omegasq_minus_v2sq));
284
285	b1_b1=4.gsq_rho_Tregge_rho_sq*Pi_omega_1_sq;
286	b2_b2=4.gsq_rho_Tregge_rho_sq*Pi_omega_2_sq;
287	b1_b2=4.gsq_rho_Tregge_rho*regge_rho
288	Pi_omega_1_sqPi_omega_2_sq
289	(1.-cos_rho)(m_omegasq_minus_v1sq*m_omegasq_minus_v2sq
290	+m_omegam_rhoGamma_rho*Gamma_rho);
291
292	a1_a1=a1_a1_rho_V_and_T_sq + (1./9.)*a1_a1_omega_V_sq
293	+ (1./6.)*a1_a1_omega_V_rho_V_and_T;
294	a2_a2=a2_a2_rho_V_and_T_sq + (1./9.)*a2_a2_omega_V_sq
295	+ (1./6.)*a2_a2_omega_V_rho_V_and_T;
296	a1_a2=a1_a2_rho_V_and_T_sq + (1./18.)*a1_a2_omega_V_sq
297	+(1./6.)*(a1_a2_omega_V_rho_V_and_T_term1
298	+a1_a2_omega_V_rho_V_and_T_term2);
299
300	b1_a1=b1_a1_rho_T_rho_V_and_T+(1./3.)*b1_a1_omega_V_rho_V_and_T;
301	b2_a2=b2_a2_rho_T_rho_V_and_T+(1./3.)*b2_a2_omega_V_rho_V_and_T;
302	b1_a2_b2_a1=b1_a2_rho_T_rho_V_and_T+b2_a1_rho_T_rho_V_and_T
303	+(1./3.)*(b1_a2_omega_V_rho_T+b2_a1_omega_V_rho_T);
304
305	}
306	else if (two_particles==(7+17)){ // pi0 eta
307	Csq=2./3.;
308	zetasq=1.;
309
310	double a1_a1_rho_V_and_T_sq=gsq_rho_V_and_Tregge_rho_sqPi_rho_1_sq;
311	double a1_a1_omega_V_sq=gsq_omega_Vregge_omega_sqPi_omega_1_sq;
312	double a1_a1_omega_V_rho_V_and_T=g_rho_V_and_Tg_omega_Vregge_rho_omega
313	Pi_rho_1_sqPi_omega_1_sq*( cos_rho_omega
314	(m_rhosq_minus_v1sqm_omegasq_minus_v1sq
315	+m_rhom_omegaGamma_rho*Gamma_omega)
316	+sin_rho_omega
317	(m_rhoGamma_rho*m_omegasq_minus_v1sq
318	-m_omegaGamma_omegam_rhosq_minus_v1sq));
319	double a2_a2_rho_V_and_T_sq=gsq_rho_V_and_Tregge_rho_sqPi_omega_2_sq;
320	double a2_a2_omega_V_sq=gsq_omega_Vregge_omega_sqPi_rho_2_sq;
321	double a2_a2_omega_V_rho_V_and_T=g_rho_V_and_Tg_omega_Vregge_rho_omega
322	Pi_rho_2_sqPi_omega_2_sq*( cos_rho_omega
323	(m_rhosq_minus_v2sqm_omegasq_minus_v2sq
324	+m_rhom_omegaGamma_rho*Gamma_omega)
325	+ sin_rho_omega
326	(m_rhoGamma_rho*m_omegasq_minus_v2sq
327	-m_omegaGamma_omegam_rhosq_minus_v2sq));
328	double a1_a2_rho_V_and_T_sq
329	=2.gsq_rho_V_and_Tregge_rhoregge_rhoPi_rho_1_sq*Pi_omega_2_sq
330	(1.-cos_rho)(m_rhosq_minus_v1sq*m_omegasq_minus_v2sq
331	+ m_rhoGamma_rhom_omega*Gamma_omega);
332	double a1_a2_omega_V_rho_V_and_T_term1=g_rho_V_and_T*g_omega_V
333	Pi_rho_1_sqPi_rho_2_sqregge_rhoregge_omega
334	*(cos_rho_omega
335	(m_rhosq_minus_v1sqm_rhosq_minus_v2sq+m_rhoGamma_rhom_rho*Gamma_rho)
336	+sin_rho_omegam_rhoGamma_rho*v1sq_minus_v2sq); // check sign!
337	double a1_a2_omega_V_rho_V_and_T_term2=g_rho_V_and_T*g_omega_V
338	Pi_omega_1_sqPi_omega_2_sqregge_rhoregge_omega
339	(cos_rho_omega(m_omegasq_minus_v1sq*m_omegasq_minus_v2sq
340	+m_omegaGamma_omegam_omega*Gamma_omega)
341	+sin_rho_omegam_omegaGamma_omega*v1sq_minus_v2sq);
342	double a1_a2_omega_V_sq
343	=gsq_omega_Vregge_omegaregge_omegaPi_omega_1_sqPi_rho_2_sq
344	2.(1.-cos_omega)(m_omegasq_minus_v1sqm_rhosq_minus_v2sq
345	+m_rhom_omegaGamma_rho*Gamma_omega);
346	double b1_a1_rho_T_rho_V_and_T
347	=-4.g_rho_Tg_rho_V_and_Tregge_rho_sqPi_rho_1_sq;
348	double b1_a1_omega_V_rho_V_and_T
349	=-g_rho_Tg_omega_VPi_rho_1_sqPi_omega_1_sqregge_omega*regge_rho
350	(cos_rho_omega(m_rhosq_minus_v1sq*m_omegasq_minus_v1sq
351	+m_rhoGamma_rhom_omega*Gamma_omega)
352	+sin_rho_omega(m_rhoGamma_rho*m_omegasq_minus_v1sq
353	-m_omegaGamma_omegam_rhosq_minus_v1sq));
354	double b2_a2_rho_T_rho_V_and_T
355	=-4.g_rho_Tg_rho_V_and_Tregge_rho_sqPi_omega_2_sq;
356	double b2_a2_omega_V_rho_V_and_T
357	=-g_rho_Tg_omega_VPi_rho_2_sqPi_omega_2_sqregge_rho*regge_omega
358	(cos_rho_omega(m_rhosq_minus_v2sq*m_omegasq_minus_v2sq
359	+m_rhoGamma_rhom_omega*Gamma_omega)
360	+sin_rho_omega(m_omegaGamma_omega*m_rhosq_minus_v2sq
361	-m_rhoGamma_rhom_omegasq_minus_v2sq));
362	double b1_a2_rho_T_rho_V_and_T
363	=-regge_rhoregge_rhog_rho_T*g_rho_V_and_T
364	Pi_rho_1_sqPi_omega_2_sq2.(1.-cos_rho)
365	(m_rhosq_minus_v1sqm_omegasq_minus_v2sq
366	+m_omegaGamma_omegam_rho*Gamma_rho);
367	double b1_a2_omega_V_rho_T=-g_rho_Tg_omega_VPi_rho_1_sq*Pi_rho_2_sq
368	regge_rhoregge_omega*(cos_rho_omega
369	(m_rhosq_minus_v1sqm_rhosq_minus_v2sq
370	+m_rhom_rhoGamma_rho*Gamma_rho)
371	-sin_rho_omega // check sign!
372	m_rhoGamma_rho*(m_rhosq_minus_v2sq
373	-m_rhosq_minus_v1sq));
374	double b2_a1_rho_T_rho_V_and_T=b1_a2_rho_T_rho_V_and_T;
375	double b2_a1_omega_V_rho_T=-g_rho_Tg_omega_VPi_omega_1_sq*Pi_omega_2_sq
376	regge_rhoregge_omega*(cos_rho_omega
377	(m_omegasq_minus_v1sqm_omegasq_minus_v2sq
378	+m_omegam_omegaGamma_omega*Gamma_omega)
379	-sin_rho_omega
380	m_omegaGamma_omega*(m_omegasq_minus_v1sq
381	-m_omegasq_minus_v2sq));
382
383	b1_b1=(4./9.)gsq_rho_Tregge_rho_sq*Pi_rho_1_sq;
384	b2_b2=(4./9.)gsq_rho_Tregge_rho_sq*Pi_omega_2_sq;
385	b1_b2=(2./9.)gsq_rho_Tregge_rho*regge_rho
386	Pi_rho_1_sqPi_omega_2_sq2.(1.-cos_rho)
387	(m_rhosq_minus_v1sqm_rhosq_minus_v2sq
388	+m_rhom_rhoGamma_rho*Gamma_rho);
389
390	a1_a1=(1./9.)*a1_a1_rho_V_and_T_sq + a1_a1_omega_V_sq
391	+ (1./6.)*a1_a1_omega_V_rho_V_and_T;
392	a2_a2=(1./9.)*a2_a2_rho_V_and_T_sq + a2_a2_omega_V_sq
393	+ (1./6.)*a2_a2_omega_V_rho_V_and_T;
394	a1_a2=(1./18.)a1_a2_rho_V_and_T_sq + (1/2.)a1_a2_omega_V_sq
395	+(1./6.)*(a1_a2_omega_V_rho_V_and_T_term1
396	+a1_a2_omega_V_rho_V_and_T_term2);
397
398	b1_a1=(1./9.)b1_a1_rho_T_rho_V_and_T+(1./3.)b1_a1_omega_V_rho_V_and_T;
399	b2_a2=(1./9.)b2_a2_rho_T_rho_V_and_T+(1./3.)b2_a2_omega_V_rho_V_and_T;
400	b1_a2_b2_a1=(1./9.)*(b1_a2_rho_T_rho_V_and_T+b2_a1_rho_T_rho_V_and_T)
401	+(1./3.)*(b1_a2_omega_V_rho_T+b2_a1_omega_V_rho_T);
402	}
403
404	double T=Csq((m_p_sq-p1_dot_p2)(a1_a1M1_M1 + a1_a2M1_M2 + a2_a2*M2_M2)
405	+ 2.(a1_a1N1_N1 + a1_a2N1_N2 + a2_a2N2_N2)
406	+ 2.m_p(b1_a1N1_N1 + b1_a2_b2_a1N1_N2 + b2_a2*N2_N2)
407	+ (m_p_sq+p1_dot_p2)(b1_b1N1_N1 + b1_b2N1_N2 + b2_b2N2_N2));
408
409	double Msq=(particles[0]+particles[1]).M2();
410	double s_minus_mp_sq=s-m_p_sq;
411	double m1=ParticleMass(particle_types[0]);
412	double m2=ParticleMass(particle_types[1]);
413	double m1_plus_m2=m1+m2;
414	double m1_minus_m2=m1-m2;
415	double k=sqrt((Msq-m1_plus_m2m1_plus_m2)(Msq-m1_minus_m2*m1_minus_m2))
416	/(2.*sqrt(Msq));
417	double hbarc_sq=389.; // Convert to micro-barns
418	double xsec=-zetasqg0_sqkThbarc_sq
419	/(256M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846s_minus_mp_sqs_minus_mp_sq);
420
421	return xsec;
422	}
423
424	double InterferenceCrossSection(TLorentzVector &q /* beam */,
425	vector<Particle_t>&particle_types,
426	vector<TLorentzVector>&particles,
427	double gR,double ReB, double ImB,
428	double gsq_rho_S,double gsq_omega_S,
429	double phase){
430	int two_particles=particle_types[0]+particle_types[1];
431
432	// Four vectors
433	TLorentzVector p1(0,0,0.,ParticleMass(Proton));
434	TLorentzVector p2=particles[2];
435	TLorentzVector p=p1-p2;
436	TLorentzVector v1=particles[0]-q;
437	TLorentzVector v2=particles[1]-q;
438
439	// Mandelstam variables
440	double s=(q+p1).M2();
441	double t=(p1-p2).M2();
442
443	// dot products
444	double p1_dot_p2=p1.Dot(p2);
445	double q_dot_p=q.Dot(p);
446	double q_dot_v1=q.Dot(v1);
447	double p_dot_v1=p.Dot(v1);
448	double q_dot_v2=q.Dot(v2);
449	double p_dot_v2=p.Dot(v2);
450	double v1sq=v1.M2();
451	double v2sq=v2.M2();
452	double psq=p.M2();
453	double b1=q_dot_pv1sq-q_dot_v1p_dot_v1;
454	double b2=q_dot_pv2sq-q_dot_v2p_dot_v2;
455	TLorentzVector c1=p_dot_v1q-q_dot_pv1;
456	TLorentzVector c2=p_dot_v2q-q_dot_pv2;
457	TLorentzVector d1=q_dot_v1v1-v1sqq;
458	TLorentzVector d2=q_dot_v2v2-v2sqq;
459	TLorentzVector N1=b1p1+p1.Dot(c1)v1+p1.Dot(d1)*p;
460	TLorentzVector N2=b2p1+p1.Dot(c2)v1+p1.Dot(d2)*p;
461	TLorentzVector N=(q_dot_p)p1-q.Dot(p1)p;
462	double N_N1=N.Dot(N1);
463	double N_N2=N.Dot(N2);
464	double M_M1=3b1q_dot_p+v1.Dot(c1)q_dot_p+p_dot_v1q.Dot(c1)
465	+p.Dot(d1)q_dot_p-psqq.Dot(d1);
466	double M_M2=3b2q_dot_p+v2.Dot(c2)q_dot_p+p_dot_v2q.Dot(c2)
467	+p.Dot(d2)q_dot_p-psqq.Dot(d2);
468
469	// Coupling constants
470	double g_omega_V=15.;
471	double gsq_omega_V=g_omega_V*g_omega_V;
472	double g_rho_V=3.4;
473	double g_rho_T=11.0; // GeV^-1
474	double gsq_rho_T=g_rho_T*g_rho_T;
475	double g_rho_V_and_T=g_rho_V+2.m_pg_rho_T;
476	double gsq_rho_V_and_T=g_rho_V_and_T*g_rho_V_and_T;
477	double g_rho_S=sqrt(gsq_rho_S);
478	double g_omega_S=sqrt(gsq_omega_S);
479
480	// Rho propagator for top exchange for double-exchange diagrams
481	double m_rho=0.7685;
482	double Gamma_rho=0.1462;
483	double m_rhosq_minus_v1sq=m_rho*m_rho-v1sq;
484	double Pi_rho_1_sq=1./(m_rhosq_minus_v1sq*m_rhosq_minus_v1sq
485	+m_rhom_rhoGamma_rho*Gamma_rho);
486	double m_rhosq_minus_v2sq=m_rho*m_rho-v2sq;
487	double Pi_rho_2_sq=1./(m_rhosq_minus_v2sq*m_rhosq_minus_v2sq
488	+m_rhom_rhoGamma_rho*Gamma_rho);
489
490	double s0=1.;
491	// Regge trajectory for rho
492	double a_rho=0.55+0.8*t;
493	double a_rho_prime=0.8;
494	double cos_rho=cos(M_PI3.14159265358979323846*a_rho);
495	double one_minus_cos_rho=1.-cos_rho;
496	double sin_rho=sin(M_PI3.14159265358979323846*a_rho);
497	double regge_rho=pow(s/s0,a_rho-1.)M_PI3.14159265358979323846a_rho_prime/(sin_rho*TMath::Gamma(a_rho)); // excluding phase factor
498	double regge_rho_sq=regge_rho*regge_rho;
499
500	// Regge cuts for rho
501	double dc=regge_cuts[0];
502	double a_rho_P=0.64+0.16*t; // Pomeron
503	double a_rho_f2=0.222+0.404*t;
504	double regge_rho_P_cut=exp(dct)pow(s/s0,a_rho_P-1.);
505	double regge_rho_f2_cut=exp(dct)pow(s/s0,a_rho_f2-1.);
506	double C_rho_P_cut=regge_cuts[3];
507	double C_rho_f2_cut=regge_cuts[4];
508
509	// omega propagator for top exchange
510	double m_omega=0.78265;
511	double Gamma_omega=0.00849;
512	double m_omegasq_minus_v1sq=m_omega*m_omega-v1sq;
513	double Pi_omega_1_sq=1./(m_omegasq_minus_v1sq*m_omegasq_minus_v1sq
514	+m_omegam_omegaGamma_omega*Gamma_omega);
515	double m_omegasq_minus_v2sq=m_omega*m_omega-v2sq;
516	double Pi_omega_2_sq=1./(m_omegasq_minus_v2sq*m_omegasq_minus_v2sq
517	+m_omegam_omegaGamma_omega*Gamma_omega);
518
519	// Regge trajectory for omega
520	double a_omega=0.44+0.9*t;
521	double a_omega_prime=0.9;
522	double cos_omega=cos(M_PI3.14159265358979323846*a_omega);
523	double one_minus_cos_omega=1.-cos_omega;
524	double sin_omega=sin(M_PI3.14159265358979323846*a_omega);
525	double regge_omega=pow(s/s0,a_omega-1.)M_PI3.14159265358979323846a_omega_prime/(sin_omega*TMath::Gamma(a_omega)); // excluding phase factor
526	double regge_omega_sq=regge_omega*regge_omega;
527
528	// Regge cuts for omega
529	double a_omega_P=0.52+0.196*t; // Pomeron
530	double a_omega_f2=0.112+0.428*t;
531	double regge_omega_P_cut=exp(dct)pow(s/s0,a_omega_P-1.);
532	double regge_omega_f2_cut=exp(dct)pow(s/s0,a_omega_f2-1.);
533	double C_omega_P_cut=regge_cuts[1];
534	double C_omega_f2_cut=regge_cuts[2];
535
536	// cut interference factors
537	double cos_rho_rho_P=cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_P));
538	double cos_rho_rho_f2=cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2));
539	double sin_rho_rho_P=sin(M_PI3.14159265358979323846(a_rho-0.5a_rho_P));
540	double sin_rho_rho_f2=sin(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2));
541	double cos_omega_rho_P=cos(M_PI3.14159265358979323846(a_omega-0.5a_rho_P));
542	double cos_omega_rho_f2=cos(M_PI3.14159265358979323846(a_omega-0.5a_rho_f2));
543	double sin_omega_rho_P=sin(M_PI3.14159265358979323846(a_omega-0.5a_rho_P));
544	double sin_omega_rho_f2=sin(M_PI3.14159265358979323846(a_omega-0.5a_rho_f2));
545	double cos_rho_omega_P=cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_P));
546	double cos_rho_omega_f2=cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2));
547	double sin_rho_omega_P=sin(M_PI3.14159265358979323846(a_rho-0.5a_omega_P));
548	double sin_rho_omega_f2=sin(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2));
549	double cos_omega_omega_P=cos(M_PI3.14159265358979323846(a_omega-0.5a_omega_P));
550	double cos_omega_omega_f2=cos(M_PI3.14159265358979323846(a_omega-0.5a_omega_f2));
551	double sin_omega_omega_P=sin(M_PI3.14159265358979323846(a_omega-0.5a_omega_P));
552	double sin_omega_omega_f2=sin(M_PI3.14159265358979323846(a_omega-0.5a_omega_f2));
553
554	// rho-omega interference
555	double cos_rho_omega_sum1=cos(M_PI3.14159265358979323846(a_rho-a_omega))-cos(M_PI3.14159265358979323846a_omega);
556	double cos_rho_omega_sum2=cos(M_PI3.14159265358979323846(a_omega-a_rho))-cos(M_PI3.14159265358979323846a_rho);
557	double sin_rho_omega_sum1=sin(M_PI3.14159265358979323846(a_rho-a_omega))+sin(M_PI3.14159265358979323846a_omega);
558	double sin_rho_omega_sum2=sin(M_PI3.14159265358979323846(a_omega-a_rho))+sin(M_PI3.14159265358979323846a_rho);
559
560	// terms for interference between resonance shape and propagator in top
561	// part of background diagrams
562	double Re_B_and_omega_1=m_omegasq_minus_v1sqReB+m_omegaGamma_omega*ImB;
563	double Im_B_and_omega_1=m_omegaGamma_omegaReB-m_omegasq_minus_v1sq*ImB;
564	double Re_B_and_omega_2=m_omegasq_minus_v2sqReB+m_omegaGamma_omega*ImB;
565	double Im_B_and_omega_2=m_omegaGamma_omegaReB-m_omegasq_minus_v2sq*ImB;
566	double Re_B_and_rho_1=m_rhosq_minus_v1sqReB+m_rhoGamma_rho*ImB;
567	double Im_B_and_rho_1=m_rhoGamma_rhoReB-m_rhosq_minus_v1sq*ImB;
568	double Re_B_and_rho_2=m_rhosq_minus_v2sqReB+m_rhoGamma_rho*ImB;
569	double Im_B_and_rho_2=m_rhoGamma_rhoReB-m_rhosq_minus_v2sq*ImB;
570
571	// terms involving complex conjugates of Regge propagators and rho/omega propagtors
572	double Re_a1_aS=0, Im_a1_aS=0, Re_a2_aS=0., Im_a2_aS=0.;
573	double Re_b1_aS=0, Im_b1_aS=0, Re_b2_aS=0., Im_b2_aS=0.;
574	double Re_a1_bS=0, Im_a1_bS=0., Re_b2_bS=0., Im_b2_bS=0.;
575	double Re_a2_bS=0, Im_a2_bS=0., Re_b1_bS=0., Im_b1_bS=0.;
576
577	// Compute imaginary and real contributions
578	double zeta=1., C=1;
579	if (two_particles==(7+7)){ // pi0 pi0
580	zeta=1./sqrt(2.);
581
582	double Re_a1_aS_rho_V_and_T_sq
583	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_1_sq
584	(Re_B_and_omega_1one_minus_cos_rho-Im_B_and_omega_1*sin_rho);
585	double Im_a1_aS_rho_V_and_T_sq
586	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_1_sq
587	(Re_B_and_omega_1sin_rho+Im_B_and_omega_1*one_minus_cos_rho);
588	double Re_a1_aS_rho_V_and_T_rho_cuts
589	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_omega_1_sq
590	(C_rho_P_cutregge_rho_P_cut
591	(cos_rho_rho_PRe_B_and_omega_1
592	+sin_rho_rho_P*Im_B_and_omega_1)
593	+ C_rho_f2_cut*regge_rho_f2_cut
594	(cos_rho_rho_f2Re_B_and_omega_1
595	+sin_rho_rho_f2*Im_B_and_omega_1));
596	double Im_a1_aS_rho_V_and_T_rho_cuts
597	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_omega_1_sq
598	(C_rho_P_cutregge_rho_P_cut
599	(cos_rho_rho_PIm_B_and_omega_1
600	-sin_rho_rho_P*Re_B_and_omega_1)
601	+ C_rho_f2_cut*regge_rho_f2_cut
602	(cos_rho_rho_f2Im_B_and_omega_1
603	-sin_rho_rho_f2*Re_B_and_omega_1));
604	double Re_a1_aS_omega_V_rho_V_and_T
605	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_omega_1_sq
606	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_omega_1
607	-sin_rho_omega_sum2*Im_B_and_omega_1);
608	double Im_a1_aS_omega_V_rho_V_and_T
609	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_omega_1_sq
610	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_omega_1
611	+sin_rho_omega_sum2*Re_B_and_omega_1);
612	double Re_a1_aS_rho_V_and_T_omega_cuts
613	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_omega_1_sq
614	(C_omega_P_cutregge_omega_P_cut
615	(cos_rho_omega_PRe_B_and_omega_1
616	+sin_rho_omega_P*Im_B_and_omega_1)
617	+ C_omega_f2_cut*regge_omega_f2_cut
618	(cos_rho_omega_f2Re_B_and_omega_1
619	+sin_rho_omega_f2*Im_B_and_omega_1));
620	double Im_a1_aS_rho_V_and_T_omega_cuts
621	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_omega_1_sq
622	(C_omega_P_cutregge_omega_P_cut
623	(cos_rho_omega_PIm_B_and_omega_1
624	-sin_rho_omega_P*Re_B_and_omega_1)
625	+ C_rho_f2_cut*regge_omega_f2_cut
626	(cos_rho_omega_f2Im_B_and_omega_1
627	-sin_rho_omega_f2*Re_B_and_omega_1));
628	double Re_a1_aS_omega_V_rho_V_and_T_rho_propagator
629	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_rho_1_sq
630	regge_rhoregge_omega(cos_rho_omega_sum1Re_B_and_rho_1
631	-sin_rho_omega_sum1*Im_B_and_rho_1);
632	double Im_a1_aS_omega_V_rho_V_and_T_rho_propagator
633	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_rho_1_sq
634	regge_rhoregge_omega(cos_rho_omega_sum1Im_B_and_rho_1
635	+sin_rho_omega_sum1*Re_B_and_rho_1);
636	double Re_a1_aS_omega_V_rho_V_rho_cuts_rho_propagator
637	=2.g_omega_Vg_rho_Vg_rho_SPi_rho_1_sq*regge_omega
638	(C_rho_P_cutregge_rho_P_cut
639	(cos_omega_rho_PRe_B_and_rho_1
640	+sin_omega_rho_P*Im_B_and_rho_1)
641	+ C_rho_f2_cut*regge_rho_f2_cut
642	(cos_omega_rho_f2Re_B_and_rho_1
643	+sin_omega_rho_f2*Im_B_and_rho_1));
644	double Im_a1_aS_omega_V_rho_V_rho_cuts_rho_propagator
645	=2.g_omega_Vg_rho_Vg_rho_SPi_rho_1_sq*regge_omega
646	(C_rho_P_cutregge_rho_P_cut
647	(cos_omega_rho_PIm_B_and_rho_1
648	-sin_omega_rho_P*Re_B_and_rho_1)
649	+ C_rho_f2_cut*regge_rho_f2_cut
650	(cos_omega_rho_f2Im_B_and_rho_1
651	-sin_omega_rho_f2*Re_B_and_rho_1));
652	double Re_a1_aS_omega_V_sq
653	=g_omega_Vg_omega_Vg_omega_SPi_rho_1_sqregge_omega_sq
654	(Re_B_and_rho_1one_minus_cos_omega-Im_B_and_rho_1*sin_omega);
655	double Im_a1_aS_omega_V_sq
656	=g_omega_Vg_omega_Vg_omega_SPi_rho_1_sqregge_omega_sq
657	(Re_B_and_rho_1sin_omega+Im_B_and_rho_1*one_minus_cos_omega);
658	double Re_a1_aS_omega_V_sq_omega_cuts
659	=2.g_omega_Vg_omega_Vg_omega_SPi_rho_1_sq*regge_omega
660	(C_omega_P_cutregge_omega_P_cut
661	(cos_omega_omega_PRe_B_and_rho_1
662	+sin_omega_omega_P*Im_B_and_rho_1)
663	+ C_omega_f2_cut*regge_omega_f2_cut
664	(cos_omega_omega_f2Re_B_and_rho_1
665	+sin_omega_omega_f2*Im_B_and_rho_1));
666	double Im_a1_aS_omega_V_sq_omega_cuts
667	=2.g_omega_Vg_omega_Vg_omega_SPi_rho_1_sq*regge_omega
668	(C_omega_P_cutregge_omega_P_cut
669	(cos_omega_omega_PIm_B_and_rho_1
670	-sin_omega_omega_P*Re_B_and_rho_1)
671	+ C_omega_f2_cut*regge_omega_f2_cut
672	(cos_omega_omega_f2Im_B_and_rho_1
673	-sin_omega_omega_f2*Re_B_and_rho_1));
674
675	double Re_a2_aS_rho_V_and_T_sq
676	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
677	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
678	double Im_a2_aS_rho_V_and_T_sq
679	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
680	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
681	double Re_a2_aS_rho_V_and_T_rho_cuts
682	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_omega_2_sq
683	(C_rho_P_cutregge_omega_P_cut
684	(cos_rho_rho_PRe_B_and_omega_2
685	+sin_rho_rho_P*Im_B_and_omega_2)
686	+ C_rho_f2_cut*regge_omega_f2_cut
687	(cos_rho_rho_f2Re_B_and_omega_2
688	+sin_rho_rho_f2*Im_B_and_omega_2));
689	double Im_a2_aS_rho_V_and_T_rho_cuts
690	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_omega_2_sq
691	(C_rho_P_cutregge_omega_P_cut
692	(cos_rho_rho_PIm_B_and_omega_2
693	-sin_rho_rho_P*Re_B_and_omega_2)
694	+ C_rho_f2_cut*regge_omega_f2_cut
695	(cos_rho_rho_f2Im_B_and_omega_2
696	-sin_rho_rho_f2*Re_B_and_omega_2));
697	double Re_a2_aS_omega_V_rho_V_and_T
698	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_omega_2_sq
699	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_omega_2
700	-sin_rho_omega_sum2*Im_B_and_omega_2);
701	double Im_a2_aS_omega_V_rho_V_and_T
702	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_omega_2_sq
703	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_omega_2
704	+sin_rho_omega_sum2*Re_B_and_omega_2);
705	double Re_a2_aS_rho_V_and_T_omega_cuts
706	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_omega_2_sq
707	(C_omega_P_cutregge_omega_P_cut
708	(cos_rho_omega_PRe_B_and_omega_2
709	+sin_rho_omega_P*Im_B_and_omega_2)
710	+ C_omega_f2_cut*regge_omega_f2_cut
711	(cos_rho_omega_f2Re_B_and_omega_2
712	+sin_rho_omega_f2*Im_B_and_omega_2));
713	double Im_a2_aS_rho_V_and_T_omega_cuts
714	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_omega_2_sq
715	(C_omega_P_cutregge_omega_P_cut
716	(cos_rho_omega_PIm_B_and_omega_2
717	-sin_rho_omega_P*Re_B_and_omega_2)
718	+ C_rho_f2_cut*regge_omega_f2_cut
719	(cos_rho_omega_f2Im_B_and_omega_2
720	-sin_rho_omega_f2*Re_B_and_omega_2));
721	double Re_a2_aS_omega_V_rho_V_and_T_rho_propagator
722	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_rho_2_sq
723	regge_rhoregge_omega(cos_rho_omega_sum1Re_B_and_rho_2
724	-sin_rho_omega_sum1*Im_B_and_rho_2);
725	double Im_a2_aS_omega_V_rho_V_and_T_rho_propagator
726	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_rho_2_sq
727	regge_rhoregge_omega(cos_rho_omega_sum1Im_B_and_rho_2
728	+sin_rho_omega_sum1*Re_B_and_rho_2);
729	double Re_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
730	=2.g_omega_Vg_rho_Vg_rho_SPi_rho_2_sq*regge_omega
731	(C_rho_P_cutregge_rho_P_cut
732	(cos_omega_rho_PRe_B_and_rho_2
733	+sin_omega_rho_P*Im_B_and_rho_2)
734	+ C_rho_f2_cut*regge_rho_f2_cut
735	(cos_omega_rho_f2Re_B_and_rho_2
736	+sin_omega_rho_f2*Im_B_and_rho_2));
737	double Im_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
738	=2.g_omega_Vg_rho_Vg_rho_SPi_rho_2_sq*regge_omega
739	(C_rho_P_cutregge_rho_P_cut
740	(cos_omega_rho_PIm_B_and_rho_2
741	-sin_omega_rho_P*Re_B_and_rho_2)
742	+ C_rho_f2_cut*regge_rho_f2_cut
743	(cos_omega_rho_f2Im_B_and_rho_2
744	-sin_omega_rho_f2*Re_B_and_rho_2));
745	double Re_a2_aS_omega_V_sq
746	=gsq_omega_Vg_omega_SPi_rho_2_sq*regge_omega_sq
747	(Re_B_and_rho_2one_minus_cos_omega-Im_B_and_rho_2*sin_omega);
748	double Im_a2_aS_omega_V_sq
749	=gsq_omega_Vg_omega_SPi_rho_2_sq*regge_omega_sq
750	(Re_B_and_rho_2sin_omega+Im_B_and_rho_2*one_minus_cos_omega);
751	double Re_a2_aS_omega_V_sq_omega_cuts
752	=2.gsq_omega_Vg_omega_SPi_rho_2_sqregge_omega
753	(C_omega_P_cutregge_omega_P_cut
754	(cos_omega_omega_PRe_B_and_rho_2
755	+sin_omega_omega_P*Im_B_and_rho_2)
756	+ C_omega_f2_cut*regge_omega_f2_cut
757	(cos_omega_omega_f2Re_B_and_rho_2
758	+sin_omega_omega_f2*Im_B_and_rho_2));
759	double Im_a2_aS_omega_V_sq_omega_cuts
760	=2.gsq_omega_Vg_omega_SPi_rho_2_sqregge_omega
761	(C_omega_P_cutregge_omega_P_cut
762	(cos_omega_omega_PIm_B_and_rho_2
763	-sin_omega_omega_P*Re_B_and_rho_2)
764	+ C_omega_f2_cut*regge_omega_f2_cut
765	(cos_omega_omega_f2Im_B_and_rho_2
766	-sin_omega_omega_f2*Re_B_and_rho_2));
767
768	double Re_b1_aS_rho_V_and_T
769	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_omega_1_sq
770	(Re_B_and_omega_1one_minus_cos_rho-Im_B_and_omega_1*sin_rho);
771	double Im_b1_aS_rho_V_and_T
772	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_omega_1_sq
773	(Re_B_and_omega_1sin_rho+Im_B_and_omega_1*one_minus_cos_rho);
774	double Re_b1_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_omega_1_sq
775	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_omega_1
776	-sin_rho_omega_sum2*Im_B_and_omega_1);
777	double Im_b1_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_omega_1_sq
778	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_omega_1
779	+sin_rho_omega_sum2*Re_B_and_omega_1);
780	double Re_b1_aS_rho_cuts
781	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_1_sq*regge_rho
782	(C_rho_P_cutregge_rho_P_cut
783	(cos_rho_rho_PRe_B_and_omega_1
784	+sin_rho_rho_P*Im_B_and_omega_1)
785	+ C_rho_f2_cut*regge_rho_f2_cut
786	(cos_rho_rho_f2Re_B_and_omega_1
787	+sin_rho_rho_f2*Im_B_and_omega_1));
788	double Im_b1_aS_rho_cuts
789	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_1_sq*regge_rho
790	(C_rho_P_cutregge_rho_P_cut
791	(cos_rho_rho_PIm_B_and_omega_1
792	-sin_rho_rho_P*Re_B_and_omega_1)
793	+ C_rho_f2_cut*regge_rho_f2_cut
794	(cos_rho_rho_f2Im_B_and_omega_1
795	-sin_rho_rho_f2*Re_B_and_omega_1));
796	double Re_b1_aS_omega_cuts
797	=-4.g_omega_Vg_omega_Sg_rho_TPi_omega_1_sq*regge_rho
798	(C_rho_P_cutregge_omega_P_cut
799	(cos_rho_omega_PRe_B_and_omega_1
800	+sin_rho_omega_P*Im_B_and_omega_1)
801	+ C_rho_f2_cut*regge_omega_f2_cut
802	(cos_rho_omega_f2Re_B_and_omega_1
803	+sin_rho_omega_f2*Im_B_and_omega_1));
804	double Im_b1_aS_omega_cuts
805	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_1_sq*regge_rho
806	(C_rho_P_cutregge_omega_P_cut
807	(cos_rho_omega_PIm_B_and_omega_1
808	-sin_rho_omega_P*Re_B_and_omega_1)
809	+ C_rho_f2_cut*regge_omega_f2_cut
810	(cos_rho_omega_f2Im_B_and_omega_1
811	-sin_rho_omega_f2*Re_B_and_omega_1));
812	double Re_a1_bS_rho_V_and_T
813	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_1_sq
814	(Re_B_and_omega_1one_minus_cos_rho-Im_B_and_omega_1*sin_rho);
815	double Im_a1_bS_rho_V_and_T
816	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_1_sq
817	(Re_B_and_omega_1sin_rho+Im_B_and_omega_1*one_minus_cos_rho);
818	double Re_a1_bS_omega_V
819	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_rho_1_sq
820	(cos_rho_omega_sum1Re_B_and_rho_1-sin_rho_omega_sum1*Im_B_and_rho_1);
821	double Im_a1_bS_omega_V
822	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_rho_1_sq
823	(cos_rho_omega_sum1Im_B_and_rho_1+sin_rho_omega_sum1*Re_B_and_rho_1);
824
825	double Re_b2_aS_rho_V_and_T
826	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_omega_2_sq
827	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
828	double Im_b2_aS_rho_V_and_T
829	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_omega_2_sq
830	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
831	double Re_b2_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_omega_2_sq
832	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_omega_2
833	-sin_rho_omega_sum2*Im_B_and_omega_2);
834	double Im_b2_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_omega_2_sq
835	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_omega_2
836	+sin_rho_omega_sum2*Re_B_and_omega_2);
837	double Re_b2_aS_rho_cuts
838	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_2_sq*regge_rho
839	(C_rho_P_cutregge_rho_P_cut
840	(cos_rho_rho_PRe_B_and_omega_2
841	+sin_rho_rho_P*Im_B_and_omega_2)
842	+ C_rho_f2_cut*regge_rho_f2_cut
843	(cos_rho_rho_f2Re_B_and_omega_2
844	+sin_rho_rho_f2*Im_B_and_omega_2));
845	double Im_b2_aS_rho_cuts
846	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_2_sq*regge_rho
847	(C_rho_P_cutregge_rho_P_cut
848	(cos_rho_rho_PIm_B_and_omega_2
849	-sin_rho_rho_P*Re_B_and_omega_2)
850	+ C_rho_f2_cut*regge_rho_f2_cut
851	(cos_rho_rho_f2Im_B_and_omega_2
852	-sin_rho_rho_f2*Re_B_and_omega_2));
853	double Re_b2_aS_omega_cuts
854	=-4.g_omega_Vg_omega_Sg_rho_TPi_omega_2_sq*regge_rho
855	(C_rho_P_cutregge_omega_P_cut
856	(cos_rho_omega_PRe_B_and_omega_2
857	+sin_rho_omega_P*Im_B_and_omega_2)
858	+ C_rho_f2_cut*regge_omega_f2_cut
859	(cos_rho_omega_f2Re_B_and_omega_2
860	+sin_rho_omega_f2*Im_B_and_omega_2));
861	double Im_b2_aS_omega_cuts
862	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_2_sq*regge_rho
863	(C_rho_P_cutregge_omega_P_cut
864	(cos_rho_omega_PIm_B_and_omega_2
865	-sin_rho_omega_P*Re_B_and_omega_2)
866	+ C_rho_f2_cut*regge_omega_f2_cut
867	(cos_rho_omega_f2Im_B_and_omega_2
868	-sin_rho_omega_f2*Re_B_and_omega_2));
869
870	double Re_a2_bS_rho_V_and_T
871	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
872	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
873	double Im_a2_bS_rho_V_and_T
874	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
875	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
876	double Re_a2_bS_omega_V
877	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_rho_2_sq
878	(cos_rho_omega_sum1Re_B_and_rho_2-sin_rho_omega_sum1*Im_B_and_rho_2);
879	double Im_a2_bS_omega_V
880	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_rho_2_sq
881	(cos_rho_omega_sum1Im_B_and_rho_2+sin_rho_omega_sum1*Re_B_and_rho_2);
882
883	Re_b1_bS=4.CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_omega_1_sq
884	(Re_B_and_omega_1one_minus_cos_rho-Im_B_and_omega_1*sin_rho);
885	Im_b1_bS=8.CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_omega_1_sq
886	(Re_B_and_omega_1sin_rho+Im_B_and_omega_1*one_minus_cos_rho);
887	Re_b2_bS=4.CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_omega_2_sq
888	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
889	Im_b2_bS=4.CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_omega_2_sq
890	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
891
892	Re_b2_aS=CgR(Re_b2_aS_rho_V_and_T + Re_b2_aS_rho_cuts
893	+ Re_b2_aS_omega_V + Re_b2_aS_omega_cuts
894	);
895	Im_b2_aS=CgR(Im_b2_aS_rho_V_and_T + Im_b2_aS_rho_cuts
896	+ Im_b2_aS_omega_V + Im_b2_aS_omega_cuts
897	);
898	Re_b1_aS=CgR(Re_b1_aS_rho_V_and_T + Re_b1_aS_rho_cuts
899	+ Re_b1_aS_omega_V + Re_b1_aS_omega_cuts
900	);
901	Im_b1_aS=CgR(Im_b1_aS_rho_V_and_T + Im_b1_aS_rho_cuts
902	+ Im_b1_aS_omega_V + Im_b1_aS_omega_cuts
903	);
904
905	Re_a1_bS=CgR(Re_a1_bS_rho_V_and_T + (1./3.)*Re_a1_bS_omega_V);
906	Im_a1_bS=CgR(Im_a1_bS_rho_V_and_T + (1./3.)*Im_a1_bS_omega_V);
907	Re_a2_bS=CgR(Re_a2_bS_rho_V_and_T + (1./3.)*Re_a2_bS_omega_V);
908	Im_a2_bS=CgR(Im_a2_bS_rho_V_and_T + (1./3.)*Im_a2_bS_omega_V);
909
910	Re_a1_aS=CgR(Re_a1_aS_rho_V_and_T_sq
911	+ Re_a1_aS_rho_V_and_T_rho_cuts
912	+ Re_a1_aS_omega_V_rho_V_and_T
913	+ Re_a1_aS_rho_V_and_T_omega_cuts
914	+ (1./3.)*Re_a1_aS_omega_V_rho_V_and_T_rho_propagator
915	+ (1./3.)*Re_a1_aS_omega_V_rho_V_rho_cuts_rho_propagator
916	+ (1./3.)*Re_a1_aS_omega_V_sq
917	+ (1./3.)*Re_a1_aS_omega_V_sq_omega_cuts
918	);
919	Im_a1_aS=CgR(Im_a1_aS_rho_V_and_T_sq
920	+ Im_a1_aS_rho_V_and_T_rho_cuts
921	+ Im_a1_aS_omega_V_rho_V_and_T
922	+ Im_a1_aS_rho_V_and_T_omega_cuts
923	+ (1./3.)*Im_a1_aS_omega_V_rho_V_and_T_rho_propagator
924	+ (1./3.)*Im_a1_aS_omega_V_rho_V_rho_cuts_rho_propagator
925	+ (1./3.)*Im_a1_aS_omega_V_sq
926	+ (1./3.)*Im_a1_aS_omega_V_sq_omega_cuts
927	);
928
929	Re_a2_aS=CgR(Re_a2_aS_rho_V_and_T_sq
930	+ Re_a2_aS_rho_V_and_T_rho_cuts
931	+ Re_a2_aS_omega_V_rho_V_and_T
932	+ Re_a2_aS_rho_V_and_T_omega_cuts
933	+ (1./3.)*Re_a2_aS_omega_V_rho_V_and_T_rho_propagator
934	+ (1./3.)*Re_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
935	+ (1./3.)*Re_a2_aS_omega_V_sq
936	+ (1./3.)*Re_a2_aS_omega_V_sq_omega_cuts
937	);
938	Im_a2_aS=CgR(Im_a2_aS_rho_V_and_T_sq
939	+ Im_a2_aS_rho_V_and_T_rho_cuts
940	+ Im_a2_aS_omega_V_rho_V_and_T
941	+ Im_a2_aS_rho_V_and_T_omega_cuts
942	+ (1./3.)*Im_a2_aS_omega_V_rho_V_and_T_rho_propagator
943	+ (1./3.)*Im_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
944	+ (1./3.)*Im_a2_aS_omega_V_sq
945	+ (1./3.)*Im_a2_aS_omega_V_sq_omega_cuts
946	);
947	}
948	if (two_particles==(7+17)){ // pi0 eta
949	C=sqrt(2./3.);
950
951	double Re_a1_aS_rho_V_and_T_sq
952	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_rho_1_sq
953	(Re_B_and_rho_1one_minus_cos_rho-Im_B_and_rho_1*sin_rho);
954	double Im_a1_aS_rho_V_and_T_sq
955	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_rho_1_sq
956	(Re_B_and_rho_1sin_rho+Im_B_and_rho_1*one_minus_cos_rho);
957	double Re_a1_aS_rho_V_and_T_rho_cuts
958	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_rho_1_sq
959	(C_rho_P_cutregge_rho_P_cut
960	(cos_rho_rho_PRe_B_and_rho_1
961	+sin_rho_rho_P*Im_B_and_rho_1)
962	+ C_rho_f2_cut*regge_rho_f2_cut
963	(cos_rho_rho_f2Re_B_and_rho_1
964	+sin_rho_rho_f2*Im_B_and_rho_1));
965	double Im_a1_aS_rho_V_and_T_rho_cuts
966	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_rho_1_sq
967	(C_rho_P_cutregge_rho_P_cut
968	(cos_rho_rho_PIm_B_and_rho_1
969	-sin_rho_rho_P*Re_B_and_rho_1)
970	+ C_rho_f2_cut*regge_rho_f2_cut
971	(cos_rho_rho_f2Im_B_and_rho_1
972	-sin_rho_rho_f2*Re_B_and_rho_1));
973	double Re_a1_aS_omega_V_rho_V_and_T
974	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_rho_1_sq
975	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_rho_1
976	-sin_rho_omega_sum2*Im_B_and_rho_1);
977	double Im_a1_aS_omega_V_rho_V_and_T
978	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_rho_1_sq
979	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_rho_1
980	+sin_rho_omega_sum2*Re_B_and_rho_1);
981	double Re_a1_aS_rho_V_and_T_omega_cuts
982	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_rho_1_sq
983	(C_omega_P_cutregge_omega_P_cut
984	(cos_rho_omega_PRe_B_and_rho_1
985	+sin_rho_omega_P*Im_B_and_rho_1)
986	+ C_omega_f2_cut*regge_omega_f2_cut
987	(cos_rho_omega_f2Re_B_and_rho_1
988	+sin_rho_omega_f2*Im_B_and_rho_1));
989	double Im_a1_aS_rho_V_and_T_omega_cuts
990	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_rho_1_sq
991	(C_omega_P_cutregge_omega_P_cut
992	(cos_rho_omega_PIm_B_and_rho_1
993	-sin_rho_omega_P*Re_B_and_rho_1)
994	+ C_rho_f2_cut*regge_omega_f2_cut
995	(cos_rho_omega_f2Im_B_and_rho_1
996	-sin_rho_omega_f2*Re_B_and_rho_1));
997	double Re_a1_aS_omega_V_rho_V_and_T_omega_propagator
998	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_omega_1_sq
999	regge_rhoregge_omega(cos_rho_omega_sum1Re_B_and_omega_1
1000	-sin_rho_omega_sum1*Im_B_and_omega_1);
1001	double Im_a1_aS_omega_V_rho_V_and_T_omega_propagator
1002	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_omega_1_sq
1003	regge_rhoregge_omega(cos_rho_omega_sum1Im_B_and_omega_1
1004	+sin_rho_omega_sum1*Re_B_and_omega_1);
1005	double Re_a1_aS_omega_V_rho_V_rho_cuts_omega_propagator
1006	=2.g_omega_Vg_rho_Vg_rho_SPi_omega_1_sq*regge_omega
1007	(C_rho_P_cutregge_rho_P_cut
1008	(cos_omega_rho_PRe_B_and_omega_1
1009	+sin_omega_rho_P*Im_B_and_omega_1)
1010	+ C_rho_f2_cut*regge_rho_f2_cut
1011	(cos_omega_rho_f2Re_B_and_omega_1
1012	+sin_omega_rho_f2*Im_B_and_omega_1));
1013	double Im_a1_aS_omega_V_rho_V_rho_cuts_omega_propagator
1014	=2.g_omega_Vg_rho_Vg_rho_SPi_omega_1_sq*regge_omega
1015	(C_rho_P_cutregge_rho_P_cut
1016	(cos_omega_rho_PIm_B_and_omega_1
1017	-sin_omega_rho_P*Re_B_and_omega_1)
1018	+ C_rho_f2_cut*regge_rho_f2_cut
1019	(cos_omega_rho_f2Im_B_and_omega_1
1020	-sin_omega_rho_f2*Re_B_and_omega_1));
1021	double Re_a1_aS_omega_V_sq
1022	=g_omega_Vg_omega_Vg_omega_SPi_omega_1_sqregge_omega_sq
1023	(Re_B_and_omega_1one_minus_cos_omega-Im_B_and_omega_1*sin_omega);
1024	double Im_a1_aS_omega_V_sq
1025	=g_omega_Vg_omega_Vg_omega_SPi_omega_1_sqregge_omega_sq
1026	(Re_B_and_omega_1sin_omega+Im_B_and_omega_1*one_minus_cos_omega);
1027	double Re_a1_aS_omega_V_sq_omega_cuts
1028	=2.g_omega_Vg_omega_Vg_omega_SPi_omega_1_sq*regge_omega
1029	(C_omega_P_cutregge_omega_P_cut
1030	(cos_omega_omega_PRe_B_and_omega_1
1031	+sin_omega_omega_P*Im_B_and_omega_1)
1032	+ C_omega_f2_cut*regge_omega_f2_cut
1033	(cos_omega_omega_f2Re_B_and_omega_1
1034	+sin_omega_omega_f2*Im_B_and_omega_1));
1035	double Im_a1_aS_omega_V_sq_omega_cuts
1036	=2.g_omega_Vg_omega_Vg_omega_SPi_omega_1_sq*regge_omega
1037	(C_omega_P_cutregge_omega_P_cut
1038	(cos_omega_omega_PIm_B_and_omega_1
1039	-sin_omega_omega_P*Re_B_and_omega_1)
1040	+ C_omega_f2_cut*regge_omega_f2_cut
1041	(cos_omega_omega_f2Im_B_and_omega_1
1042	-sin_omega_omega_f2*Re_B_and_omega_1));
1043
1044	double Re_a2_aS_rho_V_and_T_sq
1045	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
1046	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
1047	double Im_a2_aS_rho_V_and_T_sq
1048	=gsq_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
1049	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
1050	double Re_a2_aS_rho_V_and_T_rho_cuts
1051	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_omega_2_sq
1052	(C_rho_P_cutregge_rho_P_cut
1053	(cos_rho_rho_PRe_B_and_omega_2
1054	+sin_rho_rho_P*Im_B_and_omega_2)
1055	+ C_rho_f2_cut*regge_rho_f2_cut
1056	(cos_rho_rho_f2Re_B_and_omega_2
1057	+sin_rho_rho_f2*Im_B_and_omega_2));
1058	double Im_a2_aS_rho_V_and_T_rho_cuts
1059	=2.g_rho_Sg_rho_V_and_Tg_rho_Vregge_rho*Pi_omega_2_sq
1060	(C_rho_P_cutregge_rho_P_cut
1061	(cos_rho_rho_PIm_B_and_omega_2
1062	-sin_rho_rho_P*Re_B_and_omega_2)
1063	+ C_rho_f2_cut*regge_rho_f2_cut
1064	(cos_rho_rho_f2Im_B_and_omega_2
1065	-sin_rho_rho_f2*Re_B_and_omega_2));
1066	double Re_a2_aS_omega_V_rho_V_and_T
1067	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_omega_2_sq
1068	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_omega_2
1069	-sin_rho_omega_sum2*Im_B_and_omega_2);
1070	double Im_a2_aS_omega_V_rho_V_and_T
1071	=g_omega_Sg_rho_V_and_Tg_omega_V*Pi_omega_2_sq
1072	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_omega_2
1073	+sin_rho_omega_sum2*Re_B_and_omega_2);
1074	double Re_a2_aS_rho_V_and_T_omega_cuts
1075	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_omega_2_sq
1076	(C_omega_P_cutregge_omega_P_cut
1077	(cos_rho_omega_PRe_B_and_omega_2
1078	+sin_rho_omega_P*Im_B_and_omega_2)
1079	+ C_omega_f2_cut*regge_omega_f2_cut
1080	(cos_rho_omega_f2Re_B_and_omega_2
1081	+sin_rho_omega_f2*Im_B_and_omega_2));
1082	double Im_a2_aS_rho_V_and_T_omega_cuts
1083	=2.g_omega_Sg_rho_V_and_Tg_omega_Vregge_rho*Pi_omega_2_sq
1084	(C_omega_P_cutregge_omega_P_cut
1085	(cos_rho_omega_PIm_B_and_omega_2
1086	-sin_rho_omega_P*Re_B_and_omega_2)
1087	+ C_rho_f2_cut*regge_omega_f2_cut
1088	(cos_rho_omega_f2Im_B_and_omega_2
1089	-sin_rho_omega_f2*Re_B_and_omega_2));
1090	double Re_a2_aS_omega_V_rho_V_and_T_rho_propagator
1091	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_rho_2_sq
1092	regge_rhoregge_omega(cos_rho_omega_sum1Re_B_and_rho_2
1093	-sin_rho_omega_sum1*Im_B_and_rho_2);
1094	double Im_a2_aS_omega_V_rho_V_and_T_rho_propagator
1095	=g_rho_Sg_rho_V_and_Tg_omega_V*Pi_rho_2_sq
1096	regge_rhoregge_omega(cos_rho_omega_sum1Im_B_and_rho_2
1097	+sin_rho_omega_sum1*Re_B_and_rho_2);
1098	double Re_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
1099	=2.g_omega_Vg_rho_Vg_rho_SPi_rho_2_sq*regge_omega
1100	(C_rho_P_cutregge_rho_P_cut
1101	(cos_omega_rho_PRe_B_and_rho_2
1102	+sin_omega_rho_P*Im_B_and_rho_2)
1103	+ C_rho_f2_cut*regge_rho_f2_cut
1104	(cos_omega_rho_f2Re_B_and_rho_2
1105	+sin_omega_rho_f2*Im_B_and_rho_2));
1106	double Im_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
1107	=2.g_omega_Vg_rho_Vg_rho_SPi_rho_2_sq*regge_omega
1108	(C_rho_P_cutregge_rho_P_cut
1109	(cos_omega_rho_PIm_B_and_rho_2
1110	-sin_omega_rho_P*Re_B_and_rho_2)
1111	+ C_rho_f2_cut*regge_rho_f2_cut
1112	(cos_omega_rho_f2Im_B_and_rho_2
1113	-sin_omega_rho_f2*Re_B_and_rho_2));
1114	double Re_a2_aS_omega_V_sq
1115	=g_omega_Vg_omega_Vg_omega_SPi_rho_2_sqregge_omega_sq
1116	(Re_B_and_rho_2one_minus_cos_omega-Im_B_and_rho_2*sin_omega);
1117	double Im_a2_aS_omega_V_sq
1118	=g_omega_Vg_omega_Vg_omega_SPi_rho_2_sqregge_omega_sq
1119	(Re_B_and_rho_2sin_omega+Im_B_and_rho_2*one_minus_cos_omega);
1120	double Re_a2_aS_omega_V_sq_omega_cuts
1121	=2.g_omega_Vg_omega_Vg_omega_SPi_rho_2_sq*regge_omega
1122	(C_omega_P_cutregge_omega_P_cut
1123	(cos_omega_omega_PRe_B_and_rho_2
1124	+sin_omega_omega_P*Im_B_and_rho_2)
1125	+ C_omega_f2_cut*regge_omega_f2_cut
1126	(cos_omega_omega_f2Re_B_and_rho_2
1127	+sin_omega_omega_f2*Im_B_and_rho_2));
1128	double Im_a2_aS_omega_V_sq_omega_cuts
1129	=2.g_omega_Vg_omega_Vg_omega_SPi_rho_2_sq*regge_omega
1130	(C_omega_P_cutregge_omega_P_cut
1131	(cos_omega_omega_PIm_B_and_rho_2
1132	-sin_omega_omega_P*Re_B_and_rho_2)
1133	+ C_omega_f2_cut*regge_omega_f2_cut
1134	(cos_omega_omega_f2Im_B_and_rho_2
1135	-sin_omega_omega_f2*Re_B_and_rho_2));
1136
1137	double Re_b1_aS_rho_V_and_T
1138	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_rho_1_sq
1139	(Re_B_and_rho_1one_minus_cos_rho-Im_B_and_rho_1*sin_rho);
1140	double Im_b1_aS_rho_V_and_T
1141	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_rho_1_sq
1142	(Re_B_and_rho_1sin_rho+Im_B_and_rho_1*one_minus_cos_rho);
1143	double Re_b1_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_rho_1_sq
1144	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_rho_1
1145	-sin_rho_omega_sum2*Im_B_and_rho_1);
1146	double Im_b1_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_rho_1_sq
1147	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_rho_1
1148	+sin_rho_omega_sum2*Re_B_and_rho_1);
1149	double Re_b1_aS_rho_cuts
1150	=-4.g_rho_Vg_rho_Sg_rho_TPi_rho_1_sq*regge_rho
1151	(C_rho_P_cutregge_rho_P_cut
1152	(cos_rho_rho_PRe_B_and_rho_1
1153	+sin_rho_rho_P*Im_B_and_rho_1)
1154	+ C_rho_f2_cut*regge_rho_f2_cut
1155	(cos_rho_rho_f2Re_B_and_rho_1
1156	+sin_rho_rho_f2*Im_B_and_rho_1));
1157	double Im_b1_aS_rho_cuts
1158	=-4.g_rho_Vg_rho_Sg_rho_TPi_rho_1_sq*regge_rho
1159	(C_rho_P_cutregge_rho_P_cut
1160	(cos_rho_rho_PIm_B_and_rho_1
1161	-sin_rho_rho_P*Re_B_and_rho_1)
1162	+ C_rho_f2_cut*regge_rho_f2_cut
1163	(cos_rho_rho_f2Im_B_and_rho_1
1164	-sin_rho_rho_f2*Re_B_and_rho_1));
1165	double Re_b1_aS_omega_cuts
1166	=-4.g_omega_Vg_omega_Sg_rho_TPi_rho_1_sq*regge_rho
1167	(C_rho_P_cutregge_omega_P_cut
1168	(cos_rho_omega_PRe_B_and_rho_1
1169	+sin_rho_omega_P*Im_B_and_rho_1)
1170	+ C_rho_f2_cut*regge_omega_f2_cut
1171	(cos_rho_omega_f2Re_B_and_rho_1
1172	+sin_rho_omega_f2*Im_B_and_rho_1));
1173	double Im_b1_aS_omega_cuts
1174	=-4.g_rho_Vg_rho_Sg_rho_TPi_rho_1_sq*regge_rho
1175	(C_rho_P_cutregge_omega_P_cut
1176	(cos_rho_omega_PIm_B_and_rho_1
1177	-sin_rho_omega_P*Re_B_and_rho_1)
1178	+ C_rho_f2_cut*regge_omega_f2_cut
1179	(cos_rho_omega_f2Im_B_and_rho_1
1180	-sin_rho_omega_f2*Re_B_and_rho_1));
1181	double Re_b2_aS_rho_V_and_T
1182	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_omega_2_sq
1183	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
1184	double Im_b2_aS_rho_V_and_T
1185	=-2.g_rho_Sg_rho_Tg_rho_V_and_Tregge_rho_sq*Pi_omega_2_sq
1186	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
1187	double Re_b2_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_omega_2_sq
1188	regge_rhoregge_omega(cos_rho_omega_sum2Re_B_and_omega_2
1189	-sin_rho_omega_sum2*Im_B_and_omega_2);
1190	double Im_b2_aS_omega_V=-2.g_omega_Sg_rho_Vg_rho_TPi_omega_2_sq
1191	regge_rhoregge_omega(cos_rho_omega_sum2Im_B_and_omega_2
1192	+sin_rho_omega_sum2*Re_B_and_omega_2);
1193	double Re_b2_aS_rho_cuts
1194	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_2_sq*regge_rho
1195	(C_rho_P_cutregge_rho_P_cut
1196	(cos_rho_rho_PRe_B_and_omega_2
1197	+sin_rho_rho_P*Im_B_and_omega_2)
1198	+ C_rho_f2_cut*regge_rho_f2_cut
1199	(cos_rho_rho_f2Re_B_and_omega_2
1200	+sin_rho_rho_f2*Im_B_and_omega_2));
1201	double Im_b2_aS_rho_cuts
1202	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_2_sq*regge_rho
1203	(C_rho_P_cutregge_rho_P_cut
1204	(cos_rho_rho_PIm_B_and_omega_2
1205	-sin_rho_rho_P*Re_B_and_omega_2)
1206	+ C_rho_f2_cut*regge_rho_f2_cut
1207	(cos_rho_rho_f2Im_B_and_omega_2
1208	-sin_rho_rho_f2*Re_B_and_omega_2));
1209	double Re_b2_aS_omega_cuts
1210	=-4.g_omega_Vg_omega_Sg_rho_TPi_omega_2_sq*regge_rho
1211	(C_rho_P_cutregge_omega_P_cut
1212	(cos_rho_omega_PRe_B_and_omega_2
1213	+sin_rho_omega_P*Im_B_and_omega_2)
1214	+ C_rho_f2_cut*regge_omega_f2_cut
1215	(cos_rho_omega_f2Re_B_and_omega_2
1216	+sin_rho_omega_f2*Im_B_and_omega_2));
1217	double Im_b2_aS_omega_cuts
1218	=-4.g_rho_Vg_rho_Sg_rho_TPi_omega_2_sq*regge_rho
1219	(C_rho_P_cutregge_omega_P_cut
1220	(cos_rho_omega_PIm_B_and_omega_2
1221	-sin_rho_omega_P*Re_B_and_omega_2)
1222	+ C_rho_f2_cut*regge_omega_f2_cut
1223	(cos_rho_omega_f2Im_B_and_omega_2
1224	-sin_rho_omega_f2*Re_B_and_omega_2));
1225
1226	double Re_a1_bS_rho_V_and_T
1227	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_rho_1_sq
1228	(Re_B_and_rho_1one_minus_cos_rho-Im_B_and_rho_1*sin_rho);
1229	double Im_a1_bS_rho_V_and_T
1230	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_rho_1_sq
1231	(Re_B_and_rho_1sin_rho+Im_B_and_rho_1*one_minus_cos_rho);
1232	double Re_a1_bS_omega_V
1233	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_omega_1_sq
1234	(cos_rho_omega_sum1Re_B_and_omega_1-sin_rho_omega_sum1*Im_B_and_omega_1);
1235	double Im_a1_bS_omega_V
1236	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_omega_1_sq
1237	(cos_rho_omega_sum1Im_B_and_omega_1+sin_rho_omega_sum1*Re_B_and_omega_1);
1238
1239	double Re_a2_bS_rho_V_and_T
1240	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
1241	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
1242	double Im_a2_bS_rho_V_and_T
1243	=-2.g_rho_Tg_rho_V_and_Tg_rho_Sregge_rho_sq*Pi_omega_2_sq
1244	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
1245	double Re_a2_bS_omega_V
1246	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_rho_2_sq
1247	(cos_rho_omega_sum1Re_B_and_rho_2-sin_rho_omega_sum1*Im_B_and_rho_2);
1248	double Im_a2_bS_omega_V
1249	=-2.g_rho_Sg_rho_Vg_rho_Tregge_rhoregge_omegaPi_rho_2_sq
1250	(cos_rho_omega_sum1Im_B_and_rho_2+sin_rho_omega_sum1*Re_B_and_rho_2);
1251
1252
1253	Re_b1_bS=(4./3.)CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_rho_1_sq
1254	(Re_B_and_rho_1one_minus_cos_rho-Im_B_and_rho_1*sin_rho);
1255	Im_b1_bS=(4./3.)CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_rho_1_sq
1256	(Re_B_and_rho_1sin_rho+Im_B_and_rho_1*one_minus_cos_rho);
1257	Re_b2_bS=(4./3.)CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_omega_2_sq
1258	(Re_B_and_omega_2one_minus_cos_rho-Im_B_and_omega_2*sin_rho);
1259	Im_b2_bS=(4./3.)CgRgsq_rho_Tg_rho_Sregge_rho_sqPi_omega_2_sq
1260	(Re_B_and_omega_2sin_rho+Im_B_and_omega_2*one_minus_cos_rho);
1261
1262	Re_b2_aS=(1./3.)CgR*(Re_b2_aS_rho_V_and_T + Re_b2_aS_rho_cuts
1263	+ Re_b2_aS_omega_V + Re_b2_aS_omega_cuts
1264	);
1265	Im_b2_aS=(1./3.)CgR*(Im_b2_aS_rho_V_and_T + Im_b2_aS_rho_cuts
1266	+ Im_b2_aS_omega_V + Im_b2_aS_omega_cuts
1267	);
1268	Re_b1_aS=(1./3.)CgR*(Re_b1_aS_rho_V_and_T + Re_b1_aS_rho_cuts
1269	+ Re_b1_aS_omega_V + Re_b1_aS_omega_cuts
1270	);
1271	Im_b1_aS=(1./3.)CgR*(Im_b1_aS_rho_V_and_T + Im_b1_aS_rho_cuts
1272	+ Im_b1_aS_omega_V + Im_b1_aS_omega_cuts
1273	);
1274
1275	Re_a1_bS=CgR((1./3.)*Re_a1_bS_rho_V_and_T + Re_a1_bS_omega_V);
1276	Im_a1_bS=CgR((1./3.)*Im_a1_bS_rho_V_and_T + Im_a1_bS_omega_V);
1277	Re_a2_bS=CgR((1./3.)*Re_a2_bS_rho_V_and_T + Re_a2_bS_omega_V);
1278	Im_a2_bS=CgR((1./3.)*Im_a2_bS_rho_V_and_T + Im_a2_bS_omega_V);
1279
1280	Re_a1_aS=CgR((1./3.)*Re_a1_aS_rho_V_and_T_sq
1281	+ (1./3.)*Re_a1_aS_rho_V_and_T_rho_cuts
1282	+ (1./3.)*Re_a1_aS_omega_V_rho_V_and_T
1283	+ (1./3.)*Re_a1_aS_rho_V_and_T_omega_cuts
1284	+ Re_a1_aS_omega_V_rho_V_and_T_omega_propagator
1285	+ Re_a1_aS_omega_V_rho_V_rho_cuts_omega_propagator
1286	+ Re_a1_aS_omega_V_sq
1287	+ Re_a1_aS_omega_V_sq_omega_cuts
1288	);
1289	Im_a1_aS=CgR((1./3.)*Im_a1_aS_rho_V_and_T_sq
1290	+ (1./3.)*Im_a1_aS_rho_V_and_T_rho_cuts
1291	+ (1./3.)*Im_a1_aS_omega_V_rho_V_and_T
1292	+ (1./3.)*Im_a1_aS_rho_V_and_T_omega_cuts
1293	+ Im_a1_aS_omega_V_rho_V_and_T_omega_propagator
1294	+ Im_a1_aS_omega_V_rho_V_rho_cuts_omega_propagator
1295	+ Im_a1_aS_omega_V_sq
1296	+ Im_a1_aS_omega_V_sq_omega_cuts
1297	);
1298
1299	Re_a2_aS=CgR((1./3.)*Re_a2_aS_rho_V_and_T_sq
1300	+ (1./3.)*Re_a2_aS_rho_V_and_T_rho_cuts
1301	+ (1./3.)*Re_a2_aS_omega_V_rho_V_and_T
1302	+ (1./3.)*Re_a2_aS_rho_V_and_T_omega_cuts
1303	+ Re_a2_aS_omega_V_rho_V_and_T_rho_propagator
1304	+ Re_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
1305	+ Re_a2_aS_omega_V_sq
1306	+ Re_a2_aS_omega_V_sq_omega_cuts
1307	);
1308	Im_a2_aS=CgR((1./3.)*Im_a2_aS_rho_V_and_T_sq
1309	+ (1./3.)*Im_a2_aS_rho_V_and_T_rho_cuts
1310	+ (1./3.)*Im_a2_aS_omega_V_rho_V_and_T
1311	+ (1./3.)*Im_a2_aS_rho_V_and_T_omega_cuts
1312	+ Im_a2_aS_omega_V_rho_V_and_T_rho_propagator
1313	+ Im_a2_aS_omega_V_rho_V_rho_cuts_rho_propagator
1314	+ Im_a2_aS_omega_V_sq
1315	+ Im_a2_aS_omega_V_sq_omega_cuts
1316	);
1317
1318
1319	}
1320
1321	// Real and imaginary parts of the trace term
1322	double ReT=(m_p_sq-p1_dot_p2)(Re_a1_aSM_M1+Re_a2_aS*M_M2
1323	+2.Re_a1_aSN_N1+2.Re_a2_aSN_N2)
1324	+2.m_p((Re_b1_aS+Re_a1_bS)N_N1+(Re_b2_aS+Re_a2_bS)N_N2)
1325	+(m_p_sq+p1_dot_p2)(Re_b1_bSN_N1+Re_b2_bS*N_N2);
1326	double ImT=(m_p_sq-p1_dot_p2)(Im_a1_aSM_M1+Im_a2_aS*M_M2
1327	+2.Im_a1_aSN_N1+2.Im_a2_aSN_N2)
1328	+2.m_p((Im_b1_aS+Im_a1_bS)N_N1+(Im_b2_aS+Im_a2_bS)N_N2)
1329	+(m_p_sq+p1_dot_p2)(Im_b1_bSN_N1+Im_b2_bS*N_N2);
1330
1331	double Msq=(particles[0]+particles[1]).M2();
1332	double g0=sqrt(g0_sq);
1333	double s_minus_mp_sq=s-m_p_sq;
1334	double m1=ParticleMass(particle_types[0]);
1335	double m2=ParticleMass(particle_types[1]);
1336	double m1_plus_m2=m1+m2;
1337	double m1_minus_m2=m1-m2;
1338	double k=sqrt((Msq-m1_plus_m2m1_plus_m2)(Msq-m1_minus_m2*m1_minus_m2))
1339	/(2.*sqrt(Msq));
1340	double hbarc_sq=389.; // Convert to micro-barns
1341	double xsec=-zetagRg0khbarc_sq(cos(phase)ReT+sin(phase)*ImT)
1342	/(256M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846s_minus_mp_sqs_minus_mp_sq(ReBReB+ImB*ImB));
1343
1344	return xsec;
1345	}
1346
1347
1348	// Get parameters for Breit-Wigner distribution for a0(980)/f0(980) resonance shape
1349	void GetResonanceParameters(double m1,double m2, double M_sq,double M_sq_R,
1350	double &ReB,double &ImB){
1351	double m1_plus_m2=m1+m2;
1352	double m1_minus_m2=m1-m2;
1353	bool got_pipi=(fabs(m1_minus_m2)>0.01)?false:true;
1354
1355	// "No structure" model for a0(980)/f0(980)
1356	// masses
1357	double M_S=sqrt(M_sq);
1358	double MK0=ParticleMass(KShort);
1359	double MKPlus=ParticleMass(KPlus);
1360
1361	// coupling constants
1362	double gK=3.05;
1363	double g_m1m2=2.82;
1364	if (got_pipi){
1365	gK=5.006;
1366	g_m1m2=1.70315;
1367	}
1368	double gKsq=gK*gK;
1369	double g_m1m2_sq=g_m1m2*g_m1m2;
1370
1371	// kinematic factors
1372	double rhoK0sq=1.-4.MK0MK0/M_sq;
1373	double rhoKPlussq=1.-4.MKPlusMKPlus/M_sq;
1374	double rho_m1m2
1375	=sqrt((1.-m1_plus_m2m1_plus_m2/M_sq)(1-m1_minus_m2*m1_minus_m2/M_sq));
1376
1377	// Real and imaginary parts of BW amplitude
1378	ReB=M_sq_R-M_sq;
1379	if (M_S<2.*MKPlus){
1380	ReB+=gKsq/(32.M_PI3.14159265358979323846)sqrt(-rhoKPlussq);
1381	}
1382	if (M_S<2.*MK0){
1383	ReB+=gKsq/(32.M_PI3.14159265358979323846)sqrt(-rhoK0sq);
1384	}
1385	ImB=g_m1m2_sq/(16M_PI3.14159265358979323846)rho_m1m2;
1386	if (M_S>2.*MKPlus){
1387	ImB+=gKsq/(32.M_PI3.14159265358979323846)sqrt(rhoKPlussq);
1388	}
1389	if (M_S>2.*MK0){
1390	ImB+=gKsq/(32.M_PI3.14159265358979323846)sqrt(rhoK0sq);
1391	}
1392	}
1393
1394
1395	// Cross section dsigma/(dt/dM/dOmega) from Donnachie and Kalashnikova
1396	double CrossSection(double m1,double m2, double ms_sq, double s, double t,
1397	double gR,double ReB,double ImB,double gsq_rho_S,
1398	double gsq_omega_S, bool interfere=false,
1399	double gR2=0.,double ReB2=0.,double ImB2=0.,
1400	double gsq_rho_S2=0.,
1401	double gsq_omega_S2=0.,double phase=0.){
1402	// Kinematic factors
1403	double mp_sq_minus_s=m_p_sq-s;
1404	double mp_sq_plus_s=m_p_sq+s;
1405	double mp_sq_minus_s_sq=mp_sq_minus_s*mp_sq_minus_s;
1406	double temp1=ms_sq*mp_sq_plus_s-mp_sq_minus_s_sq;
1407	double temp2=mp_sq_minus_ssqrt(mp_sq_minus_s_sq-2.ms_sq*mp_sq_plus_s
1408	+ms_sq*ms_sq);
1409	double t1=(temp1+temp2)/(2.*s);
1410	double t2=(temp1-temp2)/(2.*s);
1411	double t_minus_ms_sq=t-ms_sq;
1412	double kin1=s(t-t1)(t-t2);
1413	double kin_aS_aS=0.5kin1+0.25tt_minus_ms_sqt_minus_ms_sq;
1414	double kin_aS_bS=0.5m_pkin1;
1415	double kin_bS_bS=0.125(4.m_p_sq-t)*kin1;
1416
1417	// Coupling constants
1418	double g_omega_V=15.;
1419	double gsq_omega_V=g_omega_V*g_omega_V;
1420	double g_rho_V=3.4;
1421	double gsq_rho_V=g_rho_V*g_rho_V;
1422	double g_rho_T=11.0; // GeV^-1
1423	double g_rho_V_and_T=g_rho_V+2.m_pg_rho_T;
1424	double gsq_rho_V_and_T=g_rho_V_and_T*g_rho_V_and_T;
1425	double g_rho_S=sqrt(gsq_rho_S);
1426	double g_omega_S=sqrt(gsq_omega_S);
1427	double gRsq=gR*gR;
1428
1429	// s scale for regge trajectories
1430	double s0=1.;
1431
1432	// Regge trajectory for omega
1433	double a_omega=0.44+0.9*t;
1434	double a_omega_prime=0.9;
1435	double cos_omega=cos(M_PI3.14159265358979323846*a_omega);
1436	double regge_omega=pow(s/s0,a_omega-1.)M_PI3.14159265358979323846a_omega_prime/(sin(M_PI3.14159265358979323846a_omega)TMath::Gamma(a_omega)); // excluding phase factor
1437	double regge_omega_sq=regge_omegaregge_omega0.5*(1.-cos_omega);
1438
1439	// Regge cuts for omega
1440	double a_omega_P=0.52+0.196*t; // Pomeron
1441	double a_omega_f2=0.112+0.428*t;
1442	double dc=regge_cuts[0];
1443	double regge_omega_P_cut=exp(dct)pow(s/s0,a_omega_P-1.);
1444	double regge_omega_f2_cut=exp(dct)pow(s/s0,a_omega_f2-1.);
1445	double C_omega_P_cut=regge_cuts[1];
1446	double C_omega_f2_cut=regge_cuts[2];
1447	double Csq_omega_P_cut=C_omega_P_cut*C_omega_P_cut;
1448	double Csq_omega_f2_cut=C_omega_f2_cut*C_omega_f2_cut;
1449	double C_omega_f2_C_omega_P=C_omega_f2_cut*C_omega_P_cut;
1450
1451	// Regge trajectory for rho
1452	double a_rho=0.55+0.8*t;
1453	double a_rho_prime=0.8;
1454	double cos_rho=cos(M_PI3.14159265358979323846*a_rho);
1455	double regge_rho=pow(s/s0,a_rho-1.)M_PI3.14159265358979323846a_rho_prime/(sin(M_PI3.14159265358979323846a_rho)TMath::Gamma(a_rho)); // excluding phase factor
1456	double regge_rho_sq=regge_rhoregge_rho0.5*(1.-cos_rho);
1457
1458	// Regge cuts for rho
1459	double a_rho_P=0.64+0.16*t; // Pomeron
1460	double a_rho_f2=0.222+0.404*t;
1461	double regge_rho_P_cut=exp(dct)pow(s/s0,a_rho_P-1.);
1462	double regge_rho_f2_cut=exp(dct)pow(s/s0,a_rho_f2-1.);
1463	double C_rho_P_cut=regge_cuts[3];
1464	double C_rho_f2_cut=regge_cuts[4];
1465	double Csq_rho_P_cut=C_rho_P_cut*C_rho_P_cut;
1466	double Csq_rho_f2_cut=C_rho_f2_cut*C_rho_f2_cut;
1467
1468	// rho-omega interference
1469	double cos_rho_omega_sum=cos(M_PI3.14159265358979323846(a_rho-a_omega))-cos(M_PI3.14159265358979323846a_rho)
1470	-cos(M_PI3.14159265358979323846*a_omega)+1.;
1471	// Regge propagator products
1472	double regge_omega_omega_P_cut=regge_omega*regge_omega_P_cut
1473	(cos(M_PI3.14159265358979323846(a_omega-0.5a_omega_P))-cos(M_PI_21.57079632679489661923a_omega_P)) ;
1474	double regge_omega_omega_f2_cut=regge_omega*regge_omega_f2_cut
1475	(cos(M_PI3.14159265358979323846(a_omega-0.5a_omega_f2))-cos(M_PI_21.57079632679489661923a_omega_f2));
1476	double regge_rho_omega_P_cut=regge_rho*regge_omega_P_cut
1477	(cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_P))-cos(M_PI_21.57079632679489661923a_omega_P));
1478	double regge_rho_omega_f2_cut=regge_rho*regge_omega_f2_cut
1479	(cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2))-cos(M_PI_21.57079632679489661923a_omega_f2));
1480	double regge_omega_f2_omega_P_cuts
1481	=regge_omega_f2_cutregge_omega_P_cutcos(M_PI_21.57079632679489661923*(a_omega_f2-a_omega_P));
1482	double regge_omega_P_cut_sq=regge_omega_P_cut*regge_omega_P_cut;
1483	double regge_omega_f2_cut_sq=regge_omega_f2_cut*regge_omega_f2_cut;
1484	double regge_rho_P_cut_sq=regge_rho_P_cut*regge_rho_P_cut;
1485	double regge_rho_f2_cut_sq=regge_rho_f2_cut*regge_rho_f2_cut;
1486	double regge_rho_rho_P_cut=regge_rho*regge_rho_P_cut
1487	(cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_P))-cos(M_PI_21.57079632679489661923a_rho_P));
1488	double regge_rho_rho_f2_cut=regge_rho*regge_rho_f2_cut
1489	(cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2))-cos(M_PI_21.57079632679489661923a_rho_f2));
1490	double regge_omega_rho_P_cut=regge_omega*regge_rho_P_cut
1491	(cos(M_PI3.14159265358979323846(a_omega-0.5a_rho_P))-cos(M_PI_21.57079632679489661923a_rho_P));
1492	double regge_omega_rho_f2_cut=regge_omega*regge_rho_f2_cut
1493	(cos(M_PI3.14159265358979323846(a_omega-0.5a_rho_f2))-cos(M_PI_21.57079632679489661923a_rho_f2));
1494	double regge_rho_P_omega_P_cuts=regge_rho_P_cut*regge_omega_P_cut
1495	cos(M_PI_21.57079632679489661923(a_rho_P-a_omega_P));
1496	double regge_rho_P_omega_f2_cuts=regge_rho_P_cut*regge_omega_f2_cut
1497	cos(M_PI_21.57079632679489661923(a_rho_P-a_omega_f2));
1498	double regge_rho_f2_omega_P_cuts=regge_rho_f2_cut*regge_omega_P_cut
1499	cos(M_PI_21.57079632679489661923(a_rho_f2-a_omega_P));
1500	double regge_rho_f2_omega_f2_cuts=regge_rho_f2_cut*regge_omega_f2_cut
1501	cos(M_PI_21.57079632679489661923(a_rho_f2-a_omega_f2));
1502	double regge_rho_f2_rho_P_cuts=regge_rho_f2_cut*regge_rho_P_cut
1503	cos(M_PI_21.57079632679489661923(a_rho_P-a_rho_f2));
1504
1505
1506	// Compute the square of the amplitude, including interference
1507	double T=0.;
1508	if (interfere==false){
1509	// Compute contributions to square of amplitude
1510	double Pi_R_sq=gRsq/(ReBReB+ImBImB);
1511	double aS_aS_rho_V_and_T_sq
1512	=Pi_R_sqgsq_rho_Sgsq_rho_V_and_T*regge_rho_sq;
1513	double aS_aS_omega_V_sq=Pi_R_sqgsq_omega_Vgsq_omega_S
1514	*(regge_omega_sq
1515	+ C_omega_P_cut*regge_omega_omega_P_cut
1516	+ C_omega_f2_cut*regge_omega_omega_f2_cut
1517	+ Csq_omega_P_cut*regge_omega_P_cut_sq
1518	+ Csq_omega_f2_cut*regge_omega_f2_cut_sq
1519	+ 2.C_omega_f2_C_omega_Pregge_omega_f2_omega_P_cuts);
1520
1521	double aS_aS_omega_V_rho_V_and_T
1522	=Pi_R_sqg_rho_Sg_omega_Sg_rho_V_and_Tg_omega_V
1523	(0.5regge_omegaregge_rhocos_rho_omega_sum
1524	+ C_omega_P_cut*regge_rho_omega_P_cut
1525	+ C_omega_f2_cut*regge_rho_omega_f2_cut
1526	);
1527	double aS_aS_rho_V_rho_V_and_T=Pi_R_sqgsq_rho_Sg_rho_V*g_rho_V_and_T
1528	(C_rho_P_cutregge_rho_rho_P_cut
1529	+ C_rho_f2_cut*regge_rho_rho_f2_cut);
1530	double aS_aS_rho_V_omega_V
1531	=Pi_R_sqg_rho_Sg_omega_Sg_rho_Vg_omega_V
1532	(C_rho_P_cutregge_omega_rho_P_cut
1533	+ C_rho_f2_cut*regge_omega_rho_f2_cut
1534	+ 2.C_rho_P_cutC_omega_P_cut*regge_rho_P_omega_P_cuts
1535	+ 2.C_rho_P_cutC_omega_f2_cut*regge_rho_P_omega_f2_cuts
1536	+ 2.C_rho_f2_cutC_omega_P_cut*regge_rho_f2_omega_P_cuts
1537	+ 2.C_rho_f2_cutC_omega_f2_cut*regge_rho_f2_omega_f2_cuts);
1538	double aS_aS_rho_cut_sq=Pi_R_sqgsq_rho_Sgsq_rho_V
1539	(Csq_rho_P_cutregge_rho_P_cut_sq
1540	+ Csq_rho_f2_cut*regge_rho_f2_cut_sq
1541	+ 2.C_rho_P_cutC_rho_f2_cut*regge_rho_f2_rho_P_cuts);
1542	double aS_aS=aS_aS_rho_V_and_T_sq + aS_aS_omega_V_sq + aS_aS_rho_V_rho_V_and_T
1543	+ aS_aS_omega_V_rho_V_and_T+aS_aS_rho_V_omega_V + aS_aS_rho_cut_sq;
1544
1545	double bS_bS=4Pi_R_sqgsq_rho_Sg_rho_Tg_rho_T*regge_rho_sq;
1546
1547	double aS_bS_rho_T_rho_V_and_T
1548	=-4.Pi_R_sqgsq_rho_Sg_rho_Tg_rho_V_and_T*regge_rho_sq;
1549	double aS_bS_rho_cuts=-4.Pi_R_sqgsq_rho_Sg_rho_Tg_rho_V
1550	(C_rho_P_cutregge_rho_rho_P_cut
1551	+ C_rho_f2_cut*regge_rho_rho_f2_cut);
1552	double aS_bS_omega_cuts
1553	=-4.Pi_R_sqg_rho_Sg_omega_Sg_rho_T*g_omega_V
1554	(C_omega_P_cutregge_rho_omega_P_cut
1555	+ C_omega_f2_cut*regge_rho_omega_f2_cut);
1556	double aS_bS_omega_V=-Pi_R_sqg_rho_Sg_omega_Sg_rho_Tg_omega_V
1557	regge_rhoregge_omega*cos_rho_omega_sum;
1558	double aS_bS=aS_bS_rho_T_rho_V_and_T+aS_bS_rho_cuts+aS_bS_omega_cuts+aS_bS_omega_V;
1559
1560	// Cut contribution from b1 exchange. We ignore the pole term.
1561	double Tb1=-0.5tkin1*Pi_R_sq
1562	(gsq_omega_Sgsq_omega_V(Csq_omega_P_cutregge_omega_P_cut_sq
1563	+ Csq_omega_f2_cut*regge_omega_f2_cut_sq
1564	+ 2.C_omega_P_cutC_omega_f2_cut
1565	*regge_omega_f2_omega_P_cuts)
1566	+gsq_rho_Sgsq_rho_V(Csq_rho_P_cut*regge_rho_P_cut_sq
1567	+ Csq_rho_f2_cut*regge_rho_f2_cut_sq
1568	+ 2.C_rho_P_cutC_rho_f2_cut
1569	*regge_rho_f2_rho_P_cuts)
1570	+g_rho_Vg_omega_Vg_rho_S*g_omega_S
1571	(2.C_rho_P_cutC_omega_P_cutregge_rho_P_omega_P_cuts
1572	+ 2.C_rho_f2_cutC_omega_f2_cut*regge_rho_f2_omega_f2_cuts
1573	+ 2.C_rho_P_cutC_omega_f2_cut*regge_rho_P_omega_f2_cuts
1574	+ 2.C_rho_f2_cutC_omega_P_cut*regge_rho_f2_omega_P_cuts));
1575
1576	T=aS_aSkin_aS_aS+aS_bSkin_aS_bS+bS_bS*kin_bS_bS + Tb1;
1577	}
1578	else{
1579	double bw1=gR/(ReBReB+ImBImB);
1580	double bw2=gR2/(ReB2ReB2+ImB2ImB2);
1581	double ReBWfac=ReBReB2+ImBImB2;
1582	double ImBWfac=ImBReB2-ReBImB2;
1583	double g1rho=sqrt(gsq_rho_S);
1584	double g1omega=sqrt(gsq_omega_S);
1585	double g2rho=sqrt(gsq_rho_S2);
1586	double g2omega=sqrt(gsq_omega_S2);
1587
1588	//interference factors
1589	double cos_rho_rho_P_sum=cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_P))-cos(M_PI_21.57079632679489661923*a_rho_P);
1590	double sin_rho_rho_P_sum=sin(M_PI3.14159265358979323846(a_rho-0.5a_rho_P))+sin(M_PI_21.57079632679489661923*a_rho_P);
1591	double cos_rho_rho_f2_sum=cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2))-cos(M_PI_21.57079632679489661923*a_rho_f2);
1592	double sin_rho_rho_f2_sum=sin(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2))+sin(M_PI_21.57079632679489661923*a_rho_f2);
1593	double cos_rho_omega_P_sum=cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_P))
1594	-cos(M_PI_21.57079632679489661923*a_omega_P);
1595	double sin_rho_omega_P_sum=sin(M_PI3.14159265358979323846(a_rho-0.5a_omega_P))
1596	+sin(M_PI_21.57079632679489661923*a_omega_P);
1597	double cos_rho_omega_f2_sum=cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2))
1598	-cos(M_PI_21.57079632679489661923*a_omega_f2);
1599	double sin_rho_omega_f2_sum=sin(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2))
1600	+sin(M_PI_21.57079632679489661923*a_omega_f2);
1601	double sin_rho_omega_sum=sin(M_PI3.14159265358979323846a_rho)-sin(M_PI3.14159265358979323846a_omega)
1602	-sin(M_PI3.14159265358979323846*(a_rho-a_omega));
1603	double cos_omega_rho_P_sum=cos(M_PI3.14159265358979323846(a_omega-0.5a_rho_P))-cos(M_PI_21.57079632679489661923*a_rho_P);
1604	double sin_omega_rho_P_sum=sin(M_PI3.14159265358979323846(a_omega-0.5a_rho_P))+sin(M_PI_21.57079632679489661923*a_rho_P);
1605	double cos_omega_rho_f2_sum=cos(M_PI3.14159265358979323846(a_omega-0.5a_rho_f2))-cos(M_PI_21.57079632679489661923*a_rho_f2);
1606	double sin_omega_rho_f2_sum=sin(M_PI3.14159265358979323846(a_omega-0.5a_rho_f2))+sin(M_PI_21.57079632679489661923*a_rho_f2);
1607	double cos_omega_P_rho_P=cos(M_PI_21.57079632679489661923*(a_omega_P-a_rho_P));
1608	double sin_omega_P_rho_P=sin(M_PI_21.57079632679489661923*(a_omega_P-a_rho_P));
1609	double cos_omega_P_rho_f2=cos(M_PI_21.57079632679489661923*(a_omega_P-a_rho_f2));
1610	double sin_omega_P_rho_f2=sin(M_PI_21.57079632679489661923*(a_omega_P-a_rho_f2));
1611	double cos_omega_f2_rho_P=cos(M_PI_21.57079632679489661923*(a_omega_f2-a_rho_P));
1612	double sin_omega_f2_rho_P=sin(M_PI_21.57079632679489661923*(a_omega_f2-a_rho_P));
1613	double cos_omega_f2_rho_f2=cos(M_PI_21.57079632679489661923*(a_omega_f2-a_rho_f2));
1614	double sin_omega_f2_rho_f2=sin(M_PI_21.57079632679489661923*(a_omega_f2-a_rho_f2));
1615
1616
1617	double omega_cut_interference
1618	=Csq_omega_P_cutregge_omega_P_cutregge_omega_P_cut
1619	+2.C_omega_P_cutC_omega_f2_cut*regge_omega_f2_omega_P_cuts
1620	+Csq_omega_f2_cutregge_omega_f2_cutregge_omega_f2_cut;
1621	double rho_cut_interference
1622	=Csq_rho_P_cutregge_rho_P_cutregge_rho_P_cut
1623	+2.C_rho_P_cutC_rho_f2_cut*regge_rho_f2_rho_P_cuts
1624	+Csq_rho_f2_cutregge_rho_f2_cutregge_rho_f2_cut;
1625	double Re_rho_cut_omega_cut_interference_12
1626	= 2.C_rho_P_cutC_omega_P_cutregge_rho_P_cutregge_omega_P_cut
1627	(ReBWfaccos_omega_P_rho_P-ImBWfac*sin_omega_P_rho_P)
1628	+2.C_rho_f2_cutC_omega_f2_cutregge_rho_f2_cutregge_omega_f2_cut
1629	(ReBWfaccos_omega_f2_rho_f2-ImBWfac*sin_omega_f2_rho_f2)
1630	+2.C_rho_f2_cutC_omega_P_cutregge_rho_f2_cutregge_omega_P_cut
1631	(ReBWfaccos_omega_P_rho_f2-ImBWfac*sin_omega_P_rho_f2)
1632	+2.C_rho_P_cutC_omega_f2_cutregge_rho_P_cutregge_omega_f2_cut
1633	(ReBWfaccos_omega_f2_rho_P-ImBWfac*sin_omega_f2_rho_P);
1634	double Re_rho_cut_omega_cut_interference_21
1635	=2.C_rho_P_cutC_omega_P_cutregge_rho_P_cutregge_omega_P_cut
1636	(ReBWfaccos_omega_P_rho_P+ImBWfac*sin_omega_P_rho_P)
1637	+2.C_rho_f2_cutC_omega_f2_cutregge_rho_f2_cutregge_omega_f2_cut
1638	(ReBWfaccos_omega_f2_rho_f2+ImBWfac*sin_omega_f2_rho_f2)
1639	+2.C_rho_f2_cutC_omega_P_cutregge_rho_f2_cutregge_omega_P_cut
1640	(ReBWfaccos_omega_P_rho_f2+ImBWfac*sin_omega_P_rho_f2)
1641	+2.C_rho_P_cutC_omega_f2_cutregge_rho_P_cutregge_omega_f2_cut
1642	(ReBWfaccos_omega_f2_rho_P+ImBWfac*sin_omega_f2_rho_P);
1643	double Im_rho_cut_omega_cut_interference_12
1644	=2.C_rho_P_cutC_omega_P_cutregge_rho_P_cutregge_omega_P_cut
1645	(ImBWfaccos_omega_P_rho_P+ReBWfac*sin_omega_P_rho_P)
1646	+2.C_rho_f2_cutC_omega_f2_cutregge_rho_f2_cutregge_omega_f2_cut
1647	(ImBWfaccos_omega_f2_rho_f2+ReBWfac*sin_omega_f2_rho_f2)
1648	+2.C_rho_f2_cutC_omega_P_cutregge_rho_f2_cutregge_omega_P_cut
1649	(ImBWfaccos_omega_P_rho_f2+ReBWfac*sin_omega_P_rho_f2)
1650	+2.C_rho_P_cutC_omega_f2_cutregge_rho_P_cutregge_omega_f2_cut
1651	(ImBWfaccos_omega_f2_rho_P+ReBWfac*sin_omega_f2_rho_P);
1652	double Im_rho_cut_omega_cut_interference_21
1653	=2.C_rho_P_cutC_omega_P_cutregge_rho_P_cutregge_omega_P_cut
1654	(ImBWfaccos_omega_P_rho_P-ReBWfac*sin_omega_P_rho_P)
1655	+2.C_rho_f2_cutC_omega_f2_cutregge_rho_f2_cutregge_omega_f2_cut
1656	(ImBWfaccos_omega_f2_rho_f2-ReBWfac*sin_omega_f2_rho_f2)
1657	+2.C_rho_f2_cutC_omega_P_cutregge_rho_f2_cutregge_omega_P_cut
1658	(ImBWfaccos_omega_P_rho_f2-ReBWfac*sin_omega_P_rho_f2)
1659	+2.C_rho_P_cutC_omega_f2_cutregge_rho_P_cutregge_omega_f2_cut
1660	(ImBWfaccos_omega_f2_rho_P-ReBWfac*sin_omega_f2_rho_P);
1661
1662	double aS_aS_fac1
1663	=g1rhog2rhogsq_rho_V_and_T2.regge_rho_sq
1664	+g1omegag2omegagsq_omega_V
1665	(2.regge_omega_sq + 2.*omega_cut_interference
1666	+2.C_omega_f2_cutregge_omega_omega_f2_cut
1667	+2.C_omega_P_cutregge_omega_omega_P_cut
1668	)
1669	+2.g1rhog2rhogsq_rho_Vrho_cut_interference;
1670	double aS_bS_fac1
1671	=g1rhog2rhog_rho_V_and_T4.regge_rho_sq
1672	+2.g1rhog2rhog_rho_V(C_rho_P_cut*regge_rho_rho_P_cut
1673	+C_rho_f2_cut*regge_rho_rho_f2_cut);
1674
1675	double ReTint
1676	=kin_aS_aS(aS_aS_fac1ReBWfac
1677	+2.g1rhog2rhog_rho_Vg_rho_V_and_T*regge_rho
1678	(C_rho_P_cutregge_rho_P_cut
1679	(ReBWfaccos_rho_rho_P_sum-ImBWfac*sin_rho_rho_P_sum)
1680	+C_rho_f2_cut*regge_rho_f2_cut
1681	(ReBWfaccos_rho_rho_f2_sum-ImBWfac*sin_rho_rho_f2_sum)
1682	)
1683	+g1omegag2rhog_omega_V*g_rho_V_and_T
1684	(0.5regge_rho*regge_omega
1685	(ReBWfaccos_rho_omega_sum+ImBWfac*sin_rho_omega_sum)
1686	+C_rho_P_cutregge_rhoregge_rho_P_cut
1687	(ReBWfaccos_rho_omega_P_sum-ImBWfac*sin_rho_omega_P_sum)
1688	+C_rho_f2_cutregge_rhoregge_rho_f2_cut
1689	(ReBWfaccos_rho_omega_f2_sum-ImBWfac*sin_rho_omega_f2_sum)
1690	)
1691	+g2omegag1rhog_omega_V*g_rho_V_and_T
1692	(0.5regge_rho*regge_omega
1693	(ReBWfaccos_rho_omega_sum-ImBWfac*sin_rho_omega_sum)
1694	+C_rho_P_cutregge_rhoregge_rho_P_cut
1695	(ReBWfaccos_rho_omega_P_sum+ImBWfac*sin_rho_omega_P_sum)
1696	+C_rho_f2_cutregge_rhoregge_rho_f2_cut
1697	(ReBWfaccos_rho_omega_f2_sum+ImBWfac*sin_rho_omega_f2_sum)
1698	)
1699	+g1rhog2omegag_rho_V*g_omega_V
1700	(C_rho_P_cutregge_omega*regge_rho_P_cut
1701	(ReBWfaccos_omega_rho_P_sum-ImBWfac*sin_omega_rho_P_sum)
1702	+C_rho_f2_cutregge_omegaregge_rho_P_cut
1703	(ReBWfaccos_omega_rho_f2_sum-ImBWfac*sin_omega_rho_f2_sum)
1704	+Re_rho_cut_omega_cut_interference_12
1705	)
1706	+g2rhog1omegag_rho_V*g_omega_V
1707	(C_rho_P_cutregge_omega*regge_rho_P_cut
1708	(ReBWfaccos_omega_rho_P_sum+ImBWfac*sin_omega_rho_P_sum)
1709	+C_rho_f2_cutregge_omegaregge_rho_P_cut
1710	(ReBWfaccos_omega_rho_f2_sum+ImBWfac*sin_omega_rho_f2_sum)
1711	+Re_rho_cut_omega_cut_interference_21
1712	)
1713	)
1714	+kin_aS_bS(-2.g_rho_T)
1715	(aS_bS_fac1ReBWfac
1716	+g_omega_Vg1rhog2omega
1717	(ReBWfac(0.5regge_rhoregge_omega*cos_rho_omega_sum
1718	+C_omega_P_cut*regge_rho_omega_P_cut
1719	+C_omega_f2_cut*regge_rho_omega_f2_cut)
1720	-ImBWfac(0.5regge_rhoregge_omegasin_rho_omega_sum
1721	-C_omega_P_cutregge_rhoregge_omega_P_cut
1722	*sin_rho_omega_P_sum
1723	-C_omega_f2_cutregge_rhoregge_omega_f2_cut
1724	*sin_rho_omega_f2_sum)
1725	)
1726	+g_omega_Vg2rhog1omega
1727	(ReBWfac(0.5regge_rhoregge_omega*cos_rho_omega_sum
1728	+C_omega_P_cut*regge_rho_omega_P_cut
1729	+C_omega_f2_cut*regge_rho_omega_f2_cut)
1730	+ImBWfac(0.5regge_rhoregge_omegasin_rho_omega_sum
1731	-C_omega_P_cutregge_rhoregge_omega_P_cut
1732	*sin_rho_omega_P_sum
1733	-C_omega_f2_cutregge_rhoregge_omega_f2_cut
1734	*sin_rho_omega_f2_sum)
1735	)
1736	)
1737	-0.5tkin1 // b1 contribution
1738	(2.rho_cut_interferenceg1rhog2rhog_rho_Vg_rho_V*ReBWfac
1739	+2.omega_cut_interferenceg1omegag2omegag_omega_Vg_omega_VReBWfac
1740	+g2rhog1omegag_rho_Vg_omega_VRe_rho_cut_omega_cut_interference_21
1741	+g1rhog2omegag_rho_Vg_omega_VRe_rho_cut_omega_cut_interference_12
1742	)
1743	+kin_bS_bS8.g1rhog2rhog_rho_Tg_rho_Tregge_rho_sq*ReBWfac;
1744
1745	double ImTint
1746	=kin_aS_aS(aS_aS_fac1ImBWfac
1747	+2.g1rhog2rhog_rho_Vg_rho_V_and_T*regge_rho
1748	(C_rho_P_cutregge_rho_P_cut
1749	(ImBWfaccos_rho_rho_P_sum+ReBWfac*sin_rho_rho_P_sum)
1750	+C_rho_f2_cut*regge_rho_f2_cut
1751	(ImBWfaccos_rho_rho_f2_sum+ReBWfac*sin_rho_rho_f2_sum)
1752	)
1753	+g1omegag2rhog_omega_V*g_rho_V_and_T
1754	(0.5regge_rho*regge_omega
1755	(ImBWfaccos_rho_omega_sum-ReBWfac*sin_rho_omega_sum)
1756	+C_rho_P_cutregge_rhoregge_rho_P_cut
1757	(ImBWfaccos_rho_omega_P_sum+ReBWfac*sin_rho_omega_P_sum)
1758	+C_rho_f2_cutregge_rhoregge_rho_f2_cut
1759	(ImBWfaccos_rho_omega_f2_sum+ReBWfac*sin_rho_omega_f2_sum)
1760	)
1761	+g2omegag1rhog_omega_V*g_rho_V_and_T
1762	(0.5regge_rho*regge_omega
1763	(ImBWfaccos_rho_omega_sum+ReBWfac*sin_rho_omega_sum)
1764	+C_rho_P_cutregge_rhoregge_rho_P_cut
1765	(ImBWfaccos_rho_omega_P_sum-ReBWfac*sin_rho_omega_P_sum)
1766	+C_rho_f2_cutregge_rhoregge_rho_f2_cut
1767	(ImBWfaccos_rho_omega_f2_sum-ReBWfac*sin_rho_omega_f2_sum)
1768	)
1769	+g1rhog2omegag_rho_V*g_omega_V
1770	(C_rho_P_cutregge_omega*regge_rho_P_cut
1771	(ImBWfaccos_omega_rho_P_sum+ReBWfac*sin_omega_rho_P_sum)
1772	+C_rho_f2_cutregge_omegaregge_rho_P_cut
1773	(ImBWfaccos_omega_rho_f2_sum+ReBWfac*sin_omega_rho_f2_sum)
1774	+Im_rho_cut_omega_cut_interference_12
1775	)
1776	+g2rhog1omegag_rho_V*g_omega_V
1777	(C_rho_P_cutregge_omega*regge_rho_P_cut
1778	(ImBWfaccos_omega_rho_P_sum-ReBWfac*sin_omega_rho_P_sum)
1779	+C_rho_f2_cutregge_omegaregge_rho_P_cut
1780	(ImBWfaccos_omega_rho_f2_sum-ReBWfac*sin_omega_rho_f2_sum)
1781	+Im_rho_cut_omega_cut_interference_21
1782	)
1783	)
1784	+kin_aS_bS(-2.g_rho_T)
1785	(aS_bS_fac1ImBWfac
1786	+g_omega_Vg1rhog2omega
1787	(ImBWfac(0.5regge_rhoregge_omega*cos_rho_omega_sum
1788	+C_omega_P_cut*regge_rho_omega_P_cut
1789	+C_omega_f2_cut*regge_rho_omega_f2_cut)
1790	+ReBWfac(0.5regge_rhoregge_omegasin_rho_omega_sum
1791	-C_omega_P_cutregge_rhoregge_omega_P_cut
1792	*sin_rho_omega_P_sum
1793	-C_omega_f2_cutregge_rhoregge_omega_f2_cut
1794	*sin_rho_omega_f2_sum)
1795	)
1796	+g_omega_Vg2rhog1omega
1797	(ImBWfac(0.5regge_rhoregge_omega*cos_rho_omega_sum
1798	+C_omega_P_cut*regge_rho_omega_P_cut
1799	+C_omega_f2_cut*regge_rho_omega_f2_cut)
1800	-ReBWfac(0.5regge_rhoregge_omegasin_rho_omega_sum
1801	-C_omega_P_cutregge_rhoregge_omega_P_cut
1802	*sin_rho_omega_P_sum
1803	-C_omega_f2_cutregge_rhoregge_omega_f2_cut
1804	*sin_rho_omega_f2_sum)
1805	)
1806	)
1807	-0.5tkin1 // b1 contribution
1808	(2.rho_cut_interferenceg1rhog2rhog_rho_Vg_rho_V*ImBWfac
1809	+2.omega_cut_interferenceg1omegag2omegag_omega_Vg_omega_VImBWfac
1810	+g2rhog1omegag_rho_Vg_omega_VIm_rho_cut_omega_cut_interference_21
1811	+g1rhog2omegag_rho_Vg_omega_VIm_rho_cut_omega_cut_interference_12
1812	)
1813	+kin_bS_bS8.g1rhog2rhog_rho_Tg_rho_Tregge_rho_sq*ImBWfac;
1814
1815	T=bw1bw2(ReTintcos(phase)+ImTintsin(phase));
1816	}
1817	// printf("kin %f %f %f \n",kin_aS_aS,kin_aS_bS,kin_bS_bS);
1818
1819	// Compute cross section
1820	double m1_plus_m2=m1+m2;
1821	double m1_minus_m2=m1-m2;
1822	double k=sqrt((ms_sq-m1_plus_m2m1_plus_m2)(ms_sq-m1_minus_m2*m1_minus_m2))
1823	/(2.*sqrt(ms_sq));
1824	double hbarc_sq=389.; // Convert to micro-barns
1825	double xsec=-hbarc_sqkT/(256.M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846*mp_sq_minus_s_sq);
1826
1827	return(xsec);
1828
1829
1830	}
1831
1832	double TensorCrossSection(TLorentzVector &q /* beam */,
1833	vector<Particle_t>&particle_types,
1834	vector<TLorentzVector>&particles,
1835	double gR,double ReB, double ImB){
1836	int two_particles=particle_types[0]+particle_types[1];
1837
1838	// Four vectors
1839	TLorentzVector p1(0,0,0.,ParticleMass(Proton));
1840	TLorentzVector p2=particles[2];
1841	TLorentzVector dp=p2-p1;
1842
1843	// Mandelstam variables
1844	double s=(q+p1).M2();
1845	double t=(p1-p2).M2();
1846
1847	// dot products
1848	double p1_dot_dp=p1.Dot(dp);
1849	double p2_dot_dp=p2.Dot(dp);
1850	double p1_dot_p2=p1.Dot(p2);
1851
1852	// momentum transfer compenents
1853	double dpx=dp.X();
1854	double dpy=dp.Y();
1855	double dpx2_plus_dpy2=dpxdpx+dpydpy;
1856
1857	// other constants
1858	double m_rho=0.775; // GeV
1859	double m_rho_sq=m_rho*m_rho;
1860
1861	// Coupling constants
1862	double f=100.; // scale factor to account for normalization of regge factor
1863	// to data??
1864	double gT_sq=f(2./3.)150.; // GeV^2
1865	if (two_particles==(7+17)){
1866	f=1.; // guess
1867	gT_sq=f*183.675; // GeV^2
1868	}
1869
1870	// s scale for regge trajectories
1871	double s0=1.;
1872
1873	// Regge trajectory for rho
1874	double a_rho=0.55+0.8*t;
1875	double a_rho_prime=0.8;
1876	double regge_rho=pow(s/s0,a_rho-1.)M_PI3.14159265358979323846a_rho_prime/(sin(M_PI3.14159265358979323846a_rho)TMath::Gamma(a_rho)); // excluding phase factor
1877	double regge_rho_sq=regge_rho*regge_rho;
1878
1879	// Regge trajectory for odderon
1880	double a_odderon=1.+0.25*t;
1881	double a_odderon_prime=0.25;
1882	double regge_odderon=pow(s/s0,a_odderon-1.)M_PI3.14159265358979323846a_odderon_prime
1883	/(sin(M_PI3.14159265358979323846a_odderon)TMath::Gamma(a_odderon)); // excluding phase factor
1884	double regge_odderon_sq=regge_odderon*regge_odderon;
1885
1886	// coupling constant for tensor interaction at rhoNN vertex
1887	double Kappa_rho=6.1;
1888
1889	// Amplitude
1890	double one_minus_dpx_over_m_rho=1.-dpx/m_rho;
1891	double one_minus_dpy_over_m_rho=1.-dpy/m_rho;
1892	double common_fac=(38./9.)(M_PI3.14159265358979323846/2.)(1./137.);
1893	double vector_coupling
1894	=2.dpx2_plus_dpy2/m_rho_sqp1_dot_dp*(p2_dot_dp/m_rho_sq-1.)
1895	+(m_p_sq-p1_dot_p2)(one_minus_dpx_over_m_rhoone_minus_dpx_over_m_rho
1896	+one_minus_dpy_over_m_rho*one_minus_dpy_over_m_rho);
1897	double tensor_coupling
1898	=(2.t-dpx2_plus_dpy2)(-Kappa_rho
1899	+0.25Kappa_rhoKappa_rho*(1.+p1_dot_p2/m_p_sq));
1900	double amp_sum=gT_sq(vector_coupling+tensor_coupling)regge_rho_sq;
1901	double gT_odderon=0.;//5.*f; // need better guess
1902	if (two_particles==(7+17)){
1903	//gT_odderon=5.; // need better guess
1904	}
1905	amp_sum+=gT_odderongT_odderonvector_coupling*regge_odderon_sq;
1906	amp_sum-=2.cos(M_PI3.14159265358979323846(a_odderon-a_rho))sqrt(gT_sq)gT_odderon
1907	regge_rhoregge_odderon(vector_coupling-(2.t-dpx2_plus_dpy2)*Kappa_rho);
1908	double decay_weight=1.; // take care of angular distribution of decay
1909	double T=common_facamp_sumdecay_weight;
1910
1911	// Compute cross section
1912	double m1=ParticleMass(particle_types[0]);
1913	double m2=ParticleMass(particle_types[1]);
1914	double m1_plus_m2=m1+m2;
1915	double m1_minus_m2=m1-m2;
1916	double m_sq=(particles[0]+particles[1]).M2();
1917	double k=sqrt((m_sq-m1_plus_m2m1_plus_m2)(m_sq-m1_minus_m2*m1_minus_m2))
1918	/(2.*sqrt(m_sq));
1919	double mp_sq_minus_s=m_p_sq-s;
1920	double mp_sq_minus_s_sq=mp_sq_minus_s*mp_sq_minus_s;
1921	double hbarc_sq=389.; // Convert to micro-barns
1922	double xsec=-hbarc_sq(gRgR/(ReBReB+ImBImB))Tk/(256.M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846*mp_sq_minus_s_sq);
1923
1924	return(xsec);
1925
1926	}
1927
1928	double TensorBackgroundInterference(TLorentzVector &q /* beam */,
1929	vector<Particle_t>&particle_types,
1930	vector<TLorentzVector>&particles,
1931	double gR,double ReB, double ImB,
1932	double phase){
1933	int two_particles=particle_types[0]+particle_types[1];
1934
1935	// Four vectors
1936	TLorentzVector p1(0,0,0.,ParticleMass(Proton));
1937	TLorentzVector p2=particles[2];
1938	TLorentzVector dp=p2-p1;
1939	TLorentzVector v1=particles[0]-q;
1940	TLorentzVector v2=particles[1]-q;
1941	double q_dot_dp=q.Dot(dp);
1942	double q_dot_v1=q.Dot(v1);
1943	double dp_dot_v1=dp.Dot(v1);
1944	double q_dot_v2=q.Dot(v2);
1945	double dp_dot_v2=dp.Dot(v2);
1946	double v1sq=v1.M2();
1947	double v2sq=v2.M2();
1948	double b1=-q_dot_dpv1sq+q_dot_v1dp_dot_v1;
1949	double b2=-q_dot_dpv2sq+q_dot_v2dp_dot_v2;
1950	TLorentzVector c1=-dp_dot_v1q+q_dot_dpv1;
1951	TLorentzVector c2=-dp_dot_v2q-q_dot_dpv2;
1952	TLorentzVector d1=q_dot_v1v1-v1sqq;
1953	TLorentzVector d2=q_dot_v2v2-v2sqq;
1954	double d1x_plus_d1y=d1.Vect().x()+d1.Vect().y();
1955	double d2x_plus_d2y=d2.Vect().x()+d2.Vect().y();
1956	double c1x_plus_c1y=c1.Vect().x()+c1.Vect().y();
1957	double c2x_plus_c2y=c2.Vect().x()+c2.Vect().y();
1958
1959	// Mandelstam variables
1960	double s=(q+p1).M2();
1961	double t=(p1-p2).M2();
1962
1963	// dot products
1964	double p1_dot_dp=p1.Dot(dp);
1965	double p2_dot_dp=p2.Dot(dp);
1966	double p1_dot_p2=p1.Dot(p2);
1967	double c1_dot_p1=c1.Dot(p1);
1968	double c2_dot_p1=c2.Dot(p1);
1969	double c1_dot_p2=c1.Dot(p2);
1970	double c2_dot_p2=c2.Dot(p2);
1971	double d1_dot_p1=d1.Dot(p1);
1972	double d2_dot_p1=d2.Dot(p1);
1973	double d1_dot_dp=d1.Dot(dp);
1974	double d2_dot_dp=d2.Dot(dp);
1975	double c1_dot_dp=c1.Dot(dp);
1976	double c2_dot_dp=c2.Dot(dp);
1977
1978	// momentum transfer compenents
1979	double dpx=dp.X();
1980	double dpy=dp.Y();
1981	double dpx_plus_dpy=dpx+dpy;
1982	// other momentum components
1983	double v1x_plus_v1y=v1.Vect().x()+v1.Vect().y();
1984	double v2x_plus_v2y=v2.Vect().x()+v2.Vect().y();
1985
1986	// Coupling constants
1987	double f=100.; // scale factor to account for normalization of regge factor
1988	// to data??
1989	double gT_sq=f(2./3.)150.; // GeV^2
1990	if (two_particles==(7+17)){
1991	f=1.; // guess
1992	gT_sq=f*183.675; // GeV^2
1993	}
1994	double g_omega_V=15.;
1995	double g_rho_V=3.4;
1996	double g_rho_T=11.0; // GeV^-1
1997	double g_rho_V_and_T=g_rho_V+2.m_pg_rho_T;
1998
1999	// Rho propagator for top exchange for double-exchange diagrams
2000	double m_rho=0.7685;
2001	double m_rho_sq=m_rho*m_rho;
2002	double Gamma_rho=0.1462;
2003	double m_rhosq_minus_v1sq=m_rho*m_rho-v1sq;
2004	double Pi_rho_1_sq=1./(m_rhosq_minus_v1sq*m_rhosq_minus_v1sq
2005	+m_rhom_rhoGamma_rho*Gamma_rho);
2006	double m_rhosq_minus_v2sq=m_rho*m_rho-v2sq;
2007	double Pi_rho_2_sq=1./(m_rhosq_minus_v2sq*m_rhosq_minus_v2sq
2008	+m_rhom_rhoGamma_rho*Gamma_rho);
2009
2010	// s scale for regge trajectories
2011	double s0=1.;
2012
2013	// Regge trajectory for rho
2014	double a_rho=0.55+0.8*t;
2015	double a_rho_prime=0.8;
2016	double regge_rho=pow(s/s0,a_rho-1.)M_PI3.14159265358979323846a_rho_prime/(sin(M_PI3.14159265358979323846a_rho)TMath::Gamma(a_rho)); // excluding phase factor
2017	double regge_rho_sq=regge_rho*regge_rho;
2018
2019	// Regge trajectory for odderon
2020	/*
2021	double a_odderon=1.+0.25*t;
2022	double a_odderon_prime=0.25;
2023	double regge_odderon=pow(s/s0,a_odderon-1.)M_PIa_odderon_prime
2024	/(sin(M_PIa_odderon)TMath::Gamma(a_odderon)); // excluding phase factor
2025	double regge_odderon_sq=regge_odderon*regge_odderon;
2026	*/
2027
2028	// omega propagator for top exchange
2029	double m_omega=0.78265;
2030	double Gamma_omega=0.00849;
2031	double m_omegasq_minus_v1sq=m_omega*m_omega-v1sq;
2032	double Pi_omega_1_sq=1./(m_omegasq_minus_v1sq*m_omegasq_minus_v1sq
2033	+m_omegam_omegaGamma_omega*Gamma_omega);
2034	double m_omegasq_minus_v2sq=m_omega*m_omega-v2sq;
2035	double Pi_omega_2_sq=1./(m_omegasq_minus_v2sq*m_omegasq_minus_v2sq
2036	+m_omegam_omegaGamma_omega*Gamma_omega);
2037
2038	// Regge trajectory for omega
2039	double a_omega=0.44+0.9*t;
2040	double a_omega_prime=0.9;
2041	// double cos_omega=cos(M_PI*a_omega);
2042	double sin_omega=sin(M_PI3.14159265358979323846*a_omega);
2043	double regge_omega=pow(s/s0,a_omega-1.)M_PI3.14159265358979323846a_omega_prime/(sin_omega*TMath::Gamma(a_omega)); // excluding phase factor
2044	// interference between rho and omega;
2045	double cos_rho_omega=cos(M_PI3.14159265358979323846*(a_rho-a_omega));
2046	double sin_rho_omega=sin(M_PI3.14159265358979323846*(a_rho-a_omega));
2047
2048	// coupling constant for tensor interaction at rhoNN vertex
2049	double Kappa_rho=6.1;
2050
2051	/*
2052	double gT_odderon=0.;//5.*f; // need better guess
2053	if (two_particles==(7+17)){
2054	//gT_odderon=5.; // need better guess
2055	}
2056	*/
2057
2058	// terms for interference between resonance shape and propagator in top
2059	// part of background diagrams
2060	double Re_B_and_omega_1=m_omegasq_minus_v1sqReB+m_omegaGamma_omega*ImB;
2061	double Im_B_and_omega_1=m_omegaGamma_omegaReB-m_omegasq_minus_v1sq*ImB;
2062	double Re_B_and_omega_2=m_omegasq_minus_v2sqReB+m_omegaGamma_omega*ImB;
2063	double Im_B_and_omega_2=m_omegaGamma_omegaReB-m_omegasq_minus_v2sq*ImB;
2064	double Re_B_and_rho_1=m_rhosq_minus_v1sqReB+m_rhoGamma_rho*ImB;
2065	double Im_B_and_rho_1=m_rhoGamma_rhoReB-m_rhosq_minus_v1sq*ImB;
2066	double Re_B_and_rho_2=m_rhosq_minus_v2sqReB+m_rhoGamma_rho*ImB;
2067	double Im_B_and_rho_2=m_rhoGamma_rhoReB-m_rhosq_minus_v2sq*ImB;
2068
2069
2070	double ReT=0.,ImT=0.;
2071	// Kinematic factors
2072	double temp1_1=(v1x_plus_v1yc1_dot_p1-d1_dot_p1dpx_plus_dpy)*dpx_plus_dpy;
2073	double temp2_1=(v2x_plus_v2yc2_dot_p1-d2_dot_p1dpx_plus_dpy)*dpx_plus_dpy;
2074	double temp1_2=(v1x_plus_v1yc1_dot_p2+(b1-d1_dot_p1)dpx_plus_dpy)*dpx_plus_dpy;
2075	double temp2_2=(v2x_plus_v2yc2_dot_p2+(b2-d2_dot_p1)dpx_plus_dpy)*dpx_plus_dpy;
2076	double temp3_1=((b1-d1_dot_dp)dpx_plus_dpy+v1x_plus_v1yc1.Dot(dp))*dpx_plus_dpy;
2077	double temp3_2=((b2-d2_dot_dp)dpx_plus_dpy+v2x_plus_v2yc2.Dot(dp))*dpx_plus_dpy;
2078
2079	double kin1_1=temp1_1(p2_dot_dp/m_rho_sq-1.)+temp2_1p1_dot_dp/m_rho_sq
2080	+(m_p_sq-p1_dot_p2)(temp3_1/m_rho_sq-2.b1-v1x_plus_v1y*c1x_plus_c1y
2081	+dpx_plus_dpy*d1x_plus_d1y);
2082	double kin2_1=m_ptemp1_1((p1_dot_dp+p2_dot_dp)/m_rho_sq-1.);
2083	double kin3_1=temp1_1*p1_dot_dp;
2084	double kin4_1=temp3_1*(t/m_rho_sq-1.)
2085	+t((b1-d1_dot_dp)dpx_plus_dpydpx_plus_dpy/m_rho_sq-2.b1
2086	+dpx_plus_dpy*d1x_plus_d1y
2087	+v1x_plus_v1y(dpx_plus_dpyc1_dot_dp/m_rho_sq-c1x_plus_c1y));
2088	double kin1_2=temp1_2(p2_dot_dp/m_rho_sq-1.)+temp2_2p1_dot_dp/m_rho_sq
2089	+(m_p_sq-p1_dot_p2)(temp3_2/m_rho_sq-2.b2-v2x_plus_v2y*c2x_plus_c2y
2090	+dpx_plus_dpy*d2x_plus_d2y);
2091	double kin2_2=m_ptemp1_2((p1_dot_dp+p2_dot_dp)/m_rho_sq-1.);
2092	double kin3_2=temp1_2*p1_dot_dp/m_p;
2093	double kin4_2=temp3_2*(t/m_rho_sq-1.)
2094	+t((b2-d2_dot_dp)dpx_plus_dpydpx_plus_dpy/m_rho_sq-2.b2
2095	+dpx_plus_dpy*d2x_plus_d2y
2096	+v2x_plus_v2y(dpx_plus_dpyc2_dot_dp/m_rho_sq-c2x_plus_c2y));
2097
2098	// Compute imaginary and real contributions
2099	double zeta=1., C=1;
2100	if (two_particles==(7+7)){ // pi0 pi0
2101	zeta=1./sqrt(2.);
2102
2103	ReT=(kin1_1-0.25Kappa_rhokin4_1)
2104	(2.g_rho_V_and_Tregge_rho_sqPi_omega_1_sq*Re_B_and_omega_1
2105	+(2./3.)g_omega_Vregge_rhoregge_omegaPi_rho_1_sq
2106	(Re_B_and_rho_1cos_rho_omega-Im_B_and_rho_1*sin_rho_omega)
2107	)
2108	+(kin1_2-0.25Kappa_rhokin4_2)
2109	(2.g_rho_V_and_Tregge_rho_sqPi_omega_2_sq*Re_B_and_omega_2
2110	+(2./3.)g_omega_Vregge_rhoregge_omegaPi_rho_2_sq
2111	(Re_B_and_rho_2cos_rho_omega-Im_B_and_rho_2*sin_rho_omega)
2112	)
2113	-4.(kin2_1-Kappa_rhokin3_1)regge_rho_sqg_rho_T*Pi_omega_1_sq
2114	*Re_B_and_omega_1
2115	-4.(kin2_2-Kappa_rhokin3_2)regge_rho_sqg_rho_T*Pi_omega_2_sq
2116	*Re_B_and_omega_2;
2117	ImT=(kin1_1-0.25Kappa_rhokin4_1)
2118	(2.g_rho_V_and_Tregge_rho_sqPi_omega_1_sq*Im_B_and_omega_1
2119	+(2./3.)g_omega_Vregge_rhoregge_omegaPi_rho_1_sq
2120	(Re_B_and_rho_1sin_rho_omega+Im_B_and_rho_1*cos_rho_omega)
2121	)
2122	+(kin1_2-0.25Kappa_rhokin4_2)
2123	(2.g_rho_V_and_Tregge_rho_sqPi_omega_2_sq*Im_B_and_omega_2
2124	+(2./3.)g_omega_Vregge_rhoregge_omegaPi_rho_2_sq
2125	(Re_B_and_rho_2sin_rho_omega+Im_B_and_rho_2*cos_rho_omega)
2126	)
2127	-4.(kin2_1-Kappa_rhokin3_1)regge_rho_sqg_rho_T*Pi_omega_1_sq
2128	*Im_B_and_omega_1
2129	-4.(kin2_2-Kappa_rhokin3_2)regge_rho_sqg_rho_T*Pi_omega_2_sq
2130	*Im_B_and_omega_2;
2131	}
2132	if (two_particles==(7+17)){ // pi0 eta
2133	C=sqrt(2./3.);
2134
2135	ReT=(kin1_1-0.25Kappa_rhokin4_1)
2136	((2./3.)g_rho_V_and_Tregge_rho_sqPi_rho_1_sq*Re_B_and_rho_1
2137	+2.g_omega_Vregge_rhoregge_omegaPi_omega_1_sq
2138	(Re_B_and_omega_1cos_rho_omega-Im_B_and_omega_1*sin_rho_omega)
2139	)
2140	+(kin1_2-0.25Kappa_rhokin4_2)
2141	((2./3.)g_rho_V_and_Tregge_rho_sqPi_omega_2_sq*Re_B_and_omega_2
2142	+2.g_omega_Vregge_rhoregge_omegaPi_rho_2_sq
2143	(Re_B_and_rho_2cos_rho_omega-Im_B_and_rho_2*sin_rho_omega)
2144	)
2145	-4./3.(kin2_1-Kappa_rhokin3_1)regge_rho_sqg_rho_T*Pi_rho_1_sq
2146	*Re_B_and_rho_1
2147	-4./3.(kin2_2-Kappa_rhokin3_2)regge_rho_sqg_rho_T*Pi_omega_2_sq
2148	*Re_B_and_omega_2;
2149	ImT=(kin1_1-0.25Kappa_rhokin4_1)
2150	(2./3.g_rho_V_and_Tregge_rho_sqPi_rho_1_sq*Im_B_and_rho_1
2151	+(2.)g_omega_Vregge_rhoregge_omegaPi_omega_1_sq
2152	(Re_B_and_omega_1sin_rho_omega+Im_B_and_omega_1*cos_rho_omega)
2153	)
2154	+(kin1_2-0.25Kappa_rhokin4_2)
2155	(2./3.g_rho_V_and_Tregge_rho_sqPi_omega_2_sq*Im_B_and_omega_2
2156	+(2.)g_omega_Vregge_rhoregge_omegaPi_rho_2_sq
2157	(Re_B_and_rho_2sin_rho_omega+Im_B_and_rho_2*cos_rho_omega)
2158	)
2159	-4./3.(kin2_1-Kappa_rhokin3_1)regge_rho_sqg_rho_T*Pi_rho_1_sq
2160	*Im_B_and_rho_1
2161	-4./3.(kin2_2-Kappa_rhokin3_2)regge_rho_sqg_rho_T*Pi_omega_2_sq
2162	*Im_B_and_omega_2;
2163	}
2164	double m1=ParticleMass(particle_types[0]);
2165	double m2=ParticleMass(particle_types[1]);
2166	double m1_plus_m2=m1+m2;
2167	double m1_minus_m2=m1-m2;
2168	double Msq=(particles[0]+particles[1]).M2();
2169	double k=sqrt((Msq-m1_plus_m2m1_plus_m2)(Msq-m1_minus_m2*m1_minus_m2))
2170	/(2.*sqrt(Msq));
2171	double common_fac=(5./3.)kzetaCsqrt(M_PI3.14159265358979323846/137.gT_sqg0_sq)gR/(ReBReB+ImB*ImB);
2172
2173
2174	double T=common_fac(ReTcos(phase)+ImT*sin(phase));
2175
2176	// Compute cross section
2177	double mp_sq_minus_s=m_p_sq-s;
2178	double mp_sq_minus_s_sq=mp_sq_minus_s*mp_sq_minus_s;
2179	double hbarc_sq=389.; // Convert to micro-barns
2180	double xsec=-hbarc_sqT/(256.M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846mp_sq_minus_s_sq);
2181
2182	return(xsec);
2183
2184	}
2185
2186	double TensorScalarInterference(TLorentzVector &q /* beam */,
2187	vector<Particle_t>&particle_types,
2188	vector<TLorentzVector>&particles,
2189	double gR,double ReB, double ImB,
2190	double gR_S,double ReB_S,double ImB_S,
2191	double g_omega_S,double g_rho_S,
2192	double phase){
2193	int two_particles=particle_types[0]+particle_types[1];
2194	double m_rho=0.7685;
2195	double m_rho_sq=m_rho*m_rho;
2196
2197	// Four vectors
2198	TLorentzVector p1(0,0,0.,ParticleMass(Proton));
2199	TLorentzVector p2=particles[2];
2200	TLorentzVector dp=p2-p1;
2201
2202	// Mandelstam variables
2203	double s=(q+p1).M2();
2204	double t=(p1-p2).M2();
2205	double q_dot_dp=q.Dot(dp);
2206	double q_dot_p1=q.Dot(p1);
2207
2208	// dot products
2209	double p1_dot_dp=p1.Dot(dp);
2210	double p2_dot_dp=p2.Dot(dp);
2211	double p1_dot_p2=p1.Dot(p2);
2212
2213
2214	// momentum transfer compenents
2215	double dpx=dp.X();
2216	double dpy=dp.Y();
2217	double dpx_plus_dpy=dpx+dpy;
2218
2219	// Coupling constants
2220	double f=100.; // scale factor to account for normalization of regge factor
2221	// to data??
2222	double gT_sq=f(2./3.)150.; // GeV^2
2223	if (two_particles==(7+17)){
2224	f=1.; // guess
2225	gT_sq=f*183.675; // GeV^2
2226	}
2227	double g_omega_V=15.;
2228	double g_rho_V=3.4;
2229	double g_rho_T=11.0; // GeV^-1
2230	double g_rho_V_and_T=g_rho_V+2.m_pg_rho_T;
2231
2232	// s scale for regge trajectories
2233	double s0=1.;
2234
2235	// Regge trajectory for rho
2236	double a_rho=0.55+0.8*t;
2237	double a_rho_prime=0.8;
2238	double cos_rho=cos(M_PI3.14159265358979323846*a_rho);
2239	double sin_rho=sin(M_PI3.14159265358979323846*a_rho);
2240	double regge_rho=pow(s/s0,a_rho-1.)M_PI3.14159265358979323846a_rho_prime/(sin(M_PI3.14159265358979323846a_rho)TMath::Gamma(a_rho)); // excluding phase factor
2241	double regge_rho_sq=regge_rho*regge_rho;
2242
2243	/*
2244	// Regge trajectory for odderon
2245	double a_odderon=1.+0.25*t;
2246	double a_odderon_prime=0.25;
2247	double regge_odderon=pow(s/s0,a_odderon-1.)M_PIa_odderon_prime
2248	/(sin(M_PIa_odderon)TMath::Gamma(a_odderon)); // excluding phase factor
2249	double regge_odderon_sq=regge_odderon*regge_odderon;
2250	*/
2251
2252	// Regge trajectory for omega
2253	double a_omega=0.44+0.9*t;
2254	double a_omega_prime=0.9;
2255	//double cos_omega=cos(M_PI*a_omega);
2256	double sin_omega=sin(M_PI3.14159265358979323846*a_omega);
2257	double regge_omega=pow(s/s0,a_omega-1.)M_PI3.14159265358979323846a_omega_prime/(sin_omega*TMath::Gamma(a_omega)); // excluding phase factor
2258	//double regge_omega_sq=regge_omega*regge_omega;
2259	// interference between rho and omega;
2260	double cos_rho_omega=cos(M_PI3.14159265358979323846*(a_rho-a_omega));
2261	double sin_rho_omega=sin(M_PI3.14159265358979323846*(a_rho-a_omega));
2262
2263	// Regge cuts for rho
2264	double dc=regge_cuts[0];
2265	double a_rho_P=0.64+0.16*t; // Pomeron
2266	double a_rho_f2=0.222+0.404*t;
2267	double regge_rho_P_cut=exp(dct)pow(s/s0,a_rho_P-1.);
2268	double regge_rho_f2_cut=exp(dct)pow(s/s0,a_rho_f2-1.);
2269	double C_rho_P_cut=regge_cuts[3];
2270	double C_rho_f2_cut=regge_cuts[4];
2271
2272	// Regge cuts for omega
2273	double a_omega_P=0.52+0.196*t; // Pomeron
2274	double a_omega_f2=0.112+0.428*t;
2275	double regge_omega_P_cut=exp(dct)pow(s/s0,a_omega_P-1.);
2276	double regge_omega_f2_cut=exp(dct)pow(s/s0,a_omega_f2-1.);
2277	double C_omega_P_cut=regge_cuts[1];
2278	double C_omega_f2_cut=regge_cuts[2];
2279
2280	// cut interference factors
2281	double cos_rho_rho_P=cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_P));
2282	double cos_rho_rho_f2=cos(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2));
2283	double sin_rho_rho_P=sin(M_PI3.14159265358979323846(a_rho-0.5a_rho_P));
2284	double sin_rho_rho_f2=sin(M_PI3.14159265358979323846(a_rho-0.5a_rho_f2));
2285	double cos_rho_omega_P=cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_P));
2286	double cos_rho_omega_f2=cos(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2));
2287	double sin_rho_omega_P=sin(M_PI3.14159265358979323846(a_rho-0.5a_omega_P));
2288	double sin_rho_omega_f2=sin(M_PI3.14159265358979323846(a_rho-0.5a_omega_f2));
2289
2290	// coupling constant for tensor interaction at rhoNN vertex
2291	double Kappa_rho=6.1;
2292
2293	/*
2294	double gT_odderon=0.;//5.*f; // need better guess
2295	if (two_particles==(7+17)){
2296	//gT_odderon=5.; // need better guess
2297	}
2298	*/
2299
2300	// terms for interference between resonances
2301	double Re_B_S_and_T=ReBReB_S+ImBImB_S;
2302	double Im_B_S_and_T=ReBImB_S-ReB_SImB;
2303
2304	double ReT=0.,ImT=0.;
2305	// Kinematic factors
2306	double kin_aS=(4./3.)(m_p_sq-p1_dot_p2+0.5Kappa_rhot)q_dot_dp
2307	+(5./3.)dpx_plus_dpydpx_plus_dpy
2308	q_dot_p1((p1_dot_dp+p2_dot_dp)/m_rho_sq-1.);
2309	double kin_bS=(5./3.)dpx_plus_dpydpx_plus_dpy*q_dot_p1
2310	(m_p((p1_dot_dp+p2_dot_dp)/m_rho_sq-1.)
2311	-Kappa_rho/(2.m_p)p1_dot_dp(2.p2_dot_dp/m_rho_sq-1.));
2312	double common_fac=gRgR_Ssqrt(gT_sq*M_PI3.14159265358979323846/137.)
2313	/(ReBReB+ImBImB)/(ReB_SReB_S+ImB_SImB_S);
2314
2315	ReT=kin_aS(g_rho_V_and_Tg_rho_S*regge_rho_sq
2316	(Re_B_S_and_T(1.-cos_rho)+Im_B_S_and_T*sin_rho)
2317	+2.g_rho_Vg_rho_S*regge_rho
2318	(C_rho_P_cutregge_rho_P_cut(Re_B_S_and_Tcos_rho_rho_P
2319	-Im_B_S_and_T*sin_rho_rho_P)
2320	+C_rho_f2_cutregge_rho_f2_cut(Re_B_S_and_T*cos_rho_rho_f2
2321	-Im_B_S_and_T*sin_rho_rho_f2)
2322	)
2323	+g_omega_Vg_omega_Sregge_rho
2324	(regge_omega(Re_B_S_and_T*(cos_rho_omega-cos_rho)
2325	-Im_B_S_and_T*(sin_rho_omega+sin_rho))
2326	+2.C_omega_P_cutregge_omega_P_cut
2327	(Re_B_S_and_Tcos_rho_omega_P-Im_B_S_and_T*sin_rho_omega_P)
2328	+2.C_omega_f2_cutregge_omega_f2_cut
2329	(Re_B_S_and_Tcos_rho_omega_f2-Im_B_S_and_T*sin_rho_omega_f2)
2330	)
2331	)
2332	+kin_bS(-2.g_rho_Sg_rho_Tregge_rho_sq)
2333	(Re_B_S_and_T(1.-cos_rho)+Im_B_S_and_T*sin_rho);
2334
2335	ImT=kin_aS(g_rho_V_and_Tg_rho_S*regge_rho_sq
2336	(Re_B_S_and_Tsin_rho-Im_B_S_and_T*(1.-cos_rho))
2337	+2.g_rho_Vg_rho_S*regge_rho
2338	(C_rho_P_cutregge_rho_P_cut(Re_B_S_and_Tsin_rho_rho_P
2339	+Im_B_S_and_T*cos_rho_rho_P)
2340	+C_rho_f2_cutregge_rho_f2_cut(Re_B_S_and_T*sin_rho_rho_f2
2341	+Im_B_S_and_T*cos_rho_rho_f2)
2342	)
2343	+g_omega_Vg_omega_Sregge_rho
2344	(regge_omega(Re_B_S_and_T*(sin_rho_omega+sin_rho)
2345	+Im_B_S_and_T*(cos_rho_omega-cos_rho))
2346	+2.C_omega_P_cutregge_omega_P_cut
2347	(Re_B_S_and_Tsin_rho_omega_P+Im_B_S_and_T*cos_rho_omega_P)
2348	+2.C_omega_f2_cutregge_omega_f2_cut
2349	(Re_B_S_and_Tsin_rho_omega_f2+Im_B_S_and_T*cos_rho_omega_f2)
2350	)
2351	)
2352	+kin_bS(-2.g_rho_Sg_rho_Tregge_rho_sq)
2353	(Re_B_S_and_Tsin_rho-Im_B_S_and_T*(1.-cos_rho));
2354
2355	// Sqaure of amplitudes
2356	double T=common_fac(ReTcos(phase)+ImT*sin(phase));
2357
2358	// Compute cross section
2359	double m1=ParticleMass(particle_types[0]);
2360	double m2=ParticleMass(particle_types[1]);
2361	double m1_plus_m2=m1+m2;
2362	double m1_minus_m2=m1-m2;
2363	double m_sq=(particles[0]+particles[1]).M2();
2364	double k=sqrt((m_sq-m1_plus_m2m1_plus_m2)(m_sq-m1_minus_m2*m1_minus_m2))
2365	/(2.*sqrt(m_sq));
2366	double mp_sq_minus_s=m_p_sq-s;
2367	double mp_sq_minus_s_sq=mp_sq_minus_s*mp_sq_minus_s;
2368	double hbarc_sq=389.; // Convert to micro-barns
2369	double xsec=-hbarc_sqTk/(256.M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846M_PI3.14159265358979323846*mp_sq_minus_s_sq);
2370
2371	return xsec;
2372	}
2373
2374
2375	// Put particle data into hddm format and output to file
2376	void WriteEvent(unsigned int eventNumber,TLorentzVector &beam, float vert[3],
2377	vector<Particle_t>&particle_types,
2378	vector<TLorentzVector>&particle_vectors, s_iostream_t *file){
2379	s_PhysicsEvents_t* pes;
2380	s_Reactions_t* rs;
2381	s_Target_t* ta;
2382	s_Beam_t* be;
2383	s_Vertices_t* vs;
2384	s_Origin_t* origin;
2385	s_Products_t* ps;
2386	s_HDDM_t *thisOutputEvent = make_s_HDDM();
2387	thisOutputEvent->physicsEvents = pes = make_s_PhysicsEvents(1);
2388	pes->mult = 1;
2389	pes->in[0].runNo = runNo;
2390	pes->in[0].eventNo = eventNumber;
2391	pes->in[0].reactions = rs = make_s_Reactions(1);
2392	rs->mult = 1;
2393	// Beam
2394	rs->in[0].beam = be = make_s_Beam();
2395	be->type = Gamma;
2396	be->properties = make_s_Properties();
2397	be->properties->charge = ParticleCharge(be->type);
2398	be->properties->mass = ParticleMass(be->type);
2399	be->momentum = make_s_Momentum();
2400	be->momentum->px = 0.;
2401	be->momentum->py = 0.;
2402	be->momentum->pz = beam.Pz();
2403	be->momentum->E = beam.E();
2404	// Target
2405	rs->in[0].target = ta = make_s_Target();
2406	ta->type = Proton;
2407	ta->properties = make_s_Properties();
2408	ta->properties->charge = ParticleCharge(ta->type);
2409	ta->properties->mass = ParticleMass(ta->type);
2410	ta->momentum = make_s_Momentum();
2411	ta->momentum->px = 0.;
2412	ta->momentum->py = 0.;
2413	ta->momentum->pz = 0.;
2414	ta->momentum->E = ParticleMass(ta->type);
2415	// Primary vertex
2416	rs->in[0].vertices = vs = make_s_Vertices(1);
2417	vs->mult = 1;
2418	vs->in[0].origin = origin = make_s_Origin();
2419	vs->in[0].products = ps = make_s_Products(particle_vectors.size());
2420	ps->mult = 0;
2421	origin->t = 0.0;
2422	origin->vx = vert[0];
2423	origin->vy = vert[1];
2424	origin->vz = vert[2];
2425	// Final state particles
2426	for (unsigned int i=0;i<particle_vectors.size();i++,ps->mult++){
2427	Particle_t my_particle=particle_types[i];
2428	ps->in[ps->mult].type = my_particle;
2429	ps->in[ps->mult].pdgtype = PDGtype(my_particle);
2430	ps->in[ps->mult].id = i+1; /* unique value for this particle within the event */
2431	ps->in[ps->mult].parentid = 0; /* All internally generated particles have no parent */
2432	ps->in[ps->mult].mech = 0; // ???
2433	ps->in[ps->mult].momentum = make_s_Momentum();
2434	ps->in[ps->mult].momentum->px = particle_vectors[i].Px();
2435	ps->in[ps->mult].momentum->py = particle_vectors[i].Py();
2436	ps->in[ps->mult].momentum->pz = particle_vectors[i].Pz();
2437	ps->in[ps->mult].momentum->E = particle_vectors[i].E();
2438	}
2439	flush_s_HDDM(thisOutputEvent,file);
2440	}
2441
2442	// Create some diagnostic histograms
2443	void CreateHistograms(string beamConfigFile){
2444
2445	thrown_t=new TH1D("thrown_t","Thrown t distribution",1000,0.,5);
2446	thrown_t->SetXTitle("-t [GeV^{2}]");
2447	thrown_dalitzZ=new TH1D("thrown_dalitzZ","thrown dalitz Z",110,-0.05,1.05);
2448	thrown_Egamma=new TH1D("thrown_Egamma","Thrown E_{#gamma} distribution",
2449	1000,0,12.);
2450	thrown_Egamma->SetTitle("E_{#gamma} [GeV]");
2451	thrown_mass=new TH1D("thrown_mass","Thrown mass distribution",
2452	1000,0,4.);
2453	thrown_mass->SetXTitle("mass [GeV]");
2454	thrown_dalitzXY=new TH2D("thrown_dalitzXY","Dalitz distribution Y vs X",100,-1.,1.,100,-1.,1);
2455	thrown_mass_vs_E=new TH2D("thrown_mass_vs_E","M(4#gamma) vs Ebeam",
2456	48,2.8,12.4,450,0,4.5);
2457	thrown_mass_vs_E->SetYTitle("M(4#gamma) [GeV]");
2458	thrown_mass_vs_E->SetXTitle("E(beam) [GeV]");
2459
2460	thrown_theta_vs_p=new TH2D("thrown_theta_vs_p","Proton #theta_{LAB} vs. p",
2461	200,0,2.,180,0.,90.);
2462	thrown_theta_vs_p->SetXTitle("p [GeV/c]");
2463	thrown_theta_vs_p->SetYTitle("#theta [degrees]");
2464
2465	BeamProperties beamProp(beamConfigFile);
2466	cobrems_vs_E = (TH1D*)beamProp.GetFlux();
2467	}
2468
2469
2470	// Create a graph of the cross section dsigma/dt as a function of -t
2471	void GraphCrossSection(double m1,double m2){
2472	// beam energy in lab
2473	double Egamma=cobrems_vs_E->GetBinLowEdge(1); // get from CobremsGenerator histogram;
2474	TLorentzVector beam(0,0,Egamma,Egamma);
2475	TLorentzVector target(0,0,0,m_p);
2476
2477	// CM energy
2478	double s=m_p(m_p+2.Egamma);
2479	double Ecm=sqrt(s);
2480
2481	// Coupling constants
2482	double gsq_rho_S_gamma=0.02537;
2483	double gsq_omega_S_gamma=0.2283;
2484
2485	// Parameters for integration over line shape
2486	double m1_plus_m2=m1+m2;
2487	double m1sq=m1*m1;
2488	double m2sq=m2*m2;
2489	double m1sq_plus_m2sq=m1sq+m2sq;
2490	double m_max=m_p(sqrt(1.+2.Egamma/m_p)-1.);
2491	double dmrange=m_max-m1_plus_m2;
2492	double dm=dmrange/1000.;
2493	double gR=0.,ImB=0,ReB=0; // resonance parameters
2494	bool got_pipi=(fabs(m1-m2)>0.01)?false:true;
2495	double M_sq_R=0.;
2496	if (got_pipi){ // f0(980)
2497	M_sq_R=0.9783*0.9783;
2498	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2499	double temp=4.m1sqm2sq;
2500	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2501	/(4.*M_sq_R));
2502	double partial_width=0.05; //?? // guess from note in pdg
2503	gR=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2504
2505	gsq_rho_S_gamma=0.159; // GeV^-2
2506	gsq_omega_S_gamma=(1./9.)*gsq_rho_S_gamma;
2507	}
2508	else{ // a0(980)
2509	M_sq_R=0.9825*0.9825;
2510	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2511	double temp=4.m1sqm2sq;
2512	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2513	/(4.*M_sq_R));
2514	double partial_width=0.06; //?? // guess from note in pdg
2515	gR=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2516	gsq_rho_S_gamma=0.02537;
2517	gsq_omega_S_gamma=9.*gsq_rho_S_gamma;
2518	}
2519
2520	// Momenta of incoming photon and outgoing S and proton in cm frame
2521	double p_gamma=(s-m_p_sq)/(2.*Ecm);
2522	double E_S=(s+M_sq_R-m_p_sq)/(2.*Ecm);
2523	double p_S=sqrt(E_S*E_S-M_sq_R);
2524
2525	// Momentum transfer t
2526	double p_diff=p_gamma-p_S;
2527	double t0=M_sq_RM_sq_R/(4.s)-p_diff*p_diff;
2528
2529	// differential cross section
2530	double sum=0.;
2531	double t_old=t0;
2532	double t_array[1000];
2533	double xsec_array[1000];
2534	for (unsigned int k=0;k<1000;k++){
2535	double theta_cm=M_PI3.14159265358979323846*double(k)/1000.;
2536	double sin_theta_over_2=sin(0.5*theta_cm);
2537	double t=t0-4.p_gammap_Ssin_theta_over_2sin_theta_over_2;
2538	double xsec=0.;
2539	for (unsigned int j=0;j<1000;j++){
2540	double mass=m1_plus_m2+dm*double(j);
2541	double M_sq=mass*mass;
2542	GetResonanceParameters(m1,m2,M_sq,M_sq_R,ReB,ImB);
2543	xsec+=dm*CrossSection(m1,m2,M_sq_R,s,t,gR,ReB,ImB,gsq_rho_S_gamma,gsq_omega_S_gamma);
2544	}
2545	t_array[k]=-t;
2546	xsec_array[k]=4.M_PI3.141592653589793238461000.*xsec;
2547
2548	sum-=4.M_PI3.141592653589793238461000.xsec(t-t_old);
2549	t_old=t;
2550	}
2551
2552	double m_array[1000];
2553	double xsec_array2[1000];
2554	for (unsigned int j=0;j<1000;j++){
2555	double mass=m1_plus_m2+dm*double(j);
2556	m_array[j]=mass;
2557	double M_sq=mass*mass;
2558	GetResonanceParameters(m1,m2,M_sq,M_sq_R,ReB,ImB);
2559	t_old=t0;
2560	double xsec=0.;
2561	for (unsigned int k=0;k<1000;k++){
2562	double theta_cm=M_PI3.14159265358979323846*double(k)/1000.;
2563	double sin_theta_over_2=sin(0.5*theta_cm);
2564	double t=t0-4.p_gammap_Ssin_theta_over_2sin_theta_over_2;
2565	xsec+=(t_old-t)*CrossSection(m1,m2,M_sq_R,s,t,gR,ReB,ImB,gsq_rho_S_gamma,gsq_omega_S_gamma);
2566	t_old=t;
2567	}
2568	xsec_array2[j]=4.M_PI3.141592653589793238461000.*xsec;
2569	}
2570	double xsec_array3[1000];
2571	if (got_pipi==true){
2572	gsq_rho_S_gamma=0.094;
2573	gsq_omega_S_gamma=(1./9.)*gsq_rho_S_gamma;
2574	}
2575	else{
2576	gsq_rho_S_gamma=0.0054;
2577	gsq_omega_S_gamma=9.*gsq_rho_S_gamma;
2578	}
2579	for (unsigned int j=0;j<1000;j++){
2580	double mass=m1_plus_m2+dm*double(j);
2581	m_array[j]=mass;
2582	double M_sq=mass*mass;
2583	if (got_pipi==false){ // a0(1450)
2584	double m_a1450=1.474;
2585	M_sq_R=m_a1450*m_a1450;
2586	width=0.265;
2587	ReB=M_sq_R-M_sq;
2588	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2589	double temp=4.m1sqm2sq;
2590	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2591	/(4.*M_sq_R));
2592	double partial_width=0.02*width;// estimate using pdg(2016)
2593	gR=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2594	ImB=width*sqrt(M_sq);
2595	}
2596	else{ // f0(1370)
2597	double m_f1370=1.25 ; //guess
2598	M_sq_R=m_f1370*m_f1370;
2599	width=0.5; // estimate, top of PDG range
2600	ReB=M_sq_R-M_sq;
2601	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2602	double temp=4.m1sqm2sq;
2603	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2604	/(4.*M_sq_R));
2605	double partial_width=0.26/3.*width;// estimate using pdg(2016)
2606	gR=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2607	ImB=width*sqrt(M_sq);
2608	}
2609	t_old=t0;
2610	double xsec=0.;
2611	for (unsigned int k=0;k<1000;k++){
2612	double theta_cm=M_PI3.14159265358979323846*double(k)/1000.;
2613	double sin_theta_over_2=sin(0.5*theta_cm);
2614	double t=t0-4.p_gammap_Ssin_theta_over_2sin_theta_over_2;
2615	xsec+=(t_old-t)*CrossSection(m1,m2,M_sq_R,s,t,gR,ReB,ImB,gsq_rho_S_gamma,gsq_omega_S_gamma);
2616	t_old=t;
2617	}
2618	xsec_array3[j]=4.M_PI3.141592653589793238461000.*xsec;
2619	}
2620	if (got_pipi){
2621	// f0(500)
2622	gsq_rho_S_gamma=0.1;
2623	gsq_omega_S_gamma=(1./9)*gsq_rho_S_gamma;
2624	double xsec_array4[1000];
2625	for (unsigned int j=0;j<1000;j++){
2626	double mass=m1_plus_m2+dm*double(j);
2627	m_array[j]=mass;
2628	double M_sq=mass*mass;
2629	double mf500=0.6;
2630	M_sq_R=mf500*mf500;
2631	width=1.;
2632	ReB=M_sq_R-M_sq;
2633	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2634	double temp=4.m1sqm2sq;
2635	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2636	/(4.*M_sq_R));
2637	double partial_width=width/3.;
2638	gR=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2639	ImB=width*sqrt(M_sq);
2640	t_old=t0;
2641	double xsec=0.;
2642	for (unsigned int k=0;k<1000;k++){
2643	double theta_cm=M_PI3.14159265358979323846*double(k)/1000.;
2644	double sin_theta_over_2=sin(0.5*theta_cm);
2645	double t=t0-4.p_gammap_Ssin_theta_over_2sin_theta_over_2;
2646	xsec+=(t_old-t)*CrossSection(m1,m2,M_sq_R,s,t,gR,ReB,ImB,gsq_rho_S_gamma,gsq_omega_S_gamma);
2647	t_old=t;
2648	}
2649	xsec_array4[j]=4.M_PI3.141592653589793238461000.*xsec;
2650	}
2651	TGraph *Gxsec4=new TGraph(1000,m_array,xsec_array4);
2652	Gxsec4->GetXaxis()->SetTitle("M [GeV]");
2653	Gxsec4->GetYaxis()->SetTitle("d#sigma/dM [nb/GeV]");
2654	Gxsec4->Write("Cross section d#sigma/dM");
2655	}
2656
2657	TGraph *Gxsec=new TGraph(1000,t_array,xsec_array);
2658	Gxsec->GetXaxis()->SetTitle("-t [GeV^{2}]");
2659	Gxsec->GetYaxis()->SetTitle("d#sigma/dt [nb/GeV^{2}]");
2660	Gxsec->Write("Cross section d#sigma/dt");
2661	TGraph *Gxsec2=new TGraph(1000,m_array,xsec_array2);
2662	Gxsec2->GetXaxis()->SetTitle("M [GeV]");
2663	Gxsec2->GetYaxis()->SetTitle("d#sigma/dM [nb/GeV]");
2664	Gxsec2->Write("Cross section d#sigma/dM");
2665	TGraph *Gxsec3=new TGraph(1000,m_array,xsec_array3);
2666	Gxsec3->GetXaxis()->SetTitle("M [GeV]");
2667	Gxsec3->GetYaxis()->SetTitle("d#sigma/dM [nb/GeV]");
2668	Gxsec3->Write("Cross section d#sigma/dM");
2669
2670	cout << "Total a0(980) or f0(980) cross section at " << Egamma << " GeV = "<< sum
2671	<< " nano-barns"<<endl;
2672	}
2673
2674
2675	//-----------
2676	// main
2677	//-----------
2678	int main(int narg, char *argv[])
2679	{
2680	ParseCommandLineArguments(narg, argv);
2681
2682	// open ROOT file
2683	string rootfilename="scalar_gen.root";
2684	TFile *rootfile=new TFile(rootfilename.c_str(),"RECREATE",
2685	"Produced by genScalarRegge");
2686
2687	// open HDDM file
2688	s_iostream_t *file = init_s_HDDM(output_file_name);
2689
2690
2691	// Initialize random number generator
2692	TRandom3 *myrand=new TRandom3(0);// If seed is 0, the seed is automatically computed via a TUUID object, according to TRandom3 documentation
2693
2694	// Fixed target
2695	TLorentzVector target(0.,0.,0.,m_p);
2696
2697	//----------------------------------------------------------------------------
2698	// Get production (Egamma range) and decay parameters from input file
2699	//----------------------------------------------------------------------------
2700
2701	// Start reading the input file
2702	ifstream infile(input_file_name);
2703	if (!infile.is_open()){
2704	cerr << "Input file missing! Exiting..." <<endl;
2705	exit(-1);
2706	}
2707
2708	// Get beam properties configuration file
2709	string comment_line;
2710	getline(infile,comment_line);
2711	string beamConfigFile;
2712	infile >> beamConfigFile;
2713	infile.ignore(); // ignore the '\n' at the end of this line
2714
2715	cout << "Photon beam configuration file " << beamConfigFile.data() << endl;
2716
2717	// Set number of decay particles
2718	int num_decay_particles=2;
2719	// Set up vectors of particle ids and 4-vectors
2720	int last_index=num_decay_particles;
2721	int num_final_state_particles=num_decay_particles+1;
2722	vector<TLorentzVector>particle_vectors(num_final_state_particles);
2723	vector<Particle_t>particle_types(num_final_state_particles);
2724	double *decay_masses =new double[num_decay_particles];
2725	particle_types[last_index]=Proton;
2726
2727	// GEANT ids of decay particles
2728	getline(infile,comment_line);
2729	cout << "Particle id's of decay particles =";
2730	for (int k=0;k<num_decay_particles;k++){
2731	int ipart;
2732	infile >> ipart;
2733	cout << " " << ipart;
2734	particle_types[k]=(Particle_t)ipart;
2735	decay_masses[k]=ParticleMass(particle_types[k]);
2736	}
2737	infile.ignore(); // ignore the '\n' at the end of this line
2738	cout << endl;
2739
2740	// Parameters for regge cuts
2741	getline(infile,comment_line);
2742	cout << "Regge cut parameters =";
2743	for (int k=0;k<5;k++){
2744	infile >> regge_cuts[k];
2745	cout << " " << regge_cuts[k];
2746	}
2747	infile.ignore(); // ignore the '\n' at the end of this line
2748	cout << endl;
2749
2750	// processes to simulate
2751	int num_processes=5;
2752	getline(infile,comment_line);
2753	int *generate=new int[num_processes];
2754	for (int k=0;k<num_processes;k++){
2755	infile >> generate[k];
2756	}
2757	infile.ignore(); // ignore the '\n' at the end of this line
2758
2759	// Parameters for regge cuts
2760	getline(infile,comment_line);
2761	double phase[10];
2762	cout << " Interference phases =";
2763	for (int k=0;k<9;k++){
2764	infile >> phase[k];
2765	cout << " " << phase[k];
2766	}
2767	infile.ignore(); // ignore the '\n' at the end of this line
2768	cout << endl;
2769	infile.close();
2770
2771	// Create some diagonistic histograms
2772	CreateHistograms(beamConfigFile);
2773
2774	// Make a TGraph of the cross section at a fixed beam energy
2775	GraphCrossSection(decay_masses[0],decay_masses[1]);
2776
2777	// Set up some variables for cross section calculation
2778	// masses of decay particles
2779	double m1=decay_masses[0];
2780	double m2=decay_masses[1];
2781	bool got_pipi=(fabs(m1-m2)>0.01)?false:true;
2782
2783	double mymax=0.;
2784	//----------------------------------------------------------------------------
2785	// Event generation loop
2786	//----------------------------------------------------------------------------
2787	for (int i=1;i<=Nevents;i++){
2788	// photon beam
2789	double Egamma=0.;
2790	TLorentzVector beam;
2791
2792	// Maximum value for cross section
2793	double xsec_max=(got_pipi)?10.0:3.5;
2794	double xsec=0.,xsec_test=0.;
2795
2796	// Polar angle in center of mass frame
2797	double theta_cm=0.;
2798
2799	// Scalar momentum in cm
2800	double p_S=0.;
2801
2802	// Mass squared of resonance
2803	double M_sq=0.;
2804
2805	// Transfer 4-momentum;
2806	double t=0.;
2807
2808	// vertex position at target
2809	float vert[4]={0.,0.,0.,0.};
2810
2811	// masses of decay particles
2812	double m1sq=m1*m1;
2813	double m2sq=m2*m2;
2814	double m1sq_plus_m2sq=m1sq+m2sq;
2815
2816	// Coupling constants for f0(500)
2817	double gsq_rho_f500_gamma=0.22;
2818	double gsq_omega_f500_gamma=(1./9)*gsq_rho_f500_gamma;
2819	// Coupling constants: Donnachie and Kalashnikova (2008) scenario IV
2820	double gsq_rho_S_gamma=0.02537;
2821	double gsq_omega_S_gamma=0.2283;
2822	double gsq_rho_f1370_gamma=0.094;
2823	double gsq_omega_f1370_gamma=(1./9.)*gsq_rho_f1370_gamma;
2824	double gsq_rho_a1450_gamma=0.0054;
2825	double gsq_omega_a1450_gamma=9.*gsq_rho_a1450_gamma;
2826
2827	// use the rejection method to produce S's based on the cross section
2828	do{
2829	// First generate a beam photon using bremsstrahlung spectrum
2830	Egamma = cobrems_vs_E->GetRandom();
2831
2832	// CM energy
2833	double s=m_p(m_p+2.Egamma);
2834	double Ecm=sqrt(s);
2835
2836	// Momenta of incoming photon and outgoing S and proton in cm frame
2837	double p_gamma=(s-m_p_sq)/(2.*Ecm);
2838
2839	// Mass of two-meson system
2840	double m1_plus_m2=m1+m2;
2841	double m_max=m_p(sqrt(1.+2.Egamma/m_p)-1.);
2842	double M=m1_plus_m2+myrand->Uniform(m_max-m1_plus_m2);
2843	M_sq=M*M;
2844
2845	// Momentum and energy of two-meson system
2846	double E_S=(s+M_sq-m_p_sq)/(2.*Ecm);
2847	p_S=sqrt(E_S*E_S-M_sq);
2848
2849	// Momentum transfer t
2850	double p_diff=p_gamma-p_S;
2851	double t0=M_sqM_sq/(4.s)-p_diff*p_diff;
2852	double sin_theta_over_2=0.;
2853	t=t0;
	Value stored to 't' is never read
2854
2855	// Generate cos(theta) with a uniform distribution and compute the cross
2856	// section at this value
2857	double cos_theta_cm=-1.0+myrand->Uniform(2.);
2858	theta_cm=acos(cos_theta_cm);
2859	// compute t
2860	sin_theta_over_2=sin(0.5*theta_cm);
2861	t=t0-4.p_gammap_Ssin_theta_over_2sin_theta_over_2;
2862
2863	// Generate phi using uniform distribution
2864	double phi_cm=myrand->Uniform(2.*M_PI3.14159265358979323846);
2865
2866	// beam 4-vector (ignoring px and py, which are extremely small)
2867	beam.SetXYZT(0.,0.,Egamma,Egamma);
2868
2869	// Velocity of the cm frame with respect to the lab frame
2870	TVector3 v_cm=(1./(Egamma+m_p))*beam.Vect();
2871	// Four-momentum of the S in the CM frame
2872	double pt=p_S*sin(theta_cm);
2873	TLorentzVector S4(ptcos(phi_cm),ptsin(phi_cm),p_S*cos(theta_cm),
2874	sqrt(p_S*p_S+M_sq));
2875	// S4.Print();
2876
2877	//Boost the S 4-momentum into the lab
2878	S4.Boost(v_cm);
2879	// S4.Print();
2880
2881
2882	// Compute the 4-momentum for the recoil proton
2883	TLorentzVector proton4=beam+target-S4;
2884
2885	// Generate decay of S according to phase space
2886	TGenPhaseSpace phase_space;
2887	phase_space.SetDecay(S4,num_decay_particles,decay_masses);
2888	double weight=0.,rand_weight=1.;
2889	do{
2890	weight=phase_space.Generate();
2891	rand_weight=myrand->Uniform(1.);
2892	}
2893	while (rand_weight>weight);
2894
2895	// Gather the particles in the reaction and write out event in hddm format
2896	particle_vectors[last_index]=proton4;
2897	for (int j=0;j<num_decay_particles;j++){
2898	particle_vectors[j]=*phase_space.GetDecay(j);
2899	}
2900
2901	//Resonance parameters
2902	double ReB=0.,ImB=0,gR=0.;
2903	double gR_T=0., ImB_T=0., ReB_T=0.;
2904	double gRf500=0.,ImBf500=0.,ReBf500=0.;
2905	double gRf1370=0.,ImBf1370=0.,ReBf1370=0.;
2906	double gRa1450=0.,ImBa1450=0.,ReBa1450=0.;
2907
2908	// Cross section
2909	xsec=0.;
2910
2911	// f0(600)
2912	if (got_pipi && generate[0]){
2913	double m_Sigma=0.6;
2914	double M_sq_R=m_Sigma*m_Sigma;
2915	width=1.0;
2916	ReBf500=M_sq_R-M_sq;
2917	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2918	double temp=4.m1sqm2sq;
2919	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2920	/(4.*M_sq_R));
2921	double partial_width=width/3.;
2922	gRf500=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2923	ImBf500=width*sqrt(M_sq);
2924
2925	xsec+=CrossSection(m1,m2,M_sq,s,t,gRf500,ReBf500,ImBf500,
2926	gsq_rho_f500_gamma,gsq_omega_f500_gamma);
2927
2928	}
2929	// f0(1370)
2930	if (got_pipi && generate[2]){
2931	double m_f1370=1.25; // guess
2932	double M_sq_R=m_f1370*m_f1370;
2933	width=0.5; // estimate, top of PDG range
2934	ReBf1370=M_sq_R-M_sq;
2935	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2936	double temp=4.m1sqm2sq;
2937	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2938	/(4.*M_sq_R));
2939	// partial width from Bugg(2007):arxiv.org/pdf/0706.1341.pdf, table 2
2940	double partial_width=0.325/3.;
2941	gRf1370=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2942	ImBf1370=width*sqrt(M_sq);
2943
2944	xsec+=CrossSection(m1,m2,M_sq,s,t,gRf1370,ReBf1370,ImBf1370,
2945	gsq_rho_f1370_gamma,gsq_omega_f1370_gamma);
2946
2947	}
2948	// Interference between f0(500) and f0(1370)
2949	if (got_pipi && generate[0] && generate[2]){
2950	xsec+=CrossSection(m1,m2,M_sq,s,t,gRf1370,ReBf1370,ImBf1370,
2951	gsq_rho_f1370_gamma,
2952	gsq_omega_f1370_gamma,true,gRf500,ReBf500,ImBf500,
2953	gsq_rho_f500_gamma,gsq_omega_f500_gamma,
2954	phase[9]);
2955	}
2956
2957	// a0(1450)
2958	if (got_pipi==false && generate[2]){
2959	double m_a1450=1.448;// Bugg:arXiv:0808.2706v2,Table 2
2960	double M_sq_R=m_a1450*m_a1450;
2961	width=0.192; // Bugg:arXiv:0808.2706v2,Table 2
2962	ReBa1450=M_sq_R-M_sq;
2963	double MRsq_minus_m1sq_m2sq=M_sq_R-m1sq_plus_m2sq;
2964	double temp=4.m1sqm2sq;
2965	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2966	/(4.*M_sq_R));
2967	double partial_width=0.0237;// Bugg:arXiv:0808.2706v2,Table 2
2968	gRa1450=sqrt(8.M_PI3.14159265358979323846M_sq_R*partial_width/qR);
2969	ImBa1450=width*sqrt(M_sq);
2970
2971	xsec+=CrossSection(m1,m2,M_sq,s,t,gRa1450,ReBa1450,ImBa1450,
2972	gsq_rho_a1450_gamma,gsq_omega_a1450_gamma);
2973
2974	}
2975	// f0(980)/a0(980)
2976	if (generate[1]){
2977	double my_msq_R=0.9783*0.9783;
2978	if (got_pipi){ // f0(980)
2979	double MRsq_minus_m1sq_m2sq=my_msq_R-m1sq_plus_m2sq;
2980	double temp=4.m1sqm2sq;
2981	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2982	/(4.*my_msq_R));
2983	double partial_width=0.05; //?? // guess from note in pdg
2984	gR=sqrt(8.M_PI3.14159265358979323846my_msq_R*partial_width/qR);
2985	gsq_rho_S_gamma=0.159; // GeV^-2
2986	gsq_omega_S_gamma=(1./9.)*gsq_rho_S_gamma;
2987	}
2988	else{ // a0(980)
2989	my_msq_R=0.9825*0.9825;
2990	double MRsq_minus_m1sq_m2sq=my_msq_R-m1sq_plus_m2sq;
2991	double temp=4.m1sqm2sq;
2992	double qR=sqrt((MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp)
2993	/(4.*my_msq_R));
2994	double partial_width=0.06; //?? guess from note in pdg
2995	gR=sqrt(8.M_PI3.14159265358979323846my_msq_R*partial_width/qR);
2996	gsq_rho_S_gamma=0.02537;
2997	gsq_omega_S_gamma=9.*gsq_rho_S_gamma;
2998	}
2999	GetResonanceParameters(m1,m2,M_sq,my_msq_R,ReB,ImB);
3000	xsec+=CrossSection(m1,m2,M_sq,s,t,gR,ReB,ImB,gsq_rho_S_gamma,
3001	gsq_omega_S_gamma);
3002	// Interference between f0(1370) and f0(980)
3003	if (got_pipi && generate[2]){
3004	xsec+=CrossSection(m1,m2,M_sq,s,t,gR,ReB,ImB,gsq_rho_S_gamma,
3005	gsq_omega_S_gamma,true,gRf1370,ReBf1370,ImBf1370,
3006	gsq_rho_f1370_gamma,gsq_omega_f1370_gamma,
3007	phase[1]);
3008	}
3009	// Interference between f0(500) and f0(980)
3010	if (got_pipi && generate[0]){
3011	xsec+=CrossSection(m1,m2,M_sq,s,t,gR,ReB,ImB,gsq_rho_S_gamma,
3012	gsq_omega_S_gamma,true,gRf500,ReBf500,ImBf500,
3013	gsq_rho_f500_gamma,gsq_omega_f500_gamma,
3014	phase[0]);
3015	}
3016	// Interference between a0(1450) and a0(980)
3017	if (got_pipi==false && generate[2]){
3018	xsec+=CrossSection(m1,m2,M_sq,s,t,gR,ReB,ImB,gsq_rho_S_gamma,
3019	gsq_omega_S_gamma,true,gRa1450,ReBa1450,ImBa1450,
3020	gsq_rho_a1450_gamma,gsq_omega_a1450_gamma,
3021	phase[1]);
3022	}
3023	}
3024	if (generate[4]){ // non-resonant background
3025	xsec+=BackgroundCrossSection(beam,particle_types,particle_vectors);
3026
3027	if (generate[1]){ // interference with resonant signal
3028	xsec+=InterferenceCrossSection(beam,particle_types,particle_vectors,
3029	gR,ReB,ImB,gsq_rho_S_gamma,
3030	gsq_omega_S_gamma,phase[3]);
3031	}
3032	if (got_pipi){
3033	if (generate[0]){ // interference with f0(500)
3034	xsec+=InterferenceCrossSection(beam,particle_types,particle_vectors,
3035	gRf500,ReBf500,ImBf500,
3036	gsq_rho_f500_gamma,
3037	gsq_omega_f500_gamma,phase[7]);
3038	}
3039	if (generate[2]){ // interference with f0(1370)
3040	xsec+=InterferenceCrossSection(beam,particle_types,particle_vectors,
3041	gRf1370,ReBf1370,ImBf1370,
3042	gsq_rho_f1370_gamma,
3043	gsq_omega_f1370_gamma,phase[8]);
3044	}
3045	}
3046	else{
3047	if (generate[2]){ // interference with a0(1450)
3048	xsec+=InterferenceCrossSection(beam,particle_types,particle_vectors,
3049	gRa1450,ReBa1450,ImBa1450,
3050	gsq_rho_a1450_gamma,
3051	gsq_omega_a1450_gamma,0.);
3052	}
3053
3054	}
3055	}
3056	if (generate[3]){ // Tensor background
3057	double m_T=1.275;
3058	double Gamma_T=0.185;
3059	if (!got_pipi){
3060	Gamma_T=0.107;
3061	m_T=1.3183;
3062	}
3063	double M_sq_R_T=m_T*m_T;
3064	ReB_T=M_sq_R_T-M_sq;
3065	double Msq_minus_m1sq_m2sq=M_sq-m1sq_plus_m2sq;
3066	double MRsq_minus_m1sq_m2sq=M_sq_R_T-m1sq_plus_m2sq;
3067	double temp=4.m1sqm2sq;
3068	double q_over_qR
3069	=m_T/Msqrt((Msq_minus_m1sq_m2sqMsq_minus_m1sq_m2sq-temp)
3070	/(MRsq_minus_m1sq_m2sq*MRsq_minus_m1sq_m2sq-temp));
3071	double q_over_qR_5=pow(q_over_qR,5);
3072	if (got_pipi){ // f2(1270)
3073	double partial_width=0.85(1./3.)M_sq_R_T*q_over_qR_5/M_sq;
3074	gR_T=sqrt(8.M_PI3.14159265358979323846M_sq_R_Tpartial_widthq_over_qR);
3075	ImB_T=MGamma_T(0.85q_over_qR_5M_sq_R_T/M_sq+0.15);
3076	}
3077	else { // a2(1320)
3078	double partial_width=0.155M_sq_R_Tq_over_qR_5/M_sq;
3079	gR_T=sqrt(8.M_PI3.14159265358979323846M_sq_R_Tpartial_widthq_over_qR);
3080	ImB_T=MGamma_T(0.145q_over_qR_5M_sq_R_T/M_sq+0.855);
3081	}
3082	xsec+=TensorCrossSection(beam,particle_types,particle_vectors,
3083	gR_T,ReB_T,ImB_T);
3084	if (generate[1]){ // interference with a0(980)/f0(980)
3085	xsec+=TensorScalarInterference(beam,particle_types,particle_vectors,
3086	gR_T,ReB_T,ImB_T,gR,ReB,ImB,
3087	sqrt(gsq_omega_S_gamma),
3088	sqrt(gsq_rho_S_gamma),phase[2]);
3089	}
3090	// interference with f0(600)
3091	if (got_pipi && generate[0]){
3092	xsec+=TensorScalarInterference(beam,particle_types,particle_vectors,
3093	gR_T,ReB_T,ImB_T,gRf500,ReBf500,
3094	ImBf500,
3095	sqrt(gsq_omega_f500_gamma),
3096	sqrt(gsq_rho_f500_gamma),phase[5]);
3097	}
3098	// interference with f0(1370)
3099	if (got_pipi && generate[2]){
3100	xsec+=TensorScalarInterference(beam,particle_types,particle_vectors,
3101	gR_T,ReB_T,ImB_T,gRf1370,ReBf1370,
3102	ImBf1370,
3103	sqrt(gsq_omega_f1370_gamma),
3104	sqrt(gsq_rho_f1370_gamma),phase[6]);
3105	}
3106	// interference with a0(1450)
3107	if (got_pipi==false && generate[2]){
3108	xsec+=TensorScalarInterference(beam,particle_types,particle_vectors,
3109	gR_T,ReB_T,ImB_T,gRa1450,ReBa1450,
3110	ImBa1450,
3111	sqrt(gsq_omega_a1450_gamma),
3112	sqrt(gsq_rho_a1450_gamma),phase[6]);
3113	}
3114
3115	if (generate[4]){ //interference with background wave
3116	xsec+=TensorBackgroundInterference(beam,particle_types,
3117	particle_vectors,
3118	gR_T,ReB_T,ImB_T,phase[4]);
3119	}
3120
3121	}
3122	if (xsec>mymax ) mymax=xsec;
3123
3124
3125
3126	// Generate a test value for the cross section
3127	xsec_test=myrand->Uniform(xsec_max);
3128	}
3129	while (xsec_test>xsec);
3130
3131	// Other diagnostic histograms
3132	thrown_t->Fill(-t);
3133	thrown_Egamma->Fill(Egamma);
3134	thrown_theta_vs_p->Fill(particle_vectors[last_index].P(),
3135	180./M_PI3.14159265358979323846*particle_vectors[last_index].Theta());
3136	thrown_mass->Fill((particle_vectors[0]+particle_vectors[1]).M());
3137	thrown_mass_vs_E->Fill(Egamma,(particle_vectors[0]+particle_vectors[1]).M());
3138
3139	// Randomly generate z position in target
3140	vert[2]=zmin+myrand->Uniform(zmax-zmin);
3141
3142	WriteEvent(i,beam,vert,particle_types,particle_vectors,file);
3143
3144	if ((i%(Nevents/10))==0) cout << 100.*double(i)/double(Nevents) << "\% done" << endl;
3145	}
3146
3147
3148	// Write histograms and close root file
3149	rootfile->Write();
3150	rootfile->Close();
3151
3152	// Close HDDM file
3153	close_s_HDDM(file);
3154	cout<<endl<<"Closed HDDM file"<<endl;
3155	cout<<" "<<Nevents<<" event written to "<<output_file_name<<endl;
3156
3157	cout << mymax << endl;
3158
3159	return 0;
3160	}
3161
3162