Line data    Source code

       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           0 :   : EvtPdf<EvtDalitzPoint>(), 
      23           0 :   _dp(dp), _m0(_m0), _g0(_g0), _pair(pair)
      24           0 : {}
      25             : 
      26             : 
      27             : EvtDalitzResPdf::EvtDalitzResPdf(const EvtDalitzResPdf& other)
      28           0 :   : EvtPdf<EvtDalitzPoint>(other), 
      29           0 :   _dp(other._dp),_m0(other._m0), _g0(other._g0), _pair(other._pair)
      30           0 : {}
      31             : 
      32             : 
      33             : EvtDalitzResPdf::~EvtDalitzResPdf()
      34           0 : {}
      35             : 
      36             : EvtValError EvtDalitzResPdf::compute_integral(int N) const
      37             : {
      38           0 :   assert(N != 0);
      39             :   
      40           0 :   EvtCyclic3::Pair i = _pair;
      41           0 :   EvtCyclic3::Pair j = EvtCyclic3::next(i);
      42             : 
      43             :   // Trapezoidal integral
      44             : 
      45           0 :   double dh = (_dp.qAbsMax(j) - _dp.qAbsMin(j))/((double) N);
      46             :   double sum = 0;
      47             :   
      48             :   int ii;
      49           0 :   for(ii=1;ii<N;ii++) {
      50             :     
      51           0 :     double x = _dp.qAbsMin(j) + ii*dh;
      52           0 :     double min = (_dp.qMin(i,j,x) - _m0*_m0)/_m0/_g0;
      53           0 :     double max = (_dp.qMax(i,j,x) - _m0*_m0)/_m0/_g0;
      54           0 :     double itg = 1/EvtConst::pi*(atan(max) - atan(min));
      55           0 :     sum += itg;
      56             :   }
      57           0 :   EvtValError ret(sum*dh,0.); 
      58             :   
      59             :   return ret;
      60           0 : }
      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           0 :   EvtCyclic3::Pair i = _pair;
      70           0 :   EvtCyclic3::Pair j = EvtCyclic3::next(i);  
      71           0 :   double min = 1/EvtConst::pi*atan((_dp.qAbsMin(i) - _m0*_m0)/_m0/_g0);
      72           0 :   double max = 1/EvtConst::pi*atan((_dp.qAbsMax(i) - _m0*_m0)/_m0/_g0);
      73             : 
      74             :   int n = 0;
      75           0 :   while(n++ < 1000) {
      76             : 
      77           0 :     double qj = EvtRandom::Flat(_dp.qAbsMin(j),_dp.qAbsMax(j));
      78           0 :     double r = EvtRandom::Flat(min,max);
      79           0 :     double qi = tan(EvtConst::pi*r)*_g0*_m0 + _m0*_m0;
      80           0 :     EvtDalitzCoord x(i,qi,j,qj);
      81           0 :     EvtDalitzPoint ret(_dp,x);
      82           0 :     if(ret.isValid()) return ret;
      83           0 :   }
      84             :   
      85             :   // All generated points turned out to be outside of the Dalitz plot
      86             :   // (in the outer box)
      87             :   
      88           0 :   printf("No point generated for dalitz plot after 1000 tries\n");
      89           0 :   return EvtDalitzPoint(0.,0.,0.,0.,0.,0.);
      90           0 : }
      91             : 
      92             : 
      93             : double EvtDalitzResPdf::pdf(const EvtDalitzPoint& x) const
      94             : {
      95           0 :   EvtCyclic3::Pair i = _pair;
      96           0 :   double dq = x.q(i) - _m0*_m0;
      97           0 :   return 1/EvtConst::pi*_g0*_m0/(dq*dq + _g0*_g0*_m0*_m0);
      98             : }
      99             : 
     100             : 
     101             : double EvtDalitzResPdf::pdfMaxValue() const
     102             : {
     103           0 :   return 1/(EvtConst::pi*_g0*_m0);
     104             : }
     105             :