00001 #include "EvtGenBase/EvtPatches.hh"
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include <stdio.h>
00012 #include <math.h>
00013 #include "EvtGenBase/EvtPatches.hh"
00014 #include "EvtGenBase/EvtConst.hh"
00015 #include "EvtGenBase/EvtDalitzResPdf.hh"
00016 #include "EvtGenBase/EvtDalitzCoord.hh"
00017 #include "EvtGenBase/EvtRandom.hh"
00018 using namespace EvtCyclic3;
00019
00020 EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzPlot& dp,
00021 double _m0, double _g0, EvtCyclic3::Pair pair)
00022 : EvtPdf<EvtDalitzPoint>(),
00023 _dp(dp), _m0(_m0), _g0(_g0), _pair(pair)
00024 {}
00025
00026
00027 EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzResPdf& other)
00028 : EvtPdf<EvtDalitzPoint>(other),
00029 _dp(other._dp),_m0(other._m0), _g0(other._g0), _pair(other._pair)
00030 {}
00031
00032
00033 EvtDalitzResPdf::~EvtDalitzResPdf()
00034 {}
00035
00036 EvtValError EvtDalitzResPdf::compute_integral(int N) const
00037 {
00038 assert(N != 0);
00039
00040 EvtCyclic3::Pair i = _pair;
00041 EvtCyclic3::Pair j = EvtCyclic3::next(i);
00042
00043
00044
00045 double dh = (_dp.qAbsMax(j) - _dp.qAbsMin(j))/((double) N);
00046 double sum = 0;
00047
00048 int ii;
00049 for(ii=1;ii<N;ii++) {
00050
00051 double x = _dp.qAbsMin(j) + ii*dh;
00052 double min = (_dp.qMin(i,j,x) - _m0*_m0)/_m0/_g0;
00053 double max = (_dp.qMax(i,j,x) - _m0*_m0)/_m0/_g0;
00054 double itg = 1/EvtConst::pi*(atan(max) - atan(min));
00055 sum += itg;
00056 }
00057 EvtValError ret(sum*dh,0.);
00058
00059 return ret;
00060 }
00061
00062
00063 EvtDalitzPoint EvtDalitzResPdf::randomPoint()
00064 {
00065
00066
00067
00068
00069 EvtCyclic3::Pair i = _pair;
00070 EvtCyclic3::Pair j = EvtCyclic3::next(i);
00071 double min = 1/EvtConst::pi*atan((_dp.qAbsMin(i) - _m0*_m0)/_m0/_g0);
00072 double max = 1/EvtConst::pi*atan((_dp.qAbsMax(i) - _m0*_m0)/_m0/_g0);
00073
00074 int n = 0;
00075 while(n++ < 1000) {
00076
00077 double qj = EvtRandom::Flat(_dp.qAbsMin(j),_dp.qAbsMax(j));
00078 double r = EvtRandom::Flat(min,max);
00079 double qi = tan(EvtConst::pi*r)*_g0*_m0 + _m0*_m0;
00080 EvtDalitzCoord x(i,qi,j,qj);
00081 EvtDalitzPoint ret(_dp,x);
00082 if(ret.isValid()) return ret;
00083 }
00084
00085
00086
00087
00088 printf("No point generated for dalitz plot after 1000 tries\n");
00089 assert(0);
00090 }
00091
00092
00093 double EvtDalitzResPdf::pdf(const EvtDalitzPoint& x) const
00094 {
00095 EvtCyclic3::Pair i = _pair;
00096 double dq = x.q(i) - _m0*_m0;
00097 return 1/EvtConst::pi*_g0*_m0/(dq*dq + _g0*_g0*_m0*_m0);
00098 }
00099
00100
00101 double EvtDalitzResPdf::pdfMaxValue() const
00102 {
00103 return 1/(EvtConst::pi*_g0*_m0);
00104 }
00105