4 {
5
6
7 static float srtopi = 2.0 / sqrt( 2.0 *
M_PI );
8 static float upl = 100.0;
9
10 float prob = 0.0;
11 if ( ndof <= 0 ) return prob;
12 if ( chisq < 0.0 ) return prob;
13 if ( ndof <= 60 )
14 {
15
16 if ( chisq > upl ) return prob;
17 float sum =
exp( -0.5 * chisq );
18 float term = sum;
19
20 int m = ndof / 2;
21 if ( 2 * m == ndof )
22 {
23 if ( m == 1 ) return sum;
24 for ( int i = 2; i <= m; i++ )
25 {
26 term = 0.5 * term * chisq / ( i - 1 );
27 sum += term;
28 }
29 return sum;
30
31 }
32 else
33 {
34
35 float srty = sqrt( chisq );
36 float temp = srty / M_SQRT2;
37 prob = erfc( temp );
38 if ( ndof == 1 ) return prob;
39 if ( ndof == 3 ) return ( srtopi * srty * sum + prob );
40 m--;
41 for ( int i = 1; i <= m; i++ )
42 {
43 term = term * chisq / ( 2 * i + 1 );
44 sum += term;
45 }
46 return ( srtopi * srty * sum + prob );
47
48 }
49 }
50 else
51 {
52
53 float srty = sqrt( chisq ) - sqrt( ndof - 0.5 );
54 if ( srty < 12.0 ) prob = 0.5 * erfc( srty );
55 return prob;
56
57 }
58
59
60
61}
EvtComplex exp(const EvtComplex &c)