71 {
72
73loop:
75
77 EvtVector4R pv, ps, ppr;
78
84
85 EvtHelSys angles( ppr, pv );
86 double theta = angles.getHelAng( 1 );
87
88
89
90
92 double c0, c2, c4, c6, c8, c10;
93 c0 = 0;
94 c2 = 0;
95 c4 = 0;
96 c6 = 0;
97 c8 = 0;
98 c10 = 0;
99 if ( narg == 2 )
100 {
103 }
104 else if ( narg == 3 )
105 {
109 }
110 else if ( narg == 4 )
111 {
116 }
117 else if ( narg == 5 )
118 {
124 }
125 else if ( narg == 6 )
126 {
133 }
134
135 double costheta =
cos( theta );
136 double costheta2 = costheta * costheta;
137 double costheta4 = costheta2 * costheta2;
138 double costheta6 = costheta4 * costheta2;
139 double costheta8 = costheta6 * costheta2;
140 double costheta10 = costheta8 * costheta2;
141
142 double amp1 = ( c0 + c2 * costheta2 + c4 * costheta4 + c6 * costheta6 + c8 * costheta8 +
143 c10 * costheta10 );
144 double a0, a2, a4, a6, a8, a10;
145 double ampflag;
146 if ( c0 < 0 ) {
a0 = 0; }
148 if ( c2 < 0 ) { a2 = 0; }
149 else { a2 = c2; }
150 if ( c4 < 0 ) { a4 = 0; }
151 else { a4 = c4; }
152 if ( c6 < 0 ) { a6 = 0; }
153 else { a6 = c6; }
154 if ( c8 < 0 ) { a8 = 0; }
155 else { a8 = c8; }
156 if ( c10 < 0 ) { a10 = 0; }
157 else { a10 = c10; }
158 ampflag =
a0 + a2 + a4 + a6 + a8 + a10;
159 if ( ampflag <= 0 )
160 {
162 << " The maxium value of amplitude square should be positive, but it is " << ampflag
163 << endl;
164 ::abort();
165 }
166 ampflag = amp1 / ampflag;
168 if ( rd1 >= ampflag ) goto loop;
169 return;
170}
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)