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: EvtBreitWignerPdf.cpp,v 1.3 2009-03-16 15:55:55 robbep 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           0 :   : EvtIntegPdf1D(min,max), _m0(m0), _g0(g0)
      23           0 : {}
      24             : 
      25             : 
      26             : EvtBreitWignerPdf::EvtBreitWignerPdf(const EvtBreitWignerPdf& other)
      27           0 :   : EvtIntegPdf1D(other), _m0(other._m0), _g0(other._g0)
      28           0 : {}
      29             : 
      30             : 
      31             : EvtBreitWignerPdf::~EvtBreitWignerPdf()
      32           0 : {}
      33             : 
      34             : 
      35             : double EvtBreitWignerPdf::pdf(const EvtPoint1D& x) const
      36             : {
      37           0 :   double m = x.value();
      38           0 :   if((0 == (m - _m0)) && (0. == _g0)) {
      39             : 
      40           0 :     printf("Delta function Breit-Wigner\n");
      41           0 :     assert(0);
      42             :   }
      43             :   
      44           0 :   double ret = _g0/EvtConst::twoPi/((m-_m0)*(m-_m0)+_g0*_g0/4);
      45             :   
      46           0 :   return ret;
      47             : }
      48             : 
      49             : 
      50             : double EvtBreitWignerPdf::pdfIntegral(double m) const
      51             : {
      52             :   double itg = 0;
      53           0 :   if(_g0 == 0) {
      54             : 
      55           0 :     if(m > _m0) itg = 1.;
      56             :     else
      57           0 :       if(m < _m0) itg = 0.;
      58             :       else
      59             :         itg = 0.5;
      60             :   }
      61           0 :   else itg = atan((m-_m0)/(_g0/2.))/EvtConst::pi + 0.5; 
      62             : 
      63           0 :   return itg;
      64             : }
      65             : 
      66             : 
      67             : double EvtBreitWignerPdf::pdfIntegralInverse(double x) const
      68             : {
      69           0 :   if(x < 0 || x > 1) {
      70             :     
      71           0 :     printf("Invalid integral value %f\n",x);
      72           0 :     assert(0);
      73             :   }
      74             : 
      75           0 :   double m = _m0;
      76           0 :   if(_g0 != 0) m = _m0 + (_g0/2.)*tan(EvtConst::pi*(x-0.5));
      77             :  
      78           0 :   return m;
      79             : }
      80             : 
      81             : 
      82             : 
      83             :