00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <iostream>
00024 #include "EvtGenBase/EvtGenKine.hh"
00025 #include "EvtGenBase/EvtRandom.hh"
00026 #include "EvtGenBase/EvtVector4R.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtConst.hh"
00029 #include <math.h>
00030 using std::endl;
00031
00032
00033 double EvtPawt(double a,double b,double c)
00034 {
00035 double temp=(a*a-(b+c)*(b+c))*(a*a-(b-c)*(b-c));
00036
00037 if (temp<=0) {
00038
00039
00040
00041
00042 return 0.0;
00043
00044 }
00045
00046 return sqrt(temp)/(2.0*a);
00047 }
00048
00049
00050 double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
00051 double mp )
00052
00053
00054
00055
00056
00057 {
00058
00059 double energy, p3, alpha, beta;
00060
00061 if ( ndaug == 1 ) {
00062 p4[0].set(mass[0],0.0,0.0,0.0);
00063 return 1.0;
00064 }
00065
00066
00067 if ( ndaug == 2 ) {
00068
00069
00070
00071 energy = ( mp*mp + mass[0]*mass[0] -
00072 mass[1]*mass[1] ) / ( 2.0 * mp );
00073
00074 p3 = sqrt( energy*energy - mass[0]*mass[0] );
00075
00076 p4[0].set( energy, 0.0, 0.0, p3 );
00077
00078 energy = mp - energy;
00079 p3 = -1.0*p3;
00080 p4[1].set( energy, 0.0, 0.0, p3 );
00081
00082
00083
00084 alpha = EvtRandom::Flat( EvtConst::twoPi );
00085 beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
00086
00087 p4[0].applyRotateEuler( alpha, beta, -alpha );
00088 p4[1].applyRotateEuler( alpha, beta, -alpha );
00089
00090 return 1.0;
00091 }
00092
00093 if ( ndaug != 2 ) {
00094
00095 double wtmax=0.0;
00096 double pm[5][30],to[4],pmin,pmax,psum,rnd[30];
00097 double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
00098 int i,il,ilr,i1,il1u,il1,il2r,ilu;
00099 int il2=0;
00100
00101 for(i=0;i<ndaug;i++){
00102 pm[4][i]=0.0;
00103 rnd[i]=0.0;
00104 }
00105
00106 pm[0][0]=mp;
00107 pm[1][0]=0.0;
00108 pm[2][0]=0.0;
00109 pm[3][0]=0.0;
00110 pm[4][0]=mp;
00111
00112 to[0]=mp;
00113 to[1]=0.0;
00114 to[2]=0.0;
00115 to[3]=0.0;
00116
00117 psum=0.0;
00118 for(i=1;i<ndaug+1;i++){
00119 psum=psum+mass[i-1];
00120 }
00121
00122 pm[4][ndaug-1]=mass[ndaug-1];
00123
00124 switch (ndaug) {
00125
00126 case 1:
00127 wtmax=1.0/16.0;
00128 break;
00129 case 2:
00130 wtmax=1.0/150.0;
00131 break;
00132 case 3:
00133 wtmax=1.0/2.0;
00134 break;
00135 case 4:
00136 wtmax=1.0/5.0;
00137 break;
00138 case 5:
00139 wtmax=1.0/15.0;
00140 break;
00141 case 6:
00142 wtmax=1.0/15.0;
00143 break;
00144 case 7:
00145 wtmax=1.0/15.0;
00146 break;
00147 case 8:
00148 wtmax=1.0/15.0;
00149 break;
00150 case 9:
00151 wtmax=1.0/15.0;
00152 break;
00153 case 10:
00154 wtmax=1.0/15.0;
00155 break;
00156 case 11:
00157 wtmax=1.0/15.0;
00158 break;
00159 case 12:
00160 wtmax=1.0/15.0;
00161 break;
00162 case 13:
00163 wtmax=1.0/15.0;
00164 break;
00165 case 14:
00166 wtmax=1.0/15.0;
00167 break;
00168 case 15:
00169 wtmax=1.0/15.0;
00170 break;
00171 default:
00172 report(ERROR,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
00173 break;
00174 }
00175
00176 pmax=mp-psum+mass[ndaug-1];
00177
00178 pmin=0.0;
00179
00180 for(ilr=2;ilr<ndaug+1;ilr++){
00181 il=ndaug+1-ilr;
00182 pmax=pmax+mass[il-1];
00183 pmin=pmin+mass[il+1-1];
00184 wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
00185 }
00186
00187
00188
00189 do{
00190
00191 rnd[0]=1.0;
00192 il1u=ndaug-1;
00193
00194 for (il1=2;il1<il1u+1;il1++){
00195 ran=EvtRandom::Flat();
00196 for (il2r=2;il2r<il1+1;il2r++){
00197 il2=il1+1-il2r;
00198 if (ran<=rnd[il2-1]) goto two39;
00199 rnd[il2+1-1]=rnd[il2-1];
00200 }
00201 two39:
00202 rnd[il2+1-1]=ran;
00203 }
00204
00205 rnd[ndaug-1]=0.0;
00206 wt=1.0;
00207 for(ilr=2;ilr<ndaug+1;ilr++){
00208 il=ndaug+1-ilr;
00209 pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
00210 (rnd[il-1]-rnd[il+1-1])*(mp-psum);
00211 wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
00212 }
00213 if (wt>wtmax) {
00214 report(ERROR,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
00215 << ndaug <<" daughters"<<endl;;
00216 }
00217 } while (wt<EvtRandom::Flat(wtmax));
00218
00219 ilu=ndaug-1;
00220
00221 for (il=1;il<ilu+1;il++){
00222 pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
00223 costh=EvtRandom::Flat(-1.0,1.0);
00224 sinth=sqrt(1.0-costh*costh);
00225 phi=EvtRandom::Flat(EvtConst::twoPi);
00226 p[1][il-1]=pa*sinth*cos(phi);
00227 p[2][il-1]=pa*sinth*sin(phi);
00228 p[3][il-1]=pa*costh;
00229 pm[1][il+1-1]=-p[1][il-1];
00230 pm[2][il+1-1]=-p[2][il-1];
00231 pm[3][il+1-1]=-p[3][il-1];
00232 p[0][il-1]=sqrt(pa*pa+mass[il-1]*mass[il-1]);
00233 pm[0][il+1-1]=sqrt(pa*pa+pm[4][il+1-1]*pm[4][il+1-1]);
00234 }
00235
00236 p[0][ndaug-1]=pm[0][ndaug-1];
00237 p[1][ndaug-1]=pm[1][ndaug-1];
00238 p[2][ndaug-1]=pm[2][ndaug-1];
00239 p[3][ndaug-1]=pm[3][ndaug-1];
00240
00241 for (ilr=2;ilr<ndaug+1;ilr++){
00242 il=ndaug+1-ilr;
00243 be[0]=pm[0][il-1]/pm[4][il-1];
00244 be[1]=pm[1][il-1]/pm[4][il-1];
00245 be[2]=pm[2][il-1]/pm[4][il-1];
00246 be[3]=pm[3][il-1]/pm[4][il-1];
00247
00248 for (i1=il;i1<ndaug+1;i1++){
00249 bep=be[1]*p[1][i1-1]+be[2]*p[2][i1-1]+
00250 be[3]*p[3][i1-1]+be[0]*p[0][i1-1];
00251 temp=(p[0][i1-1]+bep)/(be[0]+1.0);
00252 p[1][i1-1]=p[1][i1-1]+temp*be[1];
00253 p[2][i1-1]=p[2][i1-1]+temp*be[2];
00254 p[3][i1-1]=p[3][i1-1]+temp*be[3];
00255 p[0][i1-1]=bep;
00256
00257 }
00258 }
00259
00260 for (ilr=0;ilr<ndaug;ilr++){
00261 p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
00262 }
00263
00264 return 1.0;
00265 }
00266
00267 return 1.0;
00268
00269 }
00270
00271
00272 double EvtGenKine::PhaseSpacePole(double M, double m1, double m2, double m3,
00273 double a,EvtVector4R p4[10])
00274
00275
00276
00277 {
00278
00279
00280
00281
00282 double m12sqmax=(M-m3)*(M-m3);
00283 double m12sqmin=(m1+m2)*(m1+m2);
00284
00285 double m13sqmax=(M-m2)*(M-m2);
00286 double m13sqmin=(m1+m3)*(m1+m3);
00287
00288 double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
00289 double v2= a*(1.0/m12sqmin-1.0/m12sqmax)*(m13sqmax-m13sqmin);
00290
00291 double m12sq,m13sq;
00292
00293 double r=v1/(v1+v2);
00294
00295
00296
00297 double m13min,m13max;
00298
00299 do{
00300
00301 m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
00302
00303 if (r>EvtRandom::Flat()){
00304 m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
00305 }
00306 else{
00307 m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
00308 }
00309
00310
00311 double E3star=(M*M-m12sq-m3*m3)/sqrt(4*m12sq);
00312 double E1star=(m12sq+m1*m1-m2*m2)/sqrt(4*m12sq);
00313 double p3star=sqrt(E3star*E3star-m3*m3);
00314 double p1star=sqrt(E1star*E1star-m1*m1);
00315 m13max=(E3star+E1star)*(E3star+E1star)-
00316 (p3star-p1star)*(p3star-p1star);
00317 m13min=(E3star+E1star)*(E3star+E1star)-
00318 (p3star+p1star)*(p3star+p1star);
00319
00320 }while(m13sq<m13min||m13sq>m13max);
00321
00322
00323 double E2=(M*M+m2*m2-m13sq)/(2.0*M);
00324 double E3=(M*M+m3*m3-m12sq)/(2.0*M);
00325 double E1=M-E2-E3;
00326 double p1mom=sqrt(E1*E1-m1*m1);
00327 double p3mom=sqrt(E3*E3-m3*m3);
00328 double cost13=(2.0*E1*E3+m1*m1+m3*m3-m13sq)/(2.0*p1mom*p3mom);
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340 p4[2].set(E3,0.0,0.0,p3mom);
00341 p4[0].set(E1,p1mom*sqrt(1.0-cost13*cost13),0.0,p1mom*cost13);
00342 p4[1].set(E2,-p1mom*sqrt(1.0-cost13*cost13),0.0,-p1mom*cost13-p3mom);
00343
00344
00345
00346 double alpha = EvtRandom::Flat( EvtConst::twoPi );
00347 double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
00348 double gamma = EvtRandom::Flat( EvtConst::twoPi );
00349
00350 p4[0].applyRotateEuler( alpha, beta, gamma );
00351 p4[1].applyRotateEuler( alpha, beta, gamma );
00352 p4[2].applyRotateEuler( alpha, beta, gamma );
00353
00354 return 1.0+a/(m12sq*m12sq);
00355
00356 }
00357