00001 #include "EvtGenBase/EvtPatches.hh" 00002 /******************************************************************************* 00003 * Project: BaBar detector at the SLAC PEP-II B-factory 00004 * Package: EvtGenBase 00005 * File: $Id: EvtBreitWignerPdf.cc,v 1.1.1.2 2007/10/26 05:03:14 pingrg Exp $ 00006 * Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002 00007 * 00008 * Copyright (C) 2002 Caltech 00009 *******************************************************************************/ 00010 00011 // Breit-Wigner shape PDF. If the width is zero it degenerates into a delta 00012 // function. The integral and its inverse can be still evaluated. 00013 00014 #include <assert.h> 00015 #include <stdio.h> 00016 #include <math.h> 00017 #include "EvtGenBase/EvtPatches.hh" 00018 #include "EvtGenBase/EvtBreitWignerPdf.hh" 00019 #include "EvtGenBase/EvtConst.hh" 00020 00021 EvtBreitWignerPdf::EvtBreitWignerPdf(double min, double max, double m0, double g0) 00022 : EvtIntegPdf1D(min,max), _m0(m0), _g0(g0) 00023 {} 00024 00025 00026 EvtBreitWignerPdf::EvtBreitWignerPdf(const EvtBreitWignerPdf& other) 00027 : EvtIntegPdf1D(other), _m0(other._m0), _g0(other._g0) 00028 {} 00029 00030 00031 EvtBreitWignerPdf::~EvtBreitWignerPdf() 00032 {} 00033 00034 00035 double EvtBreitWignerPdf::pdf(const EvtPoint1D& x) const 00036 { 00037 double m = x.value(); 00038 if((0 == (m - _m0)) && (0. == _g0)) { 00039 00040 printf("Delta function Breit-Wigner\n"); 00041 assert(0); 00042 } 00043 00044 double ret = _g0/EvtConst::twoPi/((m-_m0)*(m-_m0)+_g0*_g0/4); 00045 00046 return ret; 00047 } 00048 00049 00050 double EvtBreitWignerPdf::pdfIntegral(double m) const 00051 { 00052 double itg = 0; 00053 if(_g0 == 0) { 00054 00055 if(m > _m0) itg = 1.; 00056 else 00057 if(m < _m0) itg = 0.; 00058 else 00059 itg = 0.5; 00060 } 00061 else itg = atan((m-_m0)/(_g0/2.))/EvtConst::pi + 0.5; 00062 00063 return itg; 00064 } 00065 00066 00067 double EvtBreitWignerPdf::pdfIntegralInverse(double x) const 00068 { 00069 if(x < 0 || x > 1) { 00070 00071 printf("Invalid integral value %f\n",x); 00072 assert(0); 00073 } 00074 00075 double m = _m0; 00076 if(_g0 != 0) m = _m0 + (_g0/2.)*tan(EvtConst::pi*(x-0.5)); 00077 00078 return m; 00079 } 00080 00081 00082 00083