65 {
66loop:
68
70 EvtVector4R pv, ps, ppr;
71
77
78 EvtHelSys angles( ppr, pv );
79 double theta = angles.getHelAng( 1 );
80 double phi = angles.getHelAng( 2 );
81 double gamma = 0;
82
83 double m_b0 =
86
88 double M2 = pow( m_M, 2.0 );
89 double b2_p = pow( ( m_b0 +
m_b1 ), 2.0 );
90 double b2_m = pow( (
m_b1 - m_b0 ), 2.0 );
91
92 if ( M2 - b2_p <= 0 ) goto loop;
93
94 double P_c = sqrt( ( M2 - b2_p ) * ( M2 - b2_m ) / ( 4.0 * M2 ) );
95 double Q = M2 - b2_p;
96 double F1 = -0.54 / ( m_b0 * b2_p );
97 double F2 = -0.54 * (
m_b1 - m_b0 ) / ( m_b0 * ( m_b0 +
m_b1 ) );
98 double F3 = -0.89 / m_b0;
99 double H1, H0, Hm1;
100
101
102
104 {
106 H0 = 1.0;
108 Hm1 = H1;
109 }
110 else
111 {
112 H1 = -1.15 * P_c * m_M * ( -M2 * Q * F1 - Q * F2 + ( M2 - m_b0 * ( m_b0 +
m_b1 ) ) * F3 ) /
113 ( sqrt( Q ) *
m_b1 );
114 H0 =
115 -0.82 * P_c * M2 *
116 ( -(
m_b1 - m_b0 ) * Q * F1 /
m_b1 - Q * F2 / (
m_b1 * (
m_b1 - m_b0 ) ) + 2.0 * F3 ) /
117 sqrt( Q );
118 Hm1 = 2 * P_c * m_M * ( m_b0 +
m_b1 ) * F3 / sqrt( Q );
119 alpha = ( H1 * H1 + Hm1 * Hm1 - 2 * H0 * H0 ) / ( H1 * H1 + Hm1 * Hm1 + 2 * H0 * H0 );
120 }
121
122
123
124 vertex( 0, 0, 0,
Djmn( 1, 1, 0, phi, theta, gamma ) * H0 );
125 vertex( 0, 0, 1,
Djmn( 1, 1, -1, phi, theta, gamma ) * Hm1 );
127 vertex( 0, 0, 3,
Djmn( 1, 1, 1, phi, theta, gamma ) * H1 );
128 vertex( 0, 1, 0,
Djmn( 1, 1, -1, phi, theta, gamma ) * H1 );
130 vertex( 0, 1, 2,
Djmn( 1, 1, 1, phi, theta, gamma ) * Hm1 );
131 vertex( 0, 1, 3,
Djmn( 1, 1, 0, phi, theta, gamma ) * H0 );
132
133 vertex( 1, 0, 0,
Djmn( 1, -1, 0, phi, theta, gamma ) * H0 );
134 vertex( 1, 0, 1,
Djmn( 1, -1, -1, phi, theta, gamma ) * Hm1 );
136 vertex( 1, 0, 3,
Djmn( 1, -1, 1, phi, theta, gamma ) * H1 );
137 vertex( 1, 1, 0,
Djmn( 1, -1, -1, phi, theta, gamma ) * H1 );
139 vertex( 1, 1, 2,
Djmn( 1, -1, 1, phi, theta, gamma ) * Hm1 );
140 vertex( 1, 1, 3,
Djmn( 1, -1, 0, phi, theta, gamma ) * H0 );
141
142
143 vertex( 2, 0, 0,
Djmn( 1, 0, 0, phi, theta, gamma ) * H0 );
144 vertex( 2, 0, 1,
Djmn( 1, 0, -1, phi, theta, gamma ) * Hm1 );
146 vertex( 2, 0, 3,
Djmn( 1, 0, 1, phi, theta, gamma ) * H1 );
147 vertex( 2, 1, 0,
Djmn( 1, 0, -1, phi, theta, gamma ) * H1 );
149 vertex( 2, 1, 2,
Djmn( 1, 0, 1, phi, theta, gamma ) * Hm1 );
150 vertex( 2, 1, 3,
Djmn( 1, 0, 0, phi, theta, gamma ) * H0 );
151
152 return;
153}
EvtComplex Djmn(int j, int m, int n, double phi, double theta, double gamma)
*********DOUBLE PRECISION m_pi INTEGER m_lenwt !max no of aux weights INTEGER m_phmax !maximum photon multiplicity ISR FSR *DOUBLE COMPLEX m_Pauli4 DOUBLE COMPLEX m_AmpBorn DOUBLE COMPLEX m_AmpBoxy DOUBLE COMPLEX m_AmpBorn1 DOUBLE COMPLEX m_AmpBorn2 DOUBLE COMPLEX m_AmpExpo2p DOUBLE COMPLEX m_Rmat DOUBLE COMPLEX m_BoxGZut !DOUBLE COMPLEX m_F1finPair2 !DOUBLE PRECISION m_Vcut DOUBLE PRECISION m_Alfinv DOUBLE PRECISION m_Lorin1 DOUBLE PRECISION m_Lorin2 DOUBLE PRECISION m_b1
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
void vertex(const EvtComplex &)
static double getMass(EvtId i)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)