TEvtGen/EvtGenBase/EvtBreitWignerPdf.cxx

   1 #include "EvtGenBase/EvtPatches.hh"
   2 /*******************************************************************************
   3  * Project: BaBar detector at the SLAC PEP-II B-factory
   4  * Package: EvtGenBase
   5  *    File: $Id: EvtBreitWignerPdf.cc,v 1.14 2004/12/21 19:58:41 ryd Exp $
   6  *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
   7  *
   8  * Copyright (C) 2002 Caltech
   9  *******************************************************************************/
  10
  11 // Breit-Wigner shape PDF. If the width is zero it degenerates into a delta
  12 // function. The integral and its inverse can be still evaluated.
  13
  14 #include <assert.h>
  15 #include <stdio.h>
  16 #include <math.h>
  17 #include "EvtGenBase/EvtPatches.hh"
  18 #include "EvtGenBase/EvtBreitWignerPdf.hh"
  19 #include "EvtGenBase/EvtConst.hh"
  20
  21 EvtBreitWignerPdf::EvtBreitWignerPdf(double min, double max, double m0, double g0)
  22   : EvtIntegPdf1D(min,max), _m0(m0), _g0(g0)
  23 {}
  24
  25
  26 EvtBreitWignerPdf::EvtBreitWignerPdf(const EvtBreitWignerPdf& other)
  27   : EvtIntegPdf1D(other), _m0(other._m0), _g0(other._g0)
  28 {}
  29
  30
  31 EvtBreitWignerPdf::~EvtBreitWignerPdf()
  32 {}
  33
  34
  35 double EvtBreitWignerPdf::pdf(const EvtPoint1D& x) const
  36 {
  37   double m = x.value();
  38   if((0 == (m - _m0)) && (0. == _g0)) {
  39
  40     printf("Delta function Breit-Wigner\n");
  41     assert(0);
  42   }
  43
  44   double ret = _g0/EvtConst::twoPi/((m-_m0)*(m-_m0)+_g0*_g0/4);
  45
  46   return ret;
  47 }
  48
  49
  50 double EvtBreitWignerPdf::pdfIntegral(double m) const
  51 {
  52   double itg = 0;
  53   if(_g0 == 0) {
  54
  55     if(m > _m0) itg = 1.;
  56     else
  57       if(m < _m0) itg = 0.;
  58       else
  59         itg = 0.5;
  60   }
  61   else itg = atan((m-_m0)/(_g0/2.))/EvtConst::pi + 0.5;
  62
  63   return itg;
  64 }
  65
  66
  67 double EvtBreitWignerPdf::pdfIntegralInverse(double x) const
  68 {
  69   if(x < 0 || x > 1) {
  70
  71     printf("Invalid integral value %f\n",x);
  72     assert(0);
  73   }
  74
  75   double m = _m0;
  76   if(_g0 != 0) m = _m0 + (_g0/2.)*tan(EvtConst::pi*(x-0.5));
  77
  78   return m;
  79 }
  80
  81
  82
  83