69 {
70
71loop:
73
75 EvtVector4R pv, ps, ppr;
76
82
83 EvtHelSys angles( ppr, pv );
84 double theta = angles.getHelAng( 1 );
85
86
87
88
90 double c0, c2, c4, c6, c8, c10;
91 c0 = 0;
92 c2 = 0;
93 c4 = 0;
94 c6 = 0;
95 c8 = 0;
96 c10 = 0;
97 if ( narg == 2 )
98 {
101 }
102 else if ( narg == 3 )
103 {
107 }
108 else if ( narg == 4 )
109 {
114 }
115 else if ( narg == 5 )
116 {
122 }
123 else if ( narg == 6 )
124 {
131 }
132
133 double costheta =
cos( theta );
134 double costheta2 = costheta * costheta;
135 double costheta4 = costheta2 * costheta2;
136 double costheta6 = costheta4 * costheta2;
137 double costheta8 = costheta6 * costheta2;
138 double costheta10 = costheta8 * costheta2;
139
140 double amp1 = ( c0 + c2 * costheta2 + c4 * costheta4 + c6 * costheta6 + c8 * costheta8 +
141 c10 * costheta10 );
142 double a0, a2, a4, a6, a8, a10;
143 double ampflag;
144 if ( c0 < 0 ) {
a0 = 0; }
146 if ( c2 < 0 ) { a2 = 0; }
147 else { a2 = c2; }
148 if ( c4 < 0 ) { a4 = 0; }
149 else { a4 = c4; }
150 if ( c6 < 0 ) { a6 = 0; }
151 else { a6 = c6; }
152 if ( c8 < 0 ) { a8 = 0; }
153 else { a8 = c8; }
154 if ( c10 < 0 ) { a10 = 0; }
155 else { a10 = c10; }
156 ampflag =
a0 + a2 + a4 + a6 + a8 + a10;
157 if ( ampflag <= 0 )
158 {
160 << " The maxium value of amplitude square should be positive, but it is " << ampflag
161 << endl;
162 ::abort();
163 }
164 ampflag = amp1 / ampflag;
166 if ( rd1 >= ampflag ) goto loop;
167 return;
168}
character *LEPTONflag integer iresonances real zeta5 real a0
ostream & report(Severity severity, const char *facility)
double cos(const BesAngle a)
**********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)