TEvtGen/EvtGen/EvtGenBase/EvtDalitzResPdf.cpp

   1 #include "EvtGenBase/EvtPatches.hh"
   2 /*******************************************************************************
   3  * Project: BaBar detector at the SLAC PEP-II B-factory
   4  * Package: EvtGenBase
   5  *    File: $Id: EvtDalitzResPdf.cpp,v 1.3 2009-03-16 15:54:07 robbep Exp $
   6  *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
   7  *
   8  * Copyright (C) 2002 Caltech
   9  *******************************************************************************/
  10
  11 #include <stdio.h>
  12 #include <math.h>
  13 #include "EvtGenBase/EvtPatches.hh"
  14 #include "EvtGenBase/EvtConst.hh"
  15 #include "EvtGenBase/EvtDalitzResPdf.hh"
  16 #include "EvtGenBase/EvtDalitzCoord.hh"
  17 #include "EvtGenBase/EvtRandom.hh"
  18 using namespace EvtCyclic3;
  19
  20 EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzPlot& dp,
  21                                  double _m0, double _g0, EvtCyclic3::Pair pair)
  22   : EvtPdf<EvtDalitzPoint>(),
  23   _dp(dp), _m0(_m0), _g0(_g0), _pair(pair)
  24 {}
  25
  26
  27 EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzResPdf& other)
  28   : EvtPdf<EvtDalitzPoint>(other),
  29   _dp(other._dp),_m0(other._m0), _g0(other._g0), _pair(other._pair)
  30 {}
  31
  32
  33 EvtDalitzResPdf::~EvtDalitzResPdf()
  34 {}
  35
  36 EvtValError EvtDalitzResPdf::compute_integral(int N) const
  37 {
  38   assert(N != 0);
  39
  40   EvtCyclic3::Pair i = _pair;
  41   EvtCyclic3::Pair j = EvtCyclic3::next(i);
  42
  43   // Trapezoidal integral
  44
  45   double dh = (_dp.qAbsMax(j) - _dp.qAbsMin(j))/((double) N);
  46   double sum = 0;
  47
  48   int ii;
  49   for(ii=1;ii<N;ii++) {
  50
  51     double x = _dp.qAbsMin(j) + ii*dh;
  52     double min = (_dp.qMin(i,j,x) - _m0*_m0)/_m0/_g0;
  53     double max = (_dp.qMax(i,j,x) - _m0*_m0)/_m0/_g0;
  54     double itg = 1/EvtConst::pi*(atan(max) - atan(min));
  55     sum += itg;
  56   }
  57   EvtValError ret(sum*dh,0.);
  58
  59   return ret;
  60 }
  61
  62
  63 EvtDalitzPoint EvtDalitzResPdf::randomPoint()
  64 {
  65   // Random point generation must be done in a box encompassing the
  66   // Dalitz plot
  67
  68
  69   EvtCyclic3::Pair i = _pair;
  70   EvtCyclic3::Pair j = EvtCyclic3::next(i);
  71   double min = 1/EvtConst::pi*atan((_dp.qAbsMin(i) - _m0*_m0)/_m0/_g0);
  72   double max = 1/EvtConst::pi*atan((_dp.qAbsMax(i) - _m0*_m0)/_m0/_g0);
  73
  74   int n = 0;
  75   while(n++ < 1000) {
  76
  77     double qj = EvtRandom::Flat(_dp.qAbsMin(j),_dp.qAbsMax(j));
  78     double r = EvtRandom::Flat(min,max);
  79     double qi = tan(EvtConst::pi*r)*_g0*_m0 + _m0*_m0;
  80     EvtDalitzCoord x(i,qi,j,qj);
  81     EvtDalitzPoint ret(_dp,x);
  82     if(ret.isValid()) return ret;
  83   }
  84
  85   // All generated points turned out to be outside of the Dalitz plot
  86   // (in the outer box)
  87
  88   printf("No point generated for dalitz plot after 1000 tries\n");
  89   return EvtDalitzPoint(0.,0.,0.,0.,0.,0.);
  90 }
  91
  92
  93 double EvtDalitzResPdf::pdf(const EvtDalitzPoint& x) const
  94 {
  95   EvtCyclic3::Pair i = _pair;
  96   double dq = x.q(i) - _m0*_m0;
  97   return 1/EvtConst::pi*_g0*_m0/(dq*dq + _g0*_g0*_m0*_m0);
  98 }
  99
 100
 101 double EvtDalitzResPdf::pdfMaxValue() const
 102 {
 103   return 1/(EvtConst::pi*_g0*_m0);
 104 }
 105