Logo ROOT   6.10/09
Reference Guide
RooCBShape.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitModels *
4  * @(#)root/roofit:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 /** \class RooCBShape
18  \ingroup Roofit
19 
20 P.d.f implementing the Crystal Ball line shape
21 **/
22 
23 #include "RooFit.h"
24 
25 #include "Riostream.h"
26 #include "Riostream.h"
27 #include <math.h>
28 
29 #include "RooCBShape.h"
30 #include "RooAbsReal.h"
31 #include "RooRealVar.h"
32 #include "RooMath.h"
33 #include "TMath.h"
34 
35 #include "TError.h"
36 
37 using namespace std;
38 
40 
41 ////////////////////////////////////////////////////////////////////////////////
42 
43 Double_t RooCBShape::ApproxErf(Double_t arg) const
44 {
45  static const double erflim = 5.0;
46  if( arg > erflim )
47  return 1.0;
48  if( arg < -erflim )
49  return -1.0;
50 
51  return RooMath::erf(arg);
52 }
53 
54 ////////////////////////////////////////////////////////////////////////////////
55 
56 RooCBShape::RooCBShape(const char *name, const char *title,
57  RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma,
58  RooAbsReal& _alpha, RooAbsReal& _n) :
59  RooAbsPdf(name, title),
60  m("m", "Dependent", this, _m),
61  m0("m0", "M0", this, _m0),
62  sigma("sigma", "Sigma", this, _sigma),
63  alpha("alpha", "Alpha", this, _alpha),
64  n("n", "Order", this, _n)
65 {
66 }
67 
68 ////////////////////////////////////////////////////////////////////////////////
69 
70 RooCBShape::RooCBShape(const RooCBShape& other, const char* name) :
71  RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0),
72  sigma("sigma", this, other.sigma), alpha("alpha", this, other.alpha),
73  n("n", this, other.n)
74 {
75 }
76 
77 ////////////////////////////////////////////////////////////////////////////////
78 
80  Double_t t = (m-m0)/sigma;
81  if (alpha < 0) t = -t;
82 
83  Double_t absAlpha = fabs((Double_t)alpha);
84 
85  if (t >= -absAlpha) {
86  return exp(-0.5*t*t);
87  }
88  else {
89  Double_t a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
90  Double_t b= n/absAlpha - absAlpha;
91 
92  return a/TMath::Power(b - t, n);
93  }
94 }
95 
96 ////////////////////////////////////////////////////////////////////////////////
97 
98 Int_t RooCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const
99 {
100  if( matchArgs(allVars,analVars,m) )
101  return 1 ;
102 
103  return 0;
104 }
105 
106 ////////////////////////////////////////////////////////////////////////////////
107 
108 Double_t RooCBShape::analyticalIntegral(Int_t code, const char* rangeName) const
109 {
110  static const double sqrtPiOver2 = 1.2533141373;
111  static const double sqrt2 = 1.4142135624;
112 
113  R__ASSERT(code==1);
114  double result = 0.0;
115  bool useLog = false;
116 
117  if( fabs(n-1.0) < 1.0e-05 )
118  useLog = true;
119 
120  double sig = fabs((Double_t)sigma);
121 
122  double tmin = (m.min(rangeName)-m0)/sig;
123  double tmax = (m.max(rangeName)-m0)/sig;
124 
125  if(alpha < 0) {
126  double tmp = tmin;
127  tmin = -tmax;
128  tmax = -tmp;
129  }
130 
131  double absAlpha = fabs((Double_t)alpha);
132 
133  if( tmin >= -absAlpha ) {
134  result += sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
135  - ApproxErf(tmin/sqrt2) );
136  }
137  else if( tmax <= -absAlpha ) {
138  double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
139  double b = n/absAlpha - absAlpha;
140 
141  if(useLog) {
142  result += a*sig*( log(b-tmin) - log(b-tmax) );
143  }
144  else {
145  result += a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
146  - 1.0/(TMath::Power(b-tmax,n-1.0)) );
147  }
148  }
149  else {
150  double a = TMath::Power(n/absAlpha,n)*exp(-0.5*absAlpha*absAlpha);
151  double b = n/absAlpha - absAlpha;
152 
153  double term1 = 0.0;
154  if(useLog) {
155  term1 = a*sig*( log(b-tmin) - log(n/absAlpha));
156  }
157  else {
158  term1 = a*sig/(1.0-n)*( 1.0/(TMath::Power(b-tmin,n-1.0))
159  - 1.0/(TMath::Power(n/absAlpha,n-1.0)) );
160  }
161 
162  double term2 = sig*sqrtPiOver2*( ApproxErf(tmax/sqrt2)
163  - ApproxErf(-absAlpha/sqrt2) );
164 
165 
166  result += term1 + term2;
167  }
168 
169  return result;
170 }
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 /// Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
174 
176 {
177  RooArgSet dummy ;
178 
179  if (matchArgs(vars,dummy,m)) {
180  return 1 ;
181  }
182  return 0 ;
183 }
184 
185 ////////////////////////////////////////////////////////////////////////////////
186 
188 {
189  R__ASSERT(code==1) ;
190 
191  // The maximum value for given (m0,alpha,n,sigma)
192  return 1.0 ;
193 }
RooRealProxy m
Definition: RooCBShape.h:47
Double_t evaluate() const
Definition: RooCBShape.cxx:79
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooRealProxy m0
Definition: RooCBShape.h:48
#define R__ASSERT(e)
Definition: TError.h:96
virtual Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Definition: RooCBShape.cxx:108
int Int_t
Definition: RtypesCore.h:41
P.d.f implementing the Crystal Ball line shape.
Definition: RooCBShape.h:24
TArc * a
Definition: textangle.C:12
STL namespace.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:628
virtual Int_t getMaxVal(const RooArgSet &vars) const
Advertise that we know the maximum of self for given (m0,alpha,n,sigma)
Definition: RooCBShape.cxx:175
const Double_t sigma
virtual Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
Definition: RooCBShape.cxx:98
Double_t ApproxErf(Double_t arg) const
Definition: RooCBShape.cxx:43
virtual Double_t maxVal(Int_t code) const
Return maximum value for set of observables identified by code assigned in getMaxVal.
Definition: RooCBShape.cxx:187
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
TMarker * m
Definition: textangle.C:8
#define ClassImp(name)
Definition: Rtypes.h:336
Double_t min(const char *rname=0) const
Definition: RooRealProxy.h:56
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
RooRealProxy sigma
Definition: RooCBShape.h:49
static RooMathCoreReg dummy
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
Double_t max(const char *rname=0) const
Definition: RooRealProxy.h:57
RooRealProxy alpha
Definition: RooCBShape.h:50
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
RooRealProxy n
Definition: RooCBShape.h:51
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
static std::complex< double > erf(const std::complex< double > z)
complex erf function
Definition: RooMath.cxx:580
double result[121]
double exp(double)
const Int_t n
Definition: legend1.C:16
double log(double)