]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | #include "EvtGenBase/EvtPatches.hh" |
2 | /******************************************************************************* | |
3 | * Project: BaBar detector at the SLAC PEP-II B-factory | |
4 | * Package: EvtGenBase | |
0ca57c2f | 5 | * File: $Id: EvtBreitWignerPdf.cpp,v 1.3 2009-03-16 15:55:55 robbep Exp $ |
da0e9ce3 | 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 |