79 {
80
81loop:
83
86
92
94 double theta =angles.getHelAng(1);
95
96
97
99 double c0,c2,c4,c6,c8,c10;
100 c0=0;c2=0;c4=0;c6=0;c8=0;c10=0;
101 if(narg==2) {
104 } else if (narg==3){
108 } else if (narg==4){
113 } else if (narg==5){
119 } else if (narg==6){
126 }
127
128 double costheta=
cos(theta);
129 double costheta2= costheta*costheta;
130 double costheta4= costheta2*costheta2;
131 double costheta6= costheta4*costheta2;
132 double costheta8= costheta6*costheta2;
133 double costheta10= costheta8*costheta2;
134
135 double amp1=(c0+c2*costheta2+c4*costheta4+c6*costheta6+c8*costheta8+c10*costheta10);
136 double a0,a2,a4,a6,a8,a10;
137 double ampflag;
138 if(c0<0) {a0=0;}else{a0=c0;}
139 if(c2<0) {a2=0;}else{a2=c2;}
140 if(c4<0) {a4=0;}else{a4=c4;}
141 if(c6<0) {a6=0;}else{a6=c6;}
142 if(c8<0) {a8=0;}else{a8=c8;}
143 if(c10<0) {a10=0;}else{a10=c10;}
144 ampflag=a0+a2+a4+a6+a8+a10;
145 if( ampflag<=0 ) {
146 report(
ERROR,
"EvtGen")<<
" The maxium value of amplitude square should be positive, but it is "<< ampflag<<endl;
147 ::abort();
148 }
149 ampflag = amp1/ampflag;
151 if(rd1>=ampflag) goto loop;
152 return ;
153}
double cos(const BesAngle a)
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)