BOSS 7.0.8
BESIII Offline Software System
Loading...
Searching...
No Matches
NumRecipes.cxx
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// File and Version Information:
3// $Id: NumRecipes.cxx,v 1.1.1.1 2005/04/21 01:17:17 zhangy Exp $
4//
5// Description:
6//
7// Environment:
8// Software developed for the BaBar Detector at the SLAC B-Factory.
9//
10// Author List:
11// Bob Jacobsen, Ed Iskander
12//
13// Copyright Information:
14// Copyright (C) 1996 Lawrence Berkeley Laboratory
15//
16//------------------------------------------------------------------------
17//#include "BaBar/BaBar.hh"
18
19//-----------------------
20// This Class's Header --
21//-----------------------
23
24//-------------
25// C Headers --
26//-------------
27extern "C" {
28#include <assert.h>
29#include <math.h>
30#include <stdlib.h>
31}
32
33//---------------
34// C++ Headers --
35//---------------
36#include <iostream>
37//#include "ErrLogger/ErrLog.hh"
38using std::endl;
39
40double NumRecipes::gammln(double xx)
41{
42 double x,y,tmp,ser;
43 static double cof[6]={76.18009172947146,-86.50532032941677,
44 24.01409824083091,-1.231739572450155,
45 0.1208650973866179e-2,-0.5395239384953e-5};
46 int j;
47
48 y=x=xx;
49 tmp=x+5.5;
50 tmp -= (x+0.5)*log(tmp);
51 ser=1.000000000190015;
52 for (j=0;j<=5;j++) ser += cof[j]/++y;
53 return -tmp+log(2.5066282746310005*ser/x);
54}
55
56double
57NumRecipes::gammp(double a, double x)
58{
59 double gamser,gammcf,gln;
60
61 if (x < 0.0 || a <= 0.0) std::cout<<"ErrMsg(error)" <<" Invalid arguments in routine gammp x=" << x << " a=" << a << endl;
62 if (x < (a+1.0)) {
63 gser(&gamser,a,x,&gln);
64 return gamser;
65 } else {
66 gcf(&gammcf,a,x,&gln);
67 return 1.0-gammcf;
68 }
69}
70
71double
72NumRecipes::gammq(double a, double x)
73{
74 double gamser,gammcf,gln;
75
76 if (x < 0.0 || a <= 0.0) recipesErr(" Invalid arguments in routine GAMMQ");
77 if (x < (a+1.0)) {
78 gser(&gamser,a,x,&gln);
79 return 1.0-gamser;
80 } else {
81 gcf(&gammcf,a,x,&gln);
82 return gammcf;
83 }
84}
85
86void NumRecipes::gcf(double* gammcf, double a, double x, double* gln)
87 {
88 int n;
89 double gold=0.0,g,fac=1.0,b1=1.0;
90 double b0=0.0,anf,ana,an,a1,a0=1.0;
91
92 *gln=gammln(a);
93 a1=x;
94 for (n=1;n<=NUMREC_ITMAX;n++) {
95 an=(double) n;
96 ana=an-a;
97 a0=(a1+a0*ana)*fac;
98 b0=(b1+b0*ana)*fac;
99 anf=an*fac;
100 a1=x*a0+anf*a1;
101 b1=x*b0+anf*b1;
102 if (a1) {
103 fac=1.0/a1;
104 g=b1*fac;
105 if (fabs((g-gold)/g) < NUMREC_EPS) {
106 *gammcf=exp(-x+a*log(x)-(*gln))*g;
107 return;
108 }
109 gold=g;
110 }
111 }
112 recipesErr(" a too large, NUMREC_ITMAX too small in routine GCF");
113 }
114
115void NumRecipes::gser(double* gamser, double a, double x, double* gln)
116{
117 int n;
118 double sum,del,ap;
119
120 *gln=gammln(a);
121 if (x <= 0.0) {
122 if (x < 0.0) recipesErr(" x less than 0 in routine GSER");
123 *gamser=0.0;
124 return;
125 } else {
126 ap=a;
127 del=sum=1.0/a;
128 for (n=1;n<=NUMREC_ITMAX;n++) {
129 ap += 1.0;
130 del *= x/ap;
131 sum += del;
132 if (fabs(del) < fabs(sum)*NUMREC_EPS) {
133 *gamser=sum*exp(-x+a*log(x)-(*gln));
134 return;
135 }
136 }
137 recipesErr(" a too large, NUMREC_ITMAX too small in routine GSER");
138 return;
139 }
140}
141
142void NumRecipes::recipesErr(const char* c)
143{
144 std::cout<<"ErrMsg(fatal)" << " Numerical Recipes run-time error...\n" << c
145 << "\n ...now exiting to system..." << std::endl;
146 ::abort();
147}
148
const Int_t n
Double_t x[10]
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
static double gammq(double a, double x)
Definition: NumRecipes.cxx:72
static void gser(double *gamser, double a, double x, double *gln)
Definition: NumRecipes.cxx:115
static double gammln(double x)
Definition: NumRecipes.cxx:40
static void gcf(double *gammcf, double a, double x, double *gln)
Definition: NumRecipes.cxx:86
static double gammp(double a, double x)
Definition: NumRecipes.cxx:57
double y[1000]