70 {
71
72loop:
74
75 EvtParticle *
v, *s1, *s2;
76 EvtVector4R pv, ps;
77 EvtVector4R Lpt, xp;
78
86
87
88 double costheta;
89 if ( Lpt.
d3mag() < 0.0001 ) { costheta = xp.
get( 3 ) / xp.
d3mag(); }
90 else
91 {
92 costheta = ( xp.
get( 1 ) * Lpt.
get( 1 ) + xp.
get( 2 ) * Lpt.
get( 2 ) +
93 xp.
get( 3 ) * Lpt.
get( 3 ) ) /
95 }
97 double c0, c2, c4, c6, c8, c10;
98 c0 = 0;
99 c2 = 0;
100 c4 = 0;
101 c6 = 0;
102 c8 = 0;
103 c10 = 0;
104 if ( narg == 2 )
105 {
108 }
109 else if ( narg == 3 )
110 {
114 }
115 else if ( narg == 4 )
116 {
121 }
122 else if ( narg == 5 )
123 {
129 }
130 else if ( narg == 6 )
131 {
138 }
139
140 double costheta2 = costheta * costheta;
141 double costheta4 = costheta2 * costheta2;
142 double costheta6 = costheta4 * costheta2;
143 double costheta8 = costheta6 * costheta2;
144 double costheta10 = costheta8 * costheta2;
145
146 double amp1 = ( c0 + c2 * costheta2 + c4 * costheta4 + c6 * costheta6 + c8 * costheta8 +
147 c10 * costheta10 );
148 double a0, a2, a4, a6, a8, a10;
149 double ampflag;
150 if ( c0 < 0 ) {
a0 = 0; }
152 if ( c2 < 0 ) { a2 = 0; }
153 else { a2 = c2; }
154 if ( c4 < 0 ) { a4 = 0; }
155 else { a4 = c4; }
156 if ( c6 < 0 ) { a6 = 0; }
157 else { a6 = c6; }
158 if ( c8 < 0 ) { a8 = 0; }
159 else { a8 = c8; }
160 if ( c10 < 0 ) { a10 = 0; }
161 else { a10 = c10; }
162 ampflag =
a0 + a2 + a4 + a6 + a8 + a10;
163 if ( ampflag <= 0 )
164 {
166 << " The maxium value of amplitude square should be positive, but it is " << ampflag
167 << endl;
168 ::abort();
169 }
170 ampflag = amp1 / ampflag;
172 if ( rd1 >= ampflag ) goto loop;
173 return;
174}
character *LEPTONflag integer iresonances real zeta5 real a0
ostream & report(Severity severity, const char *facility)
**********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
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
EvtVector4R cross(const EvtVector4R &v2)