NCFS-Pack
A generic (astro)particle physics analysis framework
 
Loading...
Searching...
No Matches
NcMath.cxx
Go to the documentation of this file.
1
31
33
54
55#include "NcMath.h"
56#include "Riostream.h"
57
58ClassImp(NcMath); // Class implementation to enable ROOT I/O
59
61NcMath::NcMath() : TObject()
62{
68}
69
78
79NcMath::NcMath(const NcMath& m) : TObject(m)
80{
86}
87
88Double_t NcMath::Zeta(Double_t x, Int_t nterms) const
89{
103
104 if (x<=1)
105 {
106 cout << "*Zeta(x)* Wrong argument x = " << x << endl;
107 return 0;
108 }
109
110 Double_t zeta=0;
111 Double_t r=0;
112 for (Int_t i=1; i<=nterms; i++)
113 {
114 r=double(i);
115 zeta+=1./pow(r,x);
116 }
117
118 return zeta;
119}
120
121Double_t NcMath::Gamma(Double_t z) const
122{
135
136 if (z<=0.)
137 {
138 cout << "*Gamma(z)* Wrong argument z = " << z << endl;
139 return 0;
140 }
141
142 Double_t v=LnGamma(z);
143 return exp(v);
144}
145
146Double_t NcMath::Gamma(Double_t a,Double_t x,Int_t mode) const
147{
165
166 Double_t value=0;
167
168 if (a<=0.)
169 {
170 cout << "*Gamma(a,x)* Invalid argument a = " << a << endl;
171 return 0;
172 }
173
174 if (x<=0.)
175 {
176 if (x<0) cout << "*Gamma(a,x)* Invalid argument x = " << x << endl;
177 return 0;
178 }
179
180 if (x<(a+1.))
181 {
182 value=GamSer(a,x);
183 }
184 else
185 {
186 value=GamCf(a,x);
187 }
188
189 if (mode) value=value*Gamma(a);
190 return value;
191}
192
193Double_t NcMath::LnGamma(Double_t z) const
194{
209
210 if (z<=0.)
211 {
212 cout << "*LnGamma(z)* Wrong argument z = " << z << endl;
213 return 0;
214 }
215
216 // Coefficients for the series expansion
217 Double_t c[7];
218 c[0]= 2.5066282746310005;
219 c[1]= 76.18009172947146;
220 c[2]=-86.50532032941677;
221 c[3]= 24.01409824083091;
222 c[4]= -1.231739572450155;
223 c[5]= 0.1208650973866179e-2;
224 c[6]= -0.5395239384953e-5;
225
226 Double_t x=z;
227 Double_t y=x;
228 Double_t tmp=x+5.5;
229 tmp=(x+0.5)*log(tmp)-tmp;
230 Double_t ser=1.000000000190015;
231 for (Int_t i=1; i<7; i++)
232 {
233 y+=1.;
234 ser+=c[i]/y;
235 }
236 Double_t v=tmp+log(c[0]*ser/x);
237 return v;
238}
239
240Double_t NcMath::LnGamma(Double_t a,Double_t x,Int_t mode) const
241{
256
257 Double_t value=0;
258 Double_t gammaP=0;
259
260 gammaP=Gamma(a,x,0);
261
262 if (gammaP)
263 {
264 value=log(gammaP);
265 if (mode) value=value+LnGamma(a);
266 }
267
268 return value;
269}
270
271Double_t NcMath::GamSer(Double_t a,Double_t x) const
272{
284
285 Int_t itmax=100; // Maximum number of iterations
286 Double_t eps=3.e-7; // Relative accuracy
287
288 if (a<=0.)
289 {
290 cout << "*GamSer(a,x)* Invalid argument a = " << a << endl;
291 return 0;
292 }
293
294 if (x<=0.)
295 {
296 if (x<0) cout << "*GamSer(a,x)* Invalid argument x = " << x << endl;
297 return 0;
298 }
299
300 Double_t gln=LnGamma(a);
301 Double_t ap=a;
302 Double_t sum=1./a;
303 Double_t del=sum;
304 for (Int_t n=1; n<=itmax; n++)
305 {
306 ap+=1.;
307 del=del*x/ap;
308 sum+=del;
309 if (fabs(del)<fabs(sum*eps)) break;
310 if (n==itmax) cout << "*GamSer(a,x)* a too large or itmax too small" << endl;
311 }
312 Double_t v=sum*exp(-x+a*log(x)-gln);
313 return v;
314}
315
316Double_t NcMath::GamCf(Double_t a,Double_t x) const
317{
329
330 Int_t itmax=100; // Maximum number of iterations
331 Double_t eps=3.e-7; // Relative accuracy
332 Double_t fpmin=1.e-30; // Smallest Double_t value allowed here
333
334 if (a<=0.)
335 {
336 cout << "*GamCf(a,x)* Invalid argument a = " << a << endl;
337 return 0;
338 }
339
340 if (x<=0.)
341 {
342 if (x<0) cout << "*GamCf(a,x)* Invalid argument x = " << x << endl;
343 return 0;
344 }
345
346 Double_t gln=LnGamma(a);
347 Double_t b=x+1.-a;
348 Double_t c=1./fpmin;
349 Double_t d=1./b;
350 Double_t h=d;
351 Double_t an,del;
352 for (Int_t i=1; i<=itmax; i++)
353 {
354 an=double(-i)*(double(i)-a);
355 b+=2.;
356 d=an*d+b;
357 if (fabs(d)<fpmin) d=fpmin;
358 c=b+an/c;
359 if (fabs(c)<fpmin) c=fpmin;
360 d=1./d;
361 del=d*c;
362 h=h*del;
363 if (fabs(del-1.)<eps) break;
364 if (i==itmax) cout << "*GamCf(a,x)* a too large or itmax too small" << endl;
365 }
366 Double_t v=exp(-x+a*log(x)-gln)*h;
367 return (1.-v);
368}
369
370Double_t NcMath::Erf(Double_t x) const
371{
379
380 return (1.-Erfc(x));
381}
382
383Double_t NcMath::Erfc(Double_t x) const
384{
397
398 // The parameters of the Chebyshev fit
399 const Double_t ka1=-1.26551223, ka2=1.00002368,
400 ka3= 0.37409196, ka4=0.09678418,
401 ka5=-0.18628806, ka6=0.27886807,
402 ka7=-1.13520398, ka8=1.48851587,
403 ka9=-0.82215223, ka10=0.17087277;
404
405 Double_t v=1.; // The return value
406
407 Double_t z=fabs(x);
408
409 if (z <= 0.) return v; // erfc(0)=1
410
411 Double_t t=1./(1.+0.5*z);
412
413 v=t*exp((-z*z)
414 +ka1+t*(ka2+t*(ka3+t*(ka4+t*(ka5+t*(ka6+t*(ka7+t*(ka8+t*(ka9+t*ka10)))))))));
415
416 if (x < 0.) v=2.-v; // erfc(-x)=2-erfc(x)
417
418 return v;
419}
420
421Double_t NcMath::Prob(Double_t chi2,Int_t ndf,Int_t mode) const
422{
454
455 if (ndf <= 0) return 0; // Set CL to zero in case ndf<=0
456
457 if (chi2 <= 0.)
458 {
459 if (chi2 < 0.)
460 {
461 return 0;
462 }
463 else
464 {
465 return 1;
466 }
467 }
468
469 Double_t v=-1.;
470
471 switch (mode)
472 {
473 case 1: // Exact expression for ndf=1 as alternative for the gamma function
474 if (ndf==1) v=1.-Erf(sqrt(chi2)/sqrt(2.));
475 break;
476
477 case 2: // Gaussian approximation for large ndf (i.e. ndf>30) as alternative for the gamma function
478 if (ndf>30)
479 {
480 Double_t q=sqrt(2.*chi2)-sqrt(double(2*ndf-1));
481 if (q>0.) v=0.5*(1.-Erf(q/sqrt(2.)));
482 }
483 break;
484 }
485
486 if (v<0.)
487 {
488 // Evaluate the incomplete gamma function
489 Double_t a=double(ndf)/2.;
490 Double_t x=chi2/2.;
491 v=1.-Gamma(a,x);
492 }
493
494 return v;
495}
496
497Double_t NcMath::BesselI0(Double_t x) const
498{
512
513 // Parameters of the polynomial approximation
514 const Double_t kp1=1.0, kp2=3.5156229, kp3=3.0899424,
515 kp4=1.2067492, kp5=0.2659732, kp6=3.60768e-2, kp7=4.5813e-3;
516
517 const Double_t kq1= 0.39894228, kq2= 1.328592e-2, kq3= 2.25319e-3,
518 kq4=-1.57565e-3, kq5= 9.16281e-3, kq6=-2.057706e-2,
519 kq7= 2.635537e-2, kq8=-1.647633e-2, kq9= 3.92377e-3;
520
521 Double_t ax=fabs(x);
522
523 Double_t y=0,result=0;
524
525 if (ax < 3.75)
526 {
527 y=pow(x/3.75,2);
528 result=kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7)))));
529 }
530 else
531 {
532 y=3.75/ax;
533 result=(exp(ax)/sqrt(ax))
534 *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*(kq7+y*(kq8+y*kq9))))))));
535 }
536
537 return result;
538}
539
540Double_t NcMath::BesselK0(Double_t x) const
541{
555
556 // Parameters of the polynomial approximation
557 const Double_t kp1=-0.57721566, kp2=0.42278420, kp3=0.23069756,
558 kp4= 3.488590e-2, kp5=2.62698e-3, kp6=1.0750e-4, kp7=7.4e-6;
559
560 const Double_t kq1= 1.25331414, kq2=-7.832358e-2, kq3= 2.189568e-2,
561 kq4=-1.062446e-2, kq5= 5.87872e-3, kq6=-2.51540e-3, kq7=5.3208e-4;
562
563 if (x <= 0)
564 {
565 cout << " *BesselK0* Invalid argument x = " << x << endl;
566 return 0;
567 }
568
569 Double_t y=0,result=0;
570
571 if (x <= 2)
572 {
573 y=x*x/4.;
574 result=(-log(x/2.)*BesselI0(x))
575 +(kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7))))));
576 }
577 else
578 {
579 y=2./x;
580 result=(exp(-x)/sqrt(x))
581 *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*kq7))))));
582 }
583
584 return result;
585}
586
587Double_t NcMath::BesselI1(Double_t x) const
588{
602
603 // Parameters of the polynomial approximation
604 const Double_t kp1=0.5, kp2=0.87890594, kp3=0.51498869,
605 kp4=0.15084934, kp5=2.658733e-2, kp6=3.01532e-3, kp7=3.2411e-4;
606
607 const Double_t kq1= 0.39894228, kq2=-3.988024e-2, kq3=-3.62018e-3,
608 kq4= 1.63801e-3, kq5=-1.031555e-2, kq6= 2.282967e-2,
609 kq7=-2.895312e-2, kq8= 1.787654e-2, kq9=-4.20059e-3;
610
611 Double_t ax=fabs(x);
612
613 Double_t y=0,result=0;
614
615 if (ax < 3.75)
616 {
617 y=pow(x/3.75,2);
618 result=x*(kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7))))));
619 }
620 else
621 {
622 y=3.75/ax;
623 result=(exp(ax)/sqrt(ax))
624 *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*(kq7+y*(kq8+y*kq9))))))));
625 if (x < 0) result=-result;
626 }
627
628 return result;
629}
630
631Double_t NcMath::BesselK1(Double_t x) const
632{
646
647 // Parameters of the polynomial approximation
648 const Double_t kp1= 1., kp2= 0.15443144, kp3=-0.67278579,
649 kp4=-0.18156897, kp5=-1.919402e-2, kp6=-1.10404e-3, kp7=-4.686e-5;
650
651 const Double_t kq1= 1.25331414, kq2= 0.23498619, kq3=-3.655620e-2,
652 kq4= 1.504268e-2, kq5=-7.80353e-3, kq6= 3.25614e-3, kq7=-6.8245e-4;
653
654 if (x <= 0)
655 {
656 cout << " *BesselK1* Invalid argument x = " << x << endl;
657 return 0;
658 }
659
660 Double_t y=0,result=0;
661
662 if (x <= 2)
663 {
664 y=x*x/4.;
665 result=(log(x/2.)*BesselI1(x))
666 +(1./x)*(kp1+y*(kp2+y*(kp3+y*(kp4+y*(kp5+y*(kp6+y*kp7))))));
667 }
668 else
669 {
670 y=2./x;
671 result=(exp(-x)/sqrt(x))
672 *(kq1+y*(kq2+y*(kq3+y*(kq4+y*(kq5+y*(kq6+y*kq7))))));
673 }
674
675 return result;
676}
677
678Double_t NcMath::BesselK(Int_t n,Double_t x) const
679{
694
695 if (x <= 0 || n < 0)
696 {
697 cout << " *BesselK* Invalid argument(s) (n,x) = (" << n << " , " << x << ")" << endl;
698 return 0;
699 }
700
701 if (n==0) return BesselK0(x);
702
703 if (n==1) return BesselK1(x);
704
705 // Perform upward recurrence for all x
706 Double_t tox=2./x;
707 Double_t bkm=BesselK0(x);
708 Double_t bk=BesselK1(x);
709 Double_t bkp=0;
710 for (Int_t j=1; j<n; j++)
711 {
712 bkp=bkm+double(j)*tox*bk;
713 bkm=bk;
714 bk=bkp;
715 }
716
717 return bk;
718}
719
720Double_t NcMath::BesselI(Int_t n,Double_t x) const
721{
736
737 Int_t iacc=40; // Increase to enhance accuracy
738 Double_t bigno=1.e10, bigni=1.e-10;
739
740 if (n < 0)
741 {
742 cout << " *BesselI* Invalid argument (n,x) = (" << n << " , " << x << ")" << endl;
743 return 0;
744 }
745
746 if (n==0) return BesselI0(x);
747
748 if (n==1) return BesselI1(x);
749
750 if (fabs(x) < 1.e-10) return 0;
751
752 Double_t tox=2./fabs(x);
753 Double_t bip=0,bim=0;
754 Double_t bi=1;
755 Double_t result=0;
756 Int_t m=2*((n+int(sqrt(float(iacc*n))))); // Downward recurrence from even m
757 for (Int_t j=m; j<=1; j--)
758 {
759 bim=bip+double(j)*tox*bi;
760 bip=bi;
761 bi=bim;
762 if (fabs(bi) > bigno) // Renormalise to prevent overflows
763 {
764 result*=bigni;
765 bi*=bigni;
766 bip*=bigni;
767 }
768 if (j==n) result=bip;
769 }
770
771 result*=BesselI0(x)/bi; // Normalise with I0(x)
772 if ((x < 0) && (n%2 == 1)) result=-result;
773
774 return result;
775}
776
777TF1 NcMath::Chi2Dist(Int_t ndf) const
778{
789
790 TF1 pdf("Chi2PDF","1./(TMath::Gamma([0]/2.)*pow(2,[0]/2.))*pow(x,[0]/2.-1.)*exp(-x/2.)");
791 pdf.SetParName(0,"ndf");
792 pdf.SetParameter(0,ndf);
793 TString title="#chi^{2} PDF (ndf=";
794 title+=ndf;
795 title+=");#chi^{2};p(#chi^{2}|ndf)";
796 pdf.SetTitle(title.Data());
797
798 return pdf;
799}
800
801TF1 NcMath::Chi2CDF(Int_t ndf) const
802{
812
813 TF1 cdf("Chi2CDF","1.-TMath::Prob(x,[0])");
814 cdf.SetParName(0,"ndf");
815 cdf.SetParameter(0,ndf);
816 TString title="#chi^{2} CDF (ndf=";
817 title+=ndf;
818 title+=");#chi^{2};CDF for p(#chi^{2}|ndf)";
819 cdf.SetTitle(title.Data());
820
821 return cdf;
822}
823
824TF1 NcMath::StudentDist(Double_t ndf) const
825{
845
846 TF1 pdf("StudentPDF","(TMath::Gamma(([0]+1.)/2.)/(sqrt(pi*[0])*TMath::Gamma([0]/2.)))*pow(1.+x*x/[0],-([0]+1.)/2.)");
847 pdf.SetParName(0,"ndf");
848 pdf.SetParameter(0,ndf);
849 TString title="Student's T PDF (ndf=";
850 title+=ndf;
851 title+=");T;p(T|ndf)";
852 pdf.SetTitle(title.Data());
853
854 return pdf;
855}
856
857TF1 NcMath::StudentCDF(Double_t ndf) const
858{
877
878 TF1 cdf("StudentCDF","TMath::StudentI(x,[0])");
879 cdf.SetParName(0,"ndf");
880 cdf.SetParameter(0,ndf);
881 TString title="Student's T CDF (ndf=";
882 title+=ndf;
883 title+=");T;CDF for p(T|ndf)";
884 cdf.SetTitle(title.Data());
885
886 return cdf;
887}
888
889TF1 NcMath::FratioDist(Int_t ndf1,Int_t ndf2) const
890{
909
910 TF1 pdf("FratioPDF",
911 "(TMath::Gamma(([0]+[1])/2.)/(TMath::Gamma([0]/2.)*TMath::Gamma([1]/2.)))*pow([0]/[1],[0]/2.)*pow(x,([0]-2.)/2.)/pow(1.+x*[0]/[1],([0]+[1])/2.)");
912 pdf.SetParName(0,"ndf1");
913 pdf.SetParameter(0,ndf1);
914 pdf.SetParName(1,"ndf2");
915 pdf.SetParameter(1,ndf2);
916 TString title="F(ratio) PDF (ndf1=";
917 title+=ndf1;
918 title+=" and ndf2=";
919 title+=ndf2;
920 title+=");F;p(F|ndf1,ndf2)";
921 pdf.SetTitle(title.Data());
922
923 return pdf;
924}
925
926TF1 NcMath::FratioCDF(Int_t ndf1,Int_t ndf2) const
927{
944
945 TF1 cdf("FratioCDF","TMath::FDistI(x,[0],[1])");
946 cdf.SetParName(0,"ndf1");
947 cdf.SetParameter(0,ndf1);
948 cdf.SetParName(1,"ndf2");
949 cdf.SetParameter(1,ndf2);
950 TString title="F(ratio) CDF (ndf1=";
951 title+=ndf1;
952 title+=" and ndf2=";
953 title+=ndf2;
954 title+=");F;CDF for p(F|ndf1,ndf2)";
955 cdf.SetTitle(title.Data());
956
957 return cdf;
958}
959
960TF1 NcMath::BinomialDist(Int_t n,Double_t p) const
961{
976
977 TF1 pdf("BinomialPDF","TMath::Binomial(int([0]),int(x))*pow([1],int(x))*pow(1.-[1],int([0])-int(x))");
978 pdf.SetParName(0,"n");
979 pdf.SetParameter(0,n);
980 pdf.SetParName(1,"p");
981 pdf.SetParameter(1,p);
982 TString title="Binomial PDF for n=";
983 title+=n;
984 title+=" and p=%-10.3g";
985 title+=";k;p(k|n,p)";
986 TString s=title.Format(title.Data(),p);
987 pdf.SetTitle(s.Data());
988
989 return pdf;
990}
991
992TF1 NcMath::BinomialCDF(Int_t n,Double_t p) const
993{
1006
1007 TF1 cdf("BinomialCDF","1.-TMath::BetaIncomplete([1],(x+1.),([0]-x))");
1008 cdf.SetParName(0,"n");
1009 cdf.SetParameter(0,n);
1010 cdf.SetParName(1,"p");
1011 cdf.SetParameter(1,p);
1012 TString title="Binomial CDF for n=";
1013 title+=n;
1014 title+=" and p=%-10.3g";
1015 title+=";k;CDF for p(k|n,p)";
1016 TString s=title.Format(title.Data(),p);
1017 cdf.SetTitle(s.Data());
1018
1019 return cdf;
1020}
1021
1022TF1 NcMath::NegBinomialnDist(Int_t k,Double_t p) const
1023{
1038
1039 TF1 pdf("NegBinomialnPDF","TMath::Binomial(int(x)-1,int([0])-1)*pow([1],int([0]))*pow(1.-[1],int(x)-int([0]))");
1040 pdf.SetParName(0,"k");
1041 pdf.SetParameter(0,k);
1042 pdf.SetParName(1,"p");
1043 pdf.SetParameter(1,p);
1044 TString title="Negative Binomial PDF for k=";
1045 title+=k;
1046 title+=" and p=%-10.3g";
1047 title+=";Number of trials n;p(n|k,p)";
1048 TString s=title.Format(title.Data(),p);
1049 pdf.SetTitle(s.Data());
1050
1051 return pdf;
1052}
1053
1054TF1 NcMath::NegBinomialnCDF(Int_t k,Double_t p) const
1055{
1068
1069 TF1 cdf("NegBinomialnCDF","TMath::BetaIncomplete([1],[0],x-[0]+1.)");
1070 cdf.SetParName(0,"k");
1071 cdf.SetParameter(0,k);
1072 cdf.SetParName(1,"p");
1073 cdf.SetParameter(1,p);
1074 TString title="Negative Binomial CDF for k=";
1075 title+=k;
1076 title+=" and p=%-10.3g";
1077 title+=";Number of trials n;CDF for p(n|k,p)";
1078 TString s=title.Format(title.Data(),p);
1079 cdf.SetTitle(s.Data());
1080
1081 return cdf;
1082}
1083
1084TF1 NcMath::NegBinomialxDist(Int_t k,Double_t p) const
1085{
1102
1103 TF1 pdf("NegBinomialxPDF","TMath::Binomial(int(x)+[0]-1,int([0])-1)*pow([1],int([0]))*pow(1.-[1],int(x))");
1104 pdf.SetParName(0,"k");
1105 pdf.SetParameter(0,k);
1106 pdf.SetParName(1,"p");
1107 pdf.SetParameter(1,p);
1108 TString title="Negative Binomial PDF for k=";
1109 title+=k;
1110 title+=" and p=%-10.3g";
1111 title+=";Number of failures x;p(x|k,p)";
1112 TString s=title.Format(title.Data(),p);
1113 pdf.SetTitle(s.Data());
1114
1115 return pdf;
1116}
1117
1118TF1 NcMath::NegBinomialxCDF(Int_t k,Double_t p) const
1119{
1134
1135 TF1 cdf("NegBinomialxCDF","TMath::BetaIncomplete([1],[0],x+1.)");
1136 cdf.SetParName(0,"k");
1137 cdf.SetParameter(0,k);
1138 cdf.SetParName(1,"p");
1139 cdf.SetParameter(1,p);
1140 TString title="Negative Binomial CDF for k=";
1141 title+=k;
1142 title+=" and p=%-10.3g";
1143 title+=";Number of failures x;CDF for p(x|k,p)";
1144 TString s=title.Format(title.Data(),p);
1145 cdf.SetTitle(s.Data());
1146
1147 return cdf;
1148}
1149
1150TF1 NcMath::PoissonDist(Double_t mu) const
1151{
1165
1166 TF1 pdf("PoissPDFmu","exp(-[0])*pow([0],int(x))/TMath::Factorial(int(x))");
1167 pdf.SetParName(0,"mu");
1168 pdf.SetParameter(0,mu);
1169 TString title="Poisson PDF for #mu=%-10.3g";
1170 title+=";n;p(n|#mu)";
1171 TString s=title.Format(title.Data(),mu);
1172 pdf.SetTitle(s.Data());
1173
1174 return pdf;
1175}
1176
1177TF1 NcMath::PoissonCDF(Double_t mu) const
1178{
1190
1191 TF1 cdf("PoissCDFmu","1.-TMath::Gamma(x,[0])");
1192 cdf.SetParName(0,"mu");
1193 cdf.SetParameter(0,mu);
1194 TString title="Poisson PDF for #mu=%-10.3g";
1195 title+=";n;CDF for p(n|#mu)";
1196 TString s=title.Format(title.Data(),mu);
1197 cdf.SetTitle(s.Data());
1198
1199 return cdf;
1200}
1201
1202TF1 NcMath::PoissonDist(Double_t r,Double_t dt) const
1203{
1217
1218 TF1 pdf("PoissPDFrdt","exp(-[0]*[1])*pow(([0]*[1]),int(x))/TMath::Factorial(int(x))");
1219 pdf.SetParName(0,"r");
1220 pdf.SetParameter(0,r);
1221 pdf.SetParName(1,"dt");
1222 pdf.SetParameter(1,dt);
1223 TString title="Poisson PDF for r=%-10.3g";
1224 title+=" dt=%-10.3g";
1225 title+=";n;p(n|r,dt)";
1226 TString s=title.Format(title.Data(),r,dt);
1227 pdf.SetTitle(s.Data());
1228
1229 return pdf;
1230}
1231
1232TF1 NcMath::PoissonCDF(Double_t r,Double_t dt) const
1233{
1245
1246 TF1 cdf("PoissCDFrdt","1.-TMath::Gamma(x,[0]*[1])");
1247 cdf.SetParName(0,"r");
1248 cdf.SetParameter(0,r);
1249 cdf.SetParName(1,"dt");
1250 cdf.SetParameter(1,dt);
1251 TString title="Poisson CDF for r=%-10.3g";
1252 title+=" dt=%-10.3g";
1253 title+=";n;CDF for p(n|r,dt)";
1254 TString s=title.Format(title.Data(),r,dt);
1255 cdf.SetTitle(s.Data());
1256
1257 return cdf;
1258}
1259
1260TF1 NcMath::PoissonDtDist(Double_t r,Int_t n) const
1261{
1277
1278 TF1 pdf("PoissDtPDF","exp(-[0]*x)*pow(([0]),[1])*pow(x,([1]-1.))/TMath::Factorial(int([1]-1.))");
1279 pdf.SetParName(0,"r");
1280 pdf.SetParameter(0,r);
1281 pdf.SetParName(1,"n");
1282 pdf.SetParameter(1,n);
1283 TString title="Poisson related dt PDF for n=";
1284 title+=n;
1285 title+=" and r=%-10.3g";
1286 title+=";dt;p(dt|r,n)";
1287 TString s=title.Format(title.Data(),r);
1288 pdf.SetTitle(s.Data());
1289
1290 return pdf;
1291}
1292
1293TF1 NcMath::PoissonDtCDF(Double_t r,Int_t n) const
1294{
1308
1309 TF1 cdf("PoissDtCDF","TMath::Gamma([1],[0]*x)");
1310 cdf.SetParName(0,"r");
1311 cdf.SetParameter(0,r);
1312 cdf.SetParName(1,"n");
1313 cdf.SetParameter(1,n);
1314 TString title="Poisson related dt CDF for n=";
1315 title+=n;
1316 title+=" and r=%-10.3g";
1317 title+=";dt;CDF for p(dt|r,n)";
1318 TString s=title.Format(title.Data(),r);
1319 cdf.SetTitle(s.Data());
1320
1321 return cdf;
1322}
1323
1324TF1 NcMath::GammaDtDist(Double_t r,Double_t z) const
1325{
1342
1343 TF1 pdf("GammaDtPDF","exp(-[0]*x)*pow([0],[1])*pow(x,([1]-1.))/TMath::Gamma([1])");
1344 pdf.SetParName(0,"r");
1345 pdf.SetParameter(0,r);
1346 pdf.SetParName(1,"z");
1347 pdf.SetParameter(1,z);
1348 TString title="Gamma related dt PDF for r=%-10.3g";
1349 title+=" z=%-10.3g";
1350 title+=";dt;p(dt|r,z)";
1351 TString s=title.Format(title.Data(),r,z);
1352 pdf.SetTitle(s.Data());
1353
1354 return pdf;
1355}
1356
1357TF1 NcMath::GaussDist(Double_t mu,Double_t sigma) const
1358{
1372
1373 TF1 pdf("GaussPDF","TMath::Gaus(x,[0],[1],1)");
1374 pdf.SetParName(0,"mu");
1375 pdf.SetParameter(0,mu);
1376 pdf.SetParName(1,"sigma");
1377 pdf.SetParameter(1,sigma);
1378 TString title="Gaussian PDF for #mu=%-10.3g";
1379 title+=" #sigma=%-10.3g";
1380 title+=";x;p(x|#mu,#sigma)";
1381 TString s=title.Format(title.Data(),mu,sigma);
1382 pdf.SetTitle(s.Data());
1383
1384 return pdf;
1385}
1386
1387TF1 NcMath::GaussCDF(Double_t mu,Double_t sigma) const
1388{
1400
1401 TF1 cdf("GaussCDF","0.5*(1.+TMath::Erf((x-[0])/([1]*sqrt(2.))))");
1402 cdf.SetParName(0,"mu");
1403 cdf.SetParameter(0,mu);
1404 cdf.SetParName(1,"sigma");
1405 cdf.SetParameter(1,sigma);
1406 TString title="Gaussian CDF for #mu=%-10.3g";
1407 title+=" #sigma=%-10.3g";
1408 title+=";x;CDF for p(x|#mu,#sigma)";
1409 TString s=title.Format(title.Data(),mu,sigma);
1410 cdf.SetTitle(s.Data());
1411
1412 return cdf;
1413}
1414
1415TF1 NcMath::RayleighDist(Double_t sigma) const
1416{
1457
1458 TF1 pdf("RayleighPDF","x*exp(-pow(x,2)/(2.*pow([0],2)))/pow([0],2)");
1459 pdf.SetParName(0,"sigma");
1460 pdf.SetParameter(0,sigma);
1461 TString title;
1462 title.Form("Rayleigh PDF for sigma=%-g;r;p(r|sigma)",sigma);
1463 pdf.SetTitle(title);
1464
1465 return pdf;
1466}
1467
1468TF1 NcMath::RayleighCDF(Double_t sigma) const
1469{
1479
1480 TF1 cdf("RayleighCDF","1.-exp(-pow(x,2)/(2.*pow([0],2)))");
1481 cdf.SetParName(0,"sigma");
1482 cdf.SetParameter(0,sigma);
1483 TString title;
1484 title.Form("Rayleigh CDF for sigma=%-g;r;CDF for p(r|sigma)",sigma);
1485 cdf.SetTitle(title);
1486
1487 return cdf;
1488}
1489
1490TF1 NcMath::Rayleigh3Dist(Double_t sigma) const
1491{
1530
1531 TF1 pdf("3D-RayleighPDF","2*pow(x,2)*pow(2*pi,-0.5)*exp(-pow(x,2)/(2.*pow([0],2)))/pow([0],3)");
1532 pdf.SetParName(0,"sigma");
1533 pdf.SetParameter(0,sigma);
1534 TString title;
1535 title.Form("3D-Rayleigh PDF for sigma=%-g;r;p(r|sigma)",sigma);
1536 pdf.SetTitle(title);
1537
1538 return pdf;
1539}
1540
1541Double_t NcMath::GetStatistic(TF1 f,TString name,Int_t n,Double_t vref,Int_t npx) const
1542{
1575
1576 // Set number of evaluation points to ensure accurate results
1577 f.SetNpx(npx);
1578
1579 Double_t xmin=f.GetXmin();
1580 Double_t xmax=f.GetXmax();
1581
1582 Double_t val=0;
1583
1584 if (name=="mode") val=f.GetMaximumX(xmin,xmax);
1585 if (name=="mean") val=f.Mean(xmin,xmax);
1586 if (name=="var") val=f.Variance(xmin,xmax);
1587 if (name=="std")
1588 {
1589 val=f.Variance(xmin,xmax);
1590 if (val>0)
1591 {
1592 val=sqrt(val);
1593 }
1594 else
1595 {
1596 val=0;
1597 }
1598 }
1599 if (name=="median")
1600 {
1601 Double_t p[1]={0.5};
1602 Double_t q[1];
1603 Int_t nq=f.GetQuantiles(1,q,p);
1604 if (nq) val=q[0];
1605 }
1606 if (name=="mom" || name=="cmom")
1607 {
1608 if (n<=0)
1609 {
1610 printf(" *NcMath::GetStatistic* Inconsistent input n=%-i for statistic %-s. \n",n,name.Data());
1611 return 0;
1612 }
1613 if (name=="mom") val=f.Moment(n,xmin,xmax);
1614 if (name=="cmom") val=f.CentralMoment(n,xmin,xmax);
1615 }
1616 if (name=="skew" || name=="kurt")
1617 {
1618 Double_t std=GetStatistic(f,"std");
1619 if (name=="skew" && std>0) val=f.CentralMoment(3,xmin,xmax)/pow(std,3);
1620 if (name=="kurt" && std>0) val=-3.+f.CentralMoment(4,xmin,xmax)/pow(std,4);
1621 }
1622 if (name=="spread")
1623 {
1624 if (n<0 || n>2)
1625 {
1626 printf(" *NcMath::GetStatistic* Inconsistent input n=%-i for statistic %-s. \n",n,name.Data());
1627 return 0;
1628 }
1629 if (n==0) vref=GetStatistic(f,"median");
1630 if (n==1) vref=f.Mean(xmin,xmax);
1631 Double_t step=fabs(xmax-xmin)/double(npx);
1632 Double_t x=xmin;
1633 Double_t weight=0;
1634 Double_t wn=0;
1635 val=0;
1636 while (x<=xmax)
1637 {
1638 weight=fabs(f.Eval(x));
1639 val+=weight*fabs(vref-x);
1640 x+=step;
1641 wn+=weight;
1642 }
1643 if (wn) val=val/wn;
1644 }
1645
1646 return val;
1647}
1648
1649TGraph NcMath::GetCDF(TF1 f,Int_t npx) const
1650{
1663
1664 f.SetNpx(npx);
1665
1666 TGraph gr(&f,"i");
1667
1668 TString ftitle=f.GetTitle();
1669
1670 // Modify the (possible) specifications of the axis titles of "f"
1671 TString saxes=";x;CDF";
1672 Int_t j1=ftitle.First(';');
1673 if (j1>0)
1674 {
1675 saxes=";CDF";
1676 Int_t j2=ftitle.Last(';');
1677 if (j2>j1)
1678 {
1679 ftitle.Insert(j2+1,"CDF for ");
1680 saxes="";
1681 }
1682 }
1683
1684 TString gtitle;
1685 if (ftitle !="")
1686 {
1687 gtitle.Form("Original function : %-s%-s",ftitle.Data(),saxes.Data());
1688 }
1689 else
1690 {
1691 gtitle.Form("Original function : f(x)=%-s%-s",(f.GetExpFormula("p")).Data(),saxes.Data());
1692 }
1693
1694 gr.SetTitle(gtitle);
1695 gr.SetDrawOption("AL");
1696
1697 return gr;
1698}
1699
1700Double_t NcMath::GaussProb(Double_t q,Double_t mean,Double_t sigma,Int_t isig) const
1701{
1719
1720 Double_t val=-1;
1721 if (!isig)
1722 {
1723 val=Erf(fabs(q-mean)/(sigma*sqrt(2.)));
1724 }
1725 else
1726 {
1727 val=Erf(fabs(q)/sqrt(2.));
1728 }
1729 return val;
1730}
1731
1732Double_t NcMath::GaussPvalue(Double_t q,Double_t mean,Double_t sigma,Int_t sides,Int_t isig) const
1733{
1769
1770 Double_t val=-1;
1771 if (!isig)
1772 {
1773 val=Erfc(fabs(q-mean)/(sigma*sqrt(2.)));
1774 }
1775 else
1776 {
1777 if (isig==1) val=Erfc(fabs(q)/sqrt(2.));
1778 if (isig==-1) val=(q-mean)/sigma;
1779 }
1780 if (sides==1 && isig!=-1) val/=2.;
1781 return val;
1782}
1783
1784Double_t NcMath::Chi2Pvalue(Double_t chi2,Int_t ndf,Int_t sides,Int_t sigma,Int_t mode) const
1785{
1835
1836 if (ndf<=0) return 0;
1837
1838 Double_t mean=ndf;
1839
1840 if (!sides) // Automatic one-sided test
1841 {
1842 sides=1;
1843 if (chi2<mean) sides=-1;
1844 }
1845
1846 if (sides==3) // Automatic two-sided test
1847 {
1848 sides=2;
1849 if (chi2<mean) sides=-2;
1850 }
1851
1852 Double_t val=0;
1853 if (sigma) // P-value in units of sigma
1854 {
1855 Double_t s=sqrt(double(2*ndf));
1856 val=(chi2-mean)/s;
1857 }
1858 else // P-value from tail contents
1859 {
1860 if (sides>0) // Upper tail
1861 {
1862 val=Prob(chi2,ndf,mode);
1863 }
1864 else // Lower tail
1865 {
1866 val=1.-Prob(chi2,ndf,mode);
1867 }
1868 }
1869
1870 if (abs(sides)==2) val=val*2.;
1871
1872 return val;
1873}
1874
1875Double_t NcMath::StudentPvalue(Double_t t,Double_t ndf,Int_t sides,Int_t sigma) const
1876{
1920
1921 if (ndf<=0) return 0;
1922
1923 Double_t mean=0;
1924
1925 if (!sides) // Automatic one-sided test
1926 {
1927 sides=1;
1928 if (t<mean) sides=-1;
1929 }
1930
1931 if (sides==3) // Automatic two-sided test
1932 {
1933 sides=2;
1934 if (t<mean) sides=-2;
1935 }
1936
1937 Double_t val=0;
1938 if (sigma) // P-value in units of sigma
1939 {
1940 if (ndf>2) // Sigma is only defined for ndf>2
1941 {
1942 Double_t s=sqrt(ndf/(ndf-2.));
1943 val=t/s;
1944 }
1945 }
1946 else // P-value from tail contents
1947 {
1948 if (sides>0) // Upper tail
1949 {
1950 val=1.-TMath::StudentI(t,ndf);
1951 }
1952 else // Lower tail
1953 {
1954 val=TMath::StudentI(t,ndf);
1955 }
1956 }
1957
1958 if (abs(sides)==2) val=val*2.;
1959
1960 return val;
1961}
1962
1963Double_t NcMath::FratioPvalue(Double_t f,Int_t ndf1,Int_t ndf2,Int_t sides,Int_t sigma) const
1964{
2012
2013 if (ndf1<=0 || ndf2<=0 || f<=0) return 0;
2014
2015 Double_t mean=double(ndf2/(ndf2-2));
2016
2017 if (!sides) // Automatic one-sided test
2018 {
2019 sides=1;
2020 if (f<mean) sides=-1;
2021 }
2022
2023 if (sides==3) // Automatic two-sided test
2024 {
2025 sides=2;
2026 if (f<mean) sides=-2;
2027 }
2028
2029 Double_t val=0;
2030 if (sigma) // P-value in units of sigma
2031 {
2032 // Sigma is only defined for ndf2>4
2033 if (ndf2>4)
2034 {
2035 Double_t s=sqrt(double(ndf2*ndf2*(2*ndf2+2*ndf1-4))/double(ndf1*pow(double(ndf2-1),2)*(ndf2-4)));
2036 val=(f-mean)/s;
2037 }
2038 }
2039 else // P-value from tail contents
2040 {
2041 if (sides>0) // Upper tail
2042 {
2043 val=1.-TMath::FDistI(f,ndf1,ndf2);
2044 }
2045 else // Lower tail
2046 {
2047 val=TMath::FDistI(f,ndf1,ndf2);
2048 }
2049 }
2050
2051 if (abs(sides)==2) val=val*2.;
2052
2053 return val;
2054}
2055
2056Double_t NcMath::BinomialPvalue(Int_t k,Int_t n,Double_t p,Int_t sides,Int_t sigma,Int_t mode) const
2057{
2101
2102 Double_t mean=double(n)*p;
2103
2104 if (!sides) // Automatic one-sided test
2105 {
2106 sides=1;
2107 if (k<mean) sides=-1;
2108 }
2109
2110 if (sides==3) // Automatic two-sided test
2111 {
2112 sides=2;
2113 if (k<mean) sides=-2;
2114 }
2115
2116 Double_t val=0;
2117
2118 if (sigma) // P-value in units of sigma
2119 {
2120 Double_t s=sqrt(double(n)*p*(1.-p));
2121 val=(double(k)-mean)/s;
2122 }
2123 else // P-value from tail contents
2124 {
2125 if (sides>0) // Upper tail
2126 {
2127 if (!mode) // Use Incomplete Beta function
2128 {
2129 val=TMath::BetaIncomplete(p,k,n-k+1);
2130 }
2131 else // Use straightforward summation
2132 {
2133 for (Int_t i=k; i<=n; i++)
2134 {
2135 val+=TMath::Binomial(n,i)*pow(p,i)*pow(1.-p,n-i);
2136 }
2137 }
2138 }
2139 else // Lower tail
2140 {
2141 if (!mode) // Use Incomplete Beta function
2142 {
2143 val=1.-TMath::BetaIncomplete(p,k+1,n-k);
2144 }
2145 else
2146 {
2147 for (Int_t j=0; j<=k; j++)
2148 {
2149 val+=TMath::Binomial(n,j)*pow(p,j)*pow(1.-p,n-j);
2150 }
2151 }
2152 }
2153 }
2154
2155 if (abs(sides)==2) val=val*2.;
2156
2157 return val;
2158}
2159
2160Double_t NcMath::PoissonPvalue(Int_t k,Double_t mu,Int_t sides,Int_t sigma) const
2161{
2205
2206 Double_t mean=mu;
2207
2208 if (!sides) // Automatic one-sided test
2209 {
2210 sides=1;
2211 if (k<mean) sides=-1;
2212 }
2213
2214 if (sides==3) // Automatic two-sided test
2215 {
2216 sides=2;
2217 if (k<mean) sides=-2;
2218 }
2219
2220 Double_t val=0;
2221
2222 if (sigma) // P-value in units of sigma
2223 {
2224 Double_t s=sqrt(mu);
2225 val=(double(k)-mean)/s;
2226 }
2227 else // P-value from tail contents
2228 {
2229 if (sides>0) // Upper tail
2230 {
2231 val=Gamma(k-1,mu);
2232 }
2233 else // Lower tail
2234 {
2235 val=1.-Gamma(k,mu);
2236 }
2237 }
2238
2239 if (abs(sides)==2) val=val*2.;
2240
2241 return val;
2242}
2243
2244Double_t NcMath::PoissonPvalue(Int_t k,Double_t r,Double_t dt,Int_t sides,Int_t sigma) const
2245{
2289
2290 Double_t mu=r*dt;
2291 Double_t val=PoissonPvalue(k,mu,sides,sigma);
2292 return val;
2293}
2294
2295Double_t NcMath::PoissonDtPvalue(Double_t dt,Double_t r,Int_t n,Int_t sides,Int_t sigma) const
2296{
2342
2343 Double_t val=-1;
2344
2345 if (n<=0 || r<=0) return val;
2346
2347 Double_t mean=double(n)/r;
2348 Double_t sig=sqrt(double(n)/(r*r));
2349
2350 if (!sides) // Automatic one-sided test
2351 {
2352 sides=1;
2353 if (dt<mean) sides=-1;
2354 }
2355
2356 if (sides==3) // Automatic two-sided test
2357 {
2358 sides=2;
2359 if (dt<mean) sides=-2;
2360 }
2361
2362 val=0;
2363
2364 if (sigma) // P-value in units of sigma
2365 {
2366 val=(dt-mean)/sig;
2367 }
2368 else // P-value from tail contents
2369 {
2370 if (sides>0) // Upper tail
2371 {
2372 val=1.-Gamma(n,r*dt);
2373 }
2374 else // Lower tail
2375 {
2376 val=Gamma(n,r*dt);
2377 }
2378 }
2379
2380 if (abs(sides)==2) val=val*2.;
2381
2382 return val;
2383}
2384
2385Double_t NcMath::NegBinomialnPvalue(Int_t n,Int_t k,Double_t p,Int_t sides,Int_t sigma,Int_t mode) const
2386{
2434
2435 Double_t mean=double(k)/p;
2436
2437 if (!sides) // Automatic one-sided test
2438 {
2439 sides=1;
2440 if (n<mean) sides=-1;
2441 }
2442
2443 if (sides==3) // Automatic two-sided test
2444 {
2445 sides=2;
2446 if (n<mean) sides=-2;
2447 }
2448
2449 Double_t val=0;
2450
2451 if (sigma) // P-value in units of sigma
2452 {
2453 Double_t s=sqrt(double(k)*(1.-p)/(p*p));
2454 val=(double(n)-mean)/s;
2455 }
2456 else // P-value from tail contents
2457 {
2458 if (sides>0) // Upper tail
2459 {
2460 if (!mode) // Use Incomplete Beta function
2461 {
2462 TF1 cdf=NegBinomialnCDF(k,p);
2463 val=cdf.Eval(double(n-1));
2464 }
2465 else // Use straightforward summation
2466 {
2467 val=0;
2468 for (Int_t i=1; i<n; i++)
2469 {
2470 val+=TMath::Binomial(i-1,k-1)*pow(p,k)*pow(1.-p,i-k);
2471 }
2472 }
2473 val=1.-val;
2474 }
2475 else // Lower tail
2476 {
2477 if (!mode) // Use Incomplete Beta function
2478 {
2479 TF1 cdf=NegBinomialnCDF(k,p);
2480 val=cdf.Eval(double(n));
2481 }
2482 else // Use straightforward summation
2483 {
2484 val=0;
2485 for (Int_t j=1; j<=n; j++)
2486 {
2487 val+=TMath::Binomial(j-1,k-1)*pow(p,k)*pow(1.-p,j-k);
2488 }
2489 }
2490 }
2491 }
2492
2493 if (abs(sides)==2) val=val*2.;
2494
2495 return val;
2496}
2497
2498Double_t NcMath::NegBinomialxPvalue(Int_t x,Int_t k,Double_t p,Int_t sides,Int_t sigma,Int_t mode) const
2499{
2549
2550 Double_t mean=double(x)*(1.-p)/p;
2551
2552 if (!sides) // Automatic one-sided test
2553 {
2554 sides=1;
2555 if (x<mean) sides=-1;
2556 }
2557
2558 if (sides==3) // Automatic two-sided test
2559 {
2560 sides=2;
2561 if (x<mean) sides=-2;
2562 }
2563
2564 Double_t val=0;
2565
2566 if (sigma) // P-value in units of sigma
2567 {
2568 Double_t s=sqrt(double(k)*(1.-p)/(p*p));
2569 val=(double(x)-mean)/s;
2570 }
2571 else // P-value from tail contents
2572 {
2573 if (sides>0) // Upper tail
2574 {
2575 if (!mode) // Use Incomplete Beta function
2576 {
2577 TF1 cdf=NegBinomialxCDF(k,p);
2578 val=cdf.Eval(double(x-1));
2579 }
2580 else // Use straightforward summation
2581 {
2582 val=0;
2583 for (Int_t i=0; i<x; i++)
2584 {
2585 val+=TMath::Binomial(i+k-1,k-1)*pow(p,k)*pow(1.-p,i);
2586 }
2587 }
2588 val=1.-val;
2589 }
2590 else // Lower tail
2591 {
2592 if (!mode) // Use Incomplete Beta function
2593 {
2594 TF1 cdf=NegBinomialxCDF(k,p);
2595 val=cdf.Eval(double(x));
2596 }
2597 else // Use straightforward summation
2598 {
2599 val=0;
2600 for (Int_t j=0; j<=x; j++)
2601 {
2602 val+=TMath::Binomial(j+k-1,k-1)*pow(p,k)*pow(1.-p,j);
2603 }
2604 }
2605 }
2606 }
2607
2608 if (abs(sides)==2) val=val*2.;
2609
2610 return val;
2611}
2612
2613Double_t NcMath::Log(Double_t B,Double_t x) const
2614{
2622
2623 Double_t val=0;
2624
2625 if (B<=1 || x<=0) return val;
2626
2627 val=log(x)/log(B);
2628 return val;
2629}
2630
2631Double_t NcMath::Nfac(Int_t n,Int_t mode) const
2632{
2651
2652 if (n<0) return 0;
2653 if (n==0 || n==1) return 1;
2654
2655 Double_t twopi=2.*acos(-1.);
2656 Double_t z=0;
2657 Double_t nfac=1;
2658 Int_t i=n;
2659
2660 switch (mode)
2661 {
2662 case 0: // Straightforward multiplication
2663 while (i>1)
2664 {
2665 nfac*=Double_t(i);
2666 i--;
2667 }
2668 break;
2669
2670 case 1: // Stirling's approximation
2671 z=n;
2672 nfac=sqrt(twopi)*pow(z,z+0.5)*exp(-z)*(1.+1./(12.*z));
2673 break;
2674
2675 case 2: // Use of Gamma(n+1)
2676 z=n+1;
2677 nfac=Gamma(z);
2678 break;
2679
2680 default:
2681 nfac=0;
2682 break;
2683 }
2684
2685 return nfac;
2686}
2687
2688Double_t NcMath::LnNfac(Int_t n,Int_t mode) const
2689{
2708
2709 if (n<=1) return 0;
2710
2711 Double_t twopi=2.*acos(-1.);
2712 Double_t z=0;
2713 Double_t lognfac=0;
2714
2715 switch (mode)
2716 {
2717 case 0: // Straightforward ln(n!)
2718 z=Nfac(n);
2719 lognfac=log(z);
2720 break;
2721
2722 case 1: // Stirling's approximation
2723 z=n;
2724 lognfac=0.5*log(twopi)+(z+0.5)*log(z)-z+1./(12.*z);
2725 break;
2726
2727 case 2: // Use of LnGamma(n+1)
2728 z=n+1;
2729 lognfac=LnGamma(z);
2730 break;
2731
2732 default:
2733 lognfac=0;
2734 break;
2735 }
2736
2737 return lognfac;
2738}
2739
2740Double_t NcMath::LogNfac(Int_t n,Int_t mode) const
2741{
2753
2754 if (n<=1) return 0;
2755
2756 Double_t e=exp(1.);
2757
2758 Double_t val=LnNfac(n,mode);
2759 val*=log10(e);
2760
2761 return val;
2762}
2763
2764Double_t NcMath::Rfac(Double_t r) const
2765{
2778
2779 if (r<0) return 0;
2780 if (r==0 || r==1) return 1;
2781
2782 Double_t z=r+1.;
2783 return Gamma(z);
2784}
2785
2786Double_t NcMath::LnRfac(Double_t r) const
2787{
2798
2799 if (r<=0 || r==1) return 0;
2800
2801 Double_t z=r+1.;
2802 return LnGamma(z);
2803}
2804
2805Double_t NcMath::LogRfac(Double_t r) const
2806{
2818
2819 if (r<=0 || r==1) return 0;
2820
2821 Double_t e=exp(1.);
2822
2823 Double_t val=LnRfac(r);
2824 val*=log10(e);
2825
2826 return val;
2827}
2828
2829Double_t NcMath::PsiValue(Int_t m,Int_t* n,Double_t* p,Int_t f) const
2830{
2871
2872 Double_t psi=-1;
2873
2874 if (m<=0 || !n) return psi;
2875
2876 Int_t ntot=0;
2877 for (Int_t j=0; j<m; j++)
2878 {
2879 if (n[j]>0) ntot+=n[j];
2880 }
2881
2882 psi=0;
2883 Double_t pk=1./float(m); // Prob. of getting outcome k for a uniform distr.
2884 for (Int_t i=0; i<m; i++)
2885 {
2886 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
2887 if (n[i]>0 && pk>0)
2888 {
2889 if (!f) // Exact Bayesian expression
2890 {
2891 psi+=double(n[i])*log10(pk)-LogNfac(n[i]);
2892 }
2893 else // Frequentist (Stirling) approximation
2894 {
2895 if (ntot>0) psi+=double(n[i])*log10(double(n[i])/(double(ntot)*pk));
2896 }
2897 }
2898 }
2899
2900 if (!f) // Exact Bayesian expression
2901 {
2902 psi+=LogNfac(ntot);
2903 psi*=-10.;
2904 }
2905 else // Frequentist (Stirling) approximation
2906 {
2907 psi*=10.;
2908 }
2909
2910 return psi;
2911}
2912
2913Double_t NcMath::PsiValue(Int_t m,Double_t* n,Double_t* p,Int_t f) const
2914{
2960
2961 Double_t psi=-1;
2962
2963 if (m<=0 || !n) return psi;
2964
2965 Double_t ntot=0;
2966 for (Int_t j=0; j<m; j++)
2967 {
2968 if (n[j]>0) ntot+=n[j];
2969 }
2970
2971 psi=0;
2972 Double_t pk=1./float(m); // Prob. of getting outcome k for a uniform distr.
2973 for (Int_t i=0; i<m; i++)
2974 {
2975 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
2976 if (n[i]>0 && pk>0)
2977 {
2978 if (!f) // Exact Bayesian expression
2979 {
2980 psi+=n[i]*log10(pk)-LogRfac(n[i]);
2981 }
2982 else // Frequentist (Stirling) approximation
2983 {
2984 if (ntot>0) psi+=n[i]*log10(n[i]/(ntot*pk));
2985 }
2986 }
2987 }
2988
2989 if (!f) // Exact Bayesian expression
2990 {
2991 psi+=LogRfac(ntot);
2992 psi*=-10.;
2993 }
2994 else // Frequentist (Stirling) approximation
2995 {
2996 psi*=10.;
2997 }
2998
2999 return psi;
3000}
3001
3002Double_t NcMath::PsiValue(TH1* his,TH1* hyp,TF1* pdf,Int_t f) const
3003{
3048
3049 Double_t psi=-1;
3050
3051 if (!his) return psi;
3052
3053 TAxis* xaxis=his->GetXaxis();
3054 Double_t xmin=xaxis->GetXmin();
3055 Double_t xmax=xaxis->GetXmax();
3056 Double_t range=xmax-xmin;
3057 Int_t nbins=his->GetNbinsX();
3058 Double_t nensig=his->GetSumOfWeights();
3059
3060 if (nbins<=0 || nensig<=0 || range<=0) return psi;
3061
3062 Double_t* n=new Double_t[nbins];
3063 Double_t* p=new Double_t[nbins];
3064 Double_t nk;
3065 Double_t pk;
3066
3067 // Uniform hypothesis distribution
3068 if (!hyp && !pdf)
3069 {
3070 for (Int_t i=1; i<=nbins; i++)
3071 {
3072 nk=his->GetBinContent(i);
3073 pk=his->GetBinWidth(i)/range;
3074 n[i-1]=0;
3075 p[i-1]=0;
3076 if (nk>0) n[i-1]=nk;
3077 if (pk>0) p[i-1]=pk;
3078 }
3079 psi=PsiValue(nbins,n,p,f);
3080 }
3081
3082 // Hypothesis specified via a pdf
3083 if (pdf && !hyp)
3084 {
3085 pdf->SetRange(xmin,xmax);
3086 Double_t ftot=pdf->Integral(xmin,xmax);
3087 if (ftot>0)
3088 {
3089 Double_t x1,x2;
3090 for (Int_t ipdf=1; ipdf<=nbins; ipdf++)
3091 {
3092 nk=his->GetBinContent(ipdf);
3093 x1=his->GetBinLowEdge(ipdf);
3094 x2=x1+his->GetBinWidth(ipdf);
3095 pk=pdf->Integral(x1,x2)/ftot;
3096 n[ipdf-1]=0;
3097 p[ipdf-1]=0;
3098 if (nk>0) n[ipdf-1]=nk;
3099 if (pk>0) p[ipdf-1]=pk;
3100 }
3101 psi=PsiValue(nbins,n,p,f);
3102 }
3103 }
3104
3105 // Hypothesis specified via a histogram
3106 if (hyp && !pdf)
3107 {
3108 TH1* href=(TH1*)his->Clone("href");
3109 href->Reset();
3110 Double_t nenhyp=0;
3111 Double_t x,y;
3112 for (Int_t ihyp=1; ihyp<=hyp->GetNbinsX(); ihyp++)
3113 {
3114 x=hyp->GetBinCenter(ihyp);
3115 y=hyp->GetBinContent(ihyp);
3116 href->Fill(x,y);
3117 nenhyp+=y;
3118 }
3119 for (Int_t j=1; j<=nbins; j++)
3120 {
3121 nk=his->GetBinContent(j);
3122 pk=href->GetBinContent(j)/nenhyp;
3123 n[j-1]=0;
3124 p[j-1]=0;
3125 if (nk>0) n[j-1]=nk;
3126 if (pk>0) p[j-1]=pk;
3127 }
3128 psi=PsiValue(nbins,n,p,f);
3129 delete href;
3130 }
3131
3132 delete [] n;
3133 delete [] p;
3134
3135 return psi;
3136}
3137
3138Double_t NcMath::PsiExtreme(Double_t n,Int_t m,Double_t* p,Int_t k) const
3139{
3195
3196 Double_t psi=-1;
3197
3198 if (m<=0 || n<=0 || k<-2 || k>m) return psi;
3199
3200 psi=LogRfac(n);
3201 Double_t pk=1./double(m); // Prob. of getting outcome k for a uniform distr.
3202
3204 // The user specified fixed outcome case //
3206 if (k>0)
3207 {
3208 if (p) pk=p[k-1];
3209 if (pk>0)
3210 {
3211 psi+=n*log10(pk)-LogRfac(n);
3212 psi*=-10.;
3213 }
3214 else
3215 {
3216 psi=-1;
3217 }
3218 return psi;
3219 }
3220
3222 // Determine the minimal and maximal probability of all the possible outcomes //
3224 Double_t pmin=999;
3225 Double_t pmax=-1;
3226 Int_t jmin=-1;
3227 Int_t jmax=-1;
3228 if (!p) // Uniform distribution
3229 {
3230 pmin=pk;
3231 pmax=pk;
3232 jmin=0;
3233 jmax=0;
3234 }
3235 else // User specified distribution
3236 {
3237 for (Int_t j=0; j<m; j++)
3238 {
3239 if (p[j]>0 && p[j]<pmin)
3240 {
3241 pmin=p[j];
3242 jmin=j;
3243 }
3244 if (p[j]>0 && p[j]>pmax)
3245 {
3246 pmax=p[j];
3247 jmax=j;
3248 }
3249 }
3250 }
3251
3252 // Check for validity of the encountered pmin and pmax
3253 if (jmin<0 || jmax<0) return psi;
3254
3256 // The worst matching case //
3258 if (k==-1)
3259 {
3260 if (pmin>0)
3261 {
3262 psi+=n*log10(pmin)-LogRfac(n);
3263 psi*=-10.;
3264 }
3265 else
3266 {
3267 psi=-1;
3268 }
3269 return psi;
3270 }
3271
3272 Double_t nk=0;
3273
3275 // The best matching case for fractional nk values //
3277 if (k==0)
3278 {
3279 for (Int_t i=0; i<m; i++)
3280 {
3281 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3282 if (pk>0)
3283 {
3284 nk=n*pk;
3285 psi+=nk*log10(pk)-LogRfac(nk);
3286 }
3287 }
3288 psi*=-10.;
3289 return psi;
3290 }
3291
3293 // The best matching case for integer nk values //
3295 if (k==-2)
3296 {
3297 // Determine the best matching discrete distribution
3298 // by starting from the fractional nk=n*pk distribution
3299 Double_t* narr=new Double_t[m];
3300 Double_t ndisc=0;
3301 Int_t ndiff=0;
3302 for (Int_t i=0; i<m; i++)
3303 {
3304 narr[i]=0;
3305 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3306 if (pk>0) narr[i]=n*pk;
3307 ndisc+=int(narr[i]);
3308 }
3309
3310 // Check the (integer) number of entries
3311 ndiff=int(n-ndisc);
3312
3313 // Continue to complete the best matching discrete distribution
3314 // in case the integer number of entries doesn't match the original "n".
3315 while (ndiff>0)
3316 {
3317 // For a uniform distribution just fill the "m bins" one after the other
3318 // with the remaining entries
3319 if (jmin==jmax) // Signature for a uniform distr.
3320 {
3321 Int_t im=0;
3322 for (Int_t ient=0; ient<ndiff; ient++)
3323 {
3324 narr[im]++;
3325 im++;
3326 if (im==m) im=0;
3327 }
3328 ndiff=0;
3329 }
3330 else // Increase the various "bin contents" until the total number of entries fits
3331 {
3332 ndisc=0;
3333 for (Int_t i=0; i<m; i++)
3334 {
3335 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3336 if (pk>0) narr[i]+=pk;
3337 ndisc+=int(narr[i]);
3338 }
3339 ndiff=int(n-ndisc);
3340 // Only 1 entry left to fill --> Put it at the max. probability
3341 if (ndiff==1)
3342 {
3343 narr[jmax]+=1;
3344 ndiff=0;
3345 }
3346 }
3347 } // End of while loop
3348
3349 // In case too many entries have been filled, remove (one by one) the ones
3350 // with the lowest probability.
3351 while (ndiff<0)
3352 {
3353 // Determine the filled bin with minimal pk
3354 pmin=999;
3355 jmin=0;
3356 for (Int_t i=0; i<m; i++)
3357 {
3358 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3359 if (int(narr[i])>0 && pk<pmin)
3360 {
3361 pmin=p[i];
3362 jmin=i;
3363 }
3364 }
3365 // Remove one entry
3366 narr[jmin]-=1;
3367 ndiff++;
3368 }
3369
3370 // The best matching discrete distr. is complete now
3371 for (Int_t i=0; i<m; i++)
3372 {
3373 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3374 if (pk>0)
3375 {
3376 nk=int(narr[i]);
3377 psi+=nk*log10(pk)-LogRfac(nk);
3378 }
3379 }
3380 psi*=-10.;
3381
3382 delete [] narr;
3383 }
3384 return psi;
3385}
3386
3387Double_t NcMath::PsiExtreme(TH1* his,TH1* hyp,TF1* pdf,Int_t k) const
3388{
3449
3450 Double_t psi=-1;
3451
3452 if (!his) return psi;
3453
3454 TAxis* xaxis=his->GetXaxis();
3455 Double_t xmin=xaxis->GetXmin();
3456 Double_t xmax=xaxis->GetXmax();
3457 Double_t range=xmax-xmin;
3458 Int_t nbins=his->GetNbinsX();
3459 Double_t nensig=his->GetSumOfWeights();
3460
3461 if (nbins<=0 || nensig<=0 || range<=0) return psi;
3462
3463 Double_t* p=new Double_t[nbins];
3464 Double_t pk;
3465
3466 // Uniform hypothesis distribution
3467 if (!hyp && !pdf)
3468 {
3469 for (Int_t i=1; i<=nbins; i++)
3470 {
3471 pk=his->GetBinWidth(i)/range;
3472 p[i-1]=0;
3473 if (pk>0) p[i-1]=pk;
3474 }
3475 psi=PsiExtreme(nensig,nbins,p,k);
3476 }
3477
3478 // Hypothesis specified via a pdf
3479 if (pdf && !hyp)
3480 {
3481 pdf->SetRange(xmin,xmax);
3482 Double_t ftot=pdf->Integral(xmin,xmax);
3483 if (ftot>0)
3484 {
3485 Double_t x1,x2;
3486 for (Int_t ipdf=1; ipdf<=nbins; ipdf++)
3487 {
3488 x1=his->GetBinLowEdge(ipdf);
3489 x2=x1+his->GetBinWidth(ipdf);
3490 pk=pdf->Integral(x1,x2)/ftot;
3491 p[ipdf-1]=0;
3492 if (pk>0) p[ipdf-1]=pk;
3493 }
3494 psi=PsiExtreme(nensig,nbins,p,k);
3495 }
3496 }
3497
3498 // Hypothesis specified via a histogram
3499 if (hyp && !pdf)
3500 {
3501 TH1* href=(TH1*)his->Clone("href");
3502 href->Reset();
3503 Double_t nenhyp=0;
3504 Double_t x=0;
3505 Double_t y=0;
3506 for (Int_t ihyp=1; ihyp<=hyp->GetNbinsX(); ihyp++)
3507 {
3508 x=hyp->GetBinCenter(ihyp);
3509 y=hyp->GetBinContent(ihyp);
3510 href->Fill(x,y);
3511 nenhyp+=y;
3512 }
3513 for (Int_t j=1; j<=nbins; j++)
3514 {
3515 pk=href->GetBinContent(j)/nenhyp;
3516 p[j-1]=0;
3517 if (pk>0) p[j-1]=pk;
3518 }
3519 psi=PsiExtreme(nensig,nbins,p,k);
3520 delete href;
3521 }
3522
3523 delete [] p;
3524
3525 return psi;
3526}
3527
3528Double_t NcMath::PsiPvalue(Double_t psi0,Double_t nr,Double_t n,Int_t m,Double_t* p,Int_t f,Double_t* na,TH1F* psih,Int_t ncut,Double_t* nrx,Int_t mark)
3529{
3619
3620 if (psi0<0 || nr<0 || (nr>0 && nr<2)) return -1;
3621
3622 Double_t pval=-1;
3623
3624 NcRandom q;
3625 if (na) n=-n;
3626 Double_t nrused;
3627 Double_t val=q.RanBm(nr,n,m,p,na,0,psi0,f,psih,ncut,&nrused);
3628 if (val>=0) pval=val/nrused;
3629 if (nrx) (*nrx)=nrused;
3630
3631 // Set axes titles for the "psih" histogram
3632 if (psih)
3633 {
3634 TString title="#psi value in dB";
3635 TString s="Counts after ";
3636 s+=nrused;
3637 s+=" randomisations";
3638
3639 psih->GetXaxis()->SetTitle(title.Data());
3640 psih->GetYaxis()->SetTitle(s.Data());
3641
3642 // Mark the psi0 value by a vertical line in the "psih" histogram
3643 // Also the corresponding P-value is mentioned in the legend
3644 if (mark)
3645 {
3646 Float_t x=psi0;
3647 Float_t ymin=0;
3648 Float_t ymax=psih->GetMaximum();
3649
3650 TLine* vline=new TLine(x,ymin,x,ymax);
3651 vline->SetLineStyle(2); // Dashed line
3652 vline->SetLineWidth(2);
3653 vline->SetLineColor(4); // Blue color
3654
3655 title="P-value : %-10.3g";
3656 s=title.Format(title.Data(),pval);
3657
3658 TLegend* leg=new TLegend(0.6,0.8,0.8,0.9);
3659 leg->SetFillColor(0);
3660 leg->SetHeader(s.Data());
3661 leg->AddEntry(vline,"Observed #psi","L");
3662
3663 TList* hlist=psih->GetListOfFunctions();
3664 hlist->Add(vline);
3665 hlist->Add(leg);
3666 }
3667 }
3668
3669
3670 return pval;
3671}
3672
3673Double_t NcMath::PsiPvalue(Double_t psi0,Double_t nr,TH1* his,TH1* hyp,TF1* pdf,Int_t f,Double_t* na,TH1F* psih,Int_t ncut,Double_t* nrx,Int_t mark)
3674{
3775
3776 Double_t pval=-1;
3777
3778 if (!his) return pval;
3779
3780 TAxis* xaxis=his->GetXaxis();
3781 Double_t xmin=xaxis->GetXmin();
3782 Double_t xmax=xaxis->GetXmax();
3783 Double_t range=xmax-xmin;
3784 Int_t nbins=his->GetNbinsX();
3785 Double_t nensig=his->GetSumOfWeights();
3786
3787 if (nbins<=0 || nensig<=0 || range<=0) return pval;
3788
3789 if (psi0<0) psi0=PsiValue(his,hyp,pdf,f);
3790
3791 Double_t* n=new Double_t[nbins];
3792 Double_t* p=new Double_t[nbins];
3793 Double_t nk;
3794 Double_t pk;
3795
3796 // Uniform hypothesis distribution
3797 if (!hyp && !pdf)
3798 {
3799 for (Int_t i=1; i<=nbins; i++)
3800 {
3801 nk=his->GetBinContent(i);
3802 pk=his->GetBinWidth(i)/range;
3803 n[i-1]=0;
3804 p[i-1]=0;
3805 if (nk>0) n[i-1]=nk;
3806 if (pk>0) p[i-1]=pk;
3807 }
3808 pval=PsiPvalue(psi0,nr,nensig,nbins,p,f,na,psih,ncut,nrx,mark);
3809 }
3810
3811 // Hypothesis specified via a pdf
3812 if (pdf && !hyp)
3813 {
3814 pdf->SetRange(xmin,xmax);
3815 Double_t ftot=pdf->Integral(xmin,xmax);
3816 if (ftot>0)
3817 {
3818 Double_t x1,x2;
3819 for (Int_t ipdf=1; ipdf<=nbins; ipdf++)
3820 {
3821 nk=his->GetBinContent(ipdf);
3822 x1=his->GetBinLowEdge(ipdf);
3823 x2=x1+his->GetBinWidth(ipdf);
3824 pk=pdf->Integral(x1,x2)/ftot;
3825 n[ipdf-1]=0;
3826 p[ipdf-1]=0;
3827 if (nk>0) n[ipdf-1]=nk;
3828 if (pk>0) p[ipdf-1]=pk;
3829 }
3830 pval=PsiPvalue(psi0,nr,nensig,nbins,p,f,na,psih,ncut,nrx,mark);
3831 }
3832 }
3833
3834 // Hypothesis specified via a histogram
3835 if (hyp && !pdf)
3836 {
3837 TH1* href=(TH1*)his->Clone("href");
3838 href->Reset();
3839 Double_t nenhyp=0;
3840 Double_t x,y;
3841 for (Int_t ihyp=1; ihyp<=hyp->GetNbinsX(); ihyp++)
3842 {
3843 x=hyp->GetBinCenter(ihyp);
3844 y=hyp->GetBinContent(ihyp);
3845 href->Fill(x,y);
3846 nenhyp+=y;
3847 }
3848 for (Int_t j=1; j<=nbins; j++)
3849 {
3850 nk=his->GetBinContent(j);
3851 pk=href->GetBinContent(j)/nenhyp;
3852 n[j-1]=0;
3853 p[j-1]=0;
3854 if (nk>0) n[j-1]=nk;
3855 if (pk>0) p[j-1]=pk;
3856 }
3857 pval=PsiPvalue(psi0,nr,nensig,nbins,p,f,na,psih,ncut,nrx,mark);
3858 delete href;
3859 }
3860
3861 delete [] n;
3862 delete [] p;
3863
3864 return pval;
3865}
3866
3867Double_t NcMath::Chi2Value(Int_t m,Int_t* n,Double_t* p,Int_t* ndf) const
3868{
3897
3898 Double_t chi=-1;
3899 if (ndf) *ndf=-1;
3900
3901 if (m<=0 || !n) return chi;
3902
3903 Int_t ntot=0;
3904 for (Int_t j=0; j<m; j++)
3905 {
3906 if (n[j]>0) ntot+=n[j];
3907 }
3908
3909 chi=0;
3910 Double_t pk=1./float(m); // Prob. of getting outcome k for a uniform distr.
3911 for (Int_t i=0; i<m; i++)
3912 {
3913 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3914 if (n[i]>0 && pk>0 && ntot>0)
3915 {
3916 chi+=pow(double(n[i])-double(ntot)*pk,2)/(double(ntot)*pk);
3917 if (ndf) (*ndf)=(*ndf)+1;
3918 }
3919 }
3920
3921 return chi;
3922}
3923
3924Double_t NcMath::Chi2Value(Int_t m,Double_t* n,Double_t* p,Int_t* ndf) const
3925{
3959
3960 Double_t chi=-1;
3961 if (ndf) *ndf=-1;
3962
3963 if (m<=0 || !n) return chi;
3964
3965 Double_t ntot=0;
3966 for (Int_t j=0; j<m; j++)
3967 {
3968 if (n[j]>0) ntot+=n[j];
3969 }
3970
3971 chi=0;
3972 Double_t pk=1./float(m); // Prob. of getting outcome k for a uniform distr.
3973 for (Int_t i=0; i<m; i++)
3974 {
3975 if (p) pk=p[i]; // Probabilities from a specific B_m hypothesis
3976 if (n[i]>0 && pk>0 && ntot>0)
3977 {
3978 chi+=pow(n[i]-ntot*pk,2)/(ntot*pk);
3979 if (ndf) (*ndf)=(*ndf)+1;
3980 }
3981 }
3982
3983 return chi;
3984}
3985
3986Double_t NcMath::Chi2Value(TH1* his,TH1* hyp,TF1* pdf,Int_t* ndf) const
3987{
4019
4020 Double_t chi=-1;
4021
4022 if (!his) return chi;
4023
4024 TAxis* xaxis=his->GetXaxis();
4025 Double_t xmin=xaxis->GetXmin();
4026 Double_t xmax=xaxis->GetXmax();
4027 Double_t range=xmax-xmin;
4028 Int_t nbins=his->GetNbinsX();
4029 Double_t nensig=his->GetSumOfWeights();
4030
4031 if (nbins<=0 || nensig<=0 || range<=0) return chi;
4032
4033 Double_t* n=new Double_t[nbins];
4034 Double_t* p=new Double_t[nbins];
4035 Double_t nk;
4036 Double_t pk;
4037
4038 // Uniform hypothesis distribution
4039 if (!hyp && !pdf)
4040 {
4041 for (Int_t i=1; i<=nbins; i++)
4042 {
4043 nk=his->GetBinContent(i);
4044 pk=his->GetBinWidth(i)/range;
4045 n[i-1]=0;
4046 p[i-1]=0;
4047 if (nk>0) n[i-1]=nk;
4048 if (pk>0) p[i-1]=pk;
4049 }
4050 chi=Chi2Value(nbins,n,p,ndf);
4051 }
4052
4053 // Hypothesis specified via a pdf
4054 if (pdf && !hyp)
4055 {
4056 pdf->SetRange(xmin,xmax);
4057 Double_t ftot=pdf->Integral(xmin,xmax);
4058 if (ftot>0)
4059 {
4060 Double_t x1,x2;
4061 for (Int_t ipdf=1; ipdf<=nbins; ipdf++)
4062 {
4063 nk=his->GetBinContent(ipdf);
4064 x1=his->GetBinLowEdge(ipdf);
4065 x2=x1+his->GetBinWidth(ipdf);
4066 pk=pdf->Integral(x1,x2)/ftot;
4067 n[ipdf-1]=0;
4068 p[ipdf-1]=0;
4069 if (nk>0) n[ipdf-1]=nk;
4070 if (pk>0) p[ipdf-1]=pk;
4071 }
4072 chi=Chi2Value(nbins,n,p,ndf);
4073 }
4074 }
4075
4076 // Hypothesis specified via a histogram
4077 if (hyp && !pdf)
4078 {
4079 TH1* href=(TH1*)his->Clone("href");
4080 href->Reset();
4081 Double_t nenhyp=0;
4082 Double_t x,y;
4083 for (Int_t ihyp=1; ihyp<=hyp->GetNbinsX(); ihyp++)
4084 {
4085 x=hyp->GetBinCenter(ihyp);
4086 y=hyp->GetBinContent(ihyp);
4087 href->Fill(x,y);
4088 nenhyp+=y;
4089 }
4090 for (Int_t j=1; j<=nbins; j++)
4091 {
4092 nk=his->GetBinContent(j);
4093 pk=href->GetBinContent(j)/nenhyp;
4094 n[j-1]=0;
4095 p[j-1]=0;
4096 if (nk>0) n[j-1]=nk;
4097 if (pk>0) p[j-1]=pk;
4098 }
4099 chi=Chi2Value(nbins,n,p,ndf);
4100 delete href;
4101 }
4102
4103 delete [] n;
4104 delete [] p;
4105
4106 return chi;
4107}
4108
4109Double_t NcMath::MeanMu(Double_t cl,Double_t nbkg,Int_t mode,TF1* fw,TFeldmanCousins* fc,Int_t nmax) const
4110{
4136
4137 TFeldmanCousins* f=fc;
4138 if (!f) f=new TFeldmanCousins();
4139
4140 f->SetCL(cl);
4141
4142 if (!nmax) nmax=int(nbkg+10.*sqrt(nbkg));
4143
4144 Double_t mu;
4145 Double_t muav=0;
4146
4147 // Cumulative summation for the number of observed events
4148 for (Int_t nobs=0; nobs<=nmax; nobs++)
4149 {
4150 if (mode==1)
4151 {
4152 mu=f->CalculateLowerLimit(nobs,nbkg);
4153 }
4154 else
4155 {
4156 mu=f->CalculateUpperLimit(nobs,nbkg);
4157 }
4158 if (fw)
4159 {
4160 muav+=mu*fw->Eval(nobs);
4161 }
4162 else // Poisson pdf weight
4163 {
4164 muav+=mu*pow(nbkg,nobs)*exp(-nbkg)/Rfac(nobs);
4165 }
4166 }
4167
4168 if (!fc) delete f;
4169 return muav;
4170}
4171
4172Double_t NcMath::LiMaSignificance(Double_t Non,Double_t Ton,Double_t Noff,Double_t Toff,Double_t Ra,Double_t Re) const
4173{
4200
4201 Double_t s=-1;
4202
4203 if (Non<=0 || Noff<=0 || Ton<=0 || Toff<=0 || Ra<=0 || Re<=0) return -1;
4204
4205 Double_t sum=Non+Noff;
4206
4207 Double_t a=Ra*Re*Ton/Toff;
4208
4209 s=2.*(Non*log((1.+a)*Non/(a*sum))+Noff*log((1.+a)*Noff/sum));
4210 s=sqrt(s);
4211 return s;
4212}
4213
ClassImp
ClassImp(NcMath)
NcMath
Various mathematical tools for scientific analysis.
Definition NcMath.h:26
NcMath::StudentPvalue
Double_t StudentPvalue(Double_t t, Double_t ndf, Int_t sides=0, Int_t sigma=0) const
Definition NcMath.cxx:1875
NcMath::PsiValue
Double_t PsiValue(Int_t m, Int_t *n, Double_t *p=0, Int_t f=0) const
Definition NcMath.cxx:2829
NcMath::PoissonCDF
TF1 PoissonCDF(Double_t mu) const
Definition NcMath.cxx:1177
NcMath::GetStatistic
Double_t GetStatistic(TF1 f, TString name, Int_t n=0, Double_t vref=0, Int_t npx=1000) const
Definition NcMath.cxx:1541
NcMath::BesselK
Double_t BesselK(Int_t n, Double_t x) const
Definition NcMath.cxx:678
NcMath::Log
Double_t Log(Double_t B, Double_t x) const
Definition NcMath.cxx:2613
NcMath::Chi2CDF
TF1 Chi2CDF(Int_t ndf) const
Definition NcMath.cxx:801
NcMath::NegBinomialnPvalue
Double_t NegBinomialnPvalue(Int_t n, Int_t k, Double_t p, Int_t sides=0, Int_t sigma=0, Int_t mode=0) const
Definition NcMath.cxx:2385
NcMath::StudentDist
TF1 StudentDist(Double_t ndf) const
Definition NcMath.cxx:824
NcMath::GaussPvalue
Double_t GaussPvalue(Double_t q, Double_t mean=0, Double_t sigma=1, Int_t sides=2, Int_t isig=0) const
Definition NcMath.cxx:1732
NcMath::PoissonDist
TF1 PoissonDist(Double_t mu) const
Definition NcMath.cxx:1150
NcMath::PoissonDtCDF
TF1 PoissonDtCDF(Double_t r, Int_t n) const
Definition NcMath.cxx:1293
NcMath::LnGamma
Double_t LnGamma(Double_t z) const
Definition NcMath.cxx:193
NcMath::BinomialPvalue
Double_t BinomialPvalue(Int_t k, Int_t n, Double_t p, Int_t sides=0, Int_t sigma=0, Int_t mode=0) const
Definition NcMath.cxx:2056
NcMath::Chi2Pvalue
Double_t Chi2Pvalue(Double_t chi2, Int_t ndf, Int_t sides=0, Int_t sigma=0, Int_t mode=1) const
Definition NcMath.cxx:1784
NcMath::Chi2Dist
TF1 Chi2Dist(Int_t ndf) const
Definition NcMath.cxx:777
NcMath::NegBinomialxCDF
TF1 NegBinomialxCDF(Int_t k, Double_t p) const
Definition NcMath.cxx:1118
NcMath::PoissonDtPvalue
Double_t PoissonDtPvalue(Double_t dt, Double_t r, Int_t n, Int_t sides=0, Int_t sigma=0) const
Definition NcMath.cxx:2295
NcMath::BinomialCDF
TF1 BinomialCDF(Int_t n, Double_t p) const
Definition NcMath.cxx:992
NcMath::Nfac
Double_t Nfac(Int_t n, Int_t mode=0) const
Definition NcMath.cxx:2631
NcMath::MeanMu
Double_t MeanMu(Double_t cl, Double_t nbkg, Int_t mode, TF1 *w=0, TFeldmanCousins *f=0, Int_t nmax=0) const
Definition NcMath.cxx:4109
NcMath::BesselI0
Double_t BesselI0(Double_t x) const
Definition NcMath.cxx:497
NcMath::PoissonPvalue
Double_t PoissonPvalue(Int_t k, Double_t mu, Int_t sides=0, Int_t sigma=0) const
Definition NcMath.cxx:2160
NcMath::LiMaSignificance
Double_t LiMaSignificance(Double_t Non, Double_t Ton, Double_t Noff, Double_t Toff, Double_t Ra=1, Double_t Re=1) const
Definition NcMath.cxx:4172
NcMath::Gamma
Double_t Gamma(Double_t z) const
Definition NcMath.cxx:121
NcMath::LnRfac
Double_t LnRfac(Double_t r) const
Definition NcMath.cxx:2786
NcMath::GammaDtDist
TF1 GammaDtDist(Double_t r, Double_t z) const
Definition NcMath.cxx:1324
NcMath::GaussProb
Double_t GaussProb(Double_t q, Double_t mean=0, Double_t sigma=1, Int_t isig=0) const
Definition NcMath.cxx:1700
NcMath::Erfc
Double_t Erfc(Double_t x) const
Definition NcMath.cxx:383
NcMath::BinomialDist
TF1 BinomialDist(Int_t n, Double_t p) const
Definition NcMath.cxx:960
NcMath::NcMath
NcMath()
Definition NcMath.cxx:61
NcMath::BesselI
Double_t BesselI(Int_t n, Double_t x) const
Definition NcMath.cxx:720
NcMath::FratioPvalue
Double_t FratioPvalue(Double_t f, Int_t ndf1, Int_t ndf2, Int_t sides=0, Int_t sigma=0) const
Definition NcMath.cxx:1963
NcMath::LnNfac
Double_t LnNfac(Int_t n, Int_t mode=2) const
Definition NcMath.cxx:2688
NcMath::Rfac
Double_t Rfac(Double_t r) const
Definition NcMath.cxx:2764
NcMath::Zeta
Double_t Zeta(Double_t x, Int_t nterms=100000) const
Definition NcMath.cxx:88
NcMath::Rayleigh3Dist
TF1 Rayleigh3Dist(Double_t sigma) const
Definition NcMath.cxx:1490
NcMath::PsiExtreme
Double_t PsiExtreme(Double_t n, Int_t m, Double_t *p=0, Int_t k=0) const
Definition NcMath.cxx:3138
NcMath::NegBinomialxDist
TF1 NegBinomialxDist(Int_t k, Double_t p) const
Definition NcMath.cxx:1084
NcMath::BesselK0
Double_t BesselK0(Double_t x) const
Definition NcMath.cxx:540
NcMath::RayleighCDF
TF1 RayleighCDF(Double_t sigma) const
Definition NcMath.cxx:1468
NcMath::~NcMath
virtual ~NcMath()
Definition NcMath.cxx:70
NcMath::GetCDF
TGraph GetCDF(TF1 f, Int_t npx=1000) const
Definition NcMath.cxx:1649
NcMath::LogRfac
Double_t LogRfac(Double_t r) const
Definition NcMath.cxx:2805
NcMath::Prob
Double_t Prob(Double_t chi2, Int_t ndf, Int_t mode=1) const
Definition NcMath.cxx:421
NcMath::GaussCDF
TF1 GaussCDF(Double_t mu, Double_t sigma) const
Definition NcMath.cxx:1387
NcMath::NegBinomialxPvalue
Double_t NegBinomialxPvalue(Int_t x, Int_t k, Double_t p, Int_t sides=0, Int_t sigma=0, Int_t mode=0) const
Definition NcMath.cxx:2498
NcMath::BesselK1
Double_t BesselK1(Double_t x) const
Definition NcMath.cxx:631
NcMath::Erf
Double_t Erf(Double_t x) const
Definition NcMath.cxx:370
NcMath::RayleighDist
TF1 RayleighDist(Double_t sigma) const
Definition NcMath.cxx:1415
NcMath::GaussDist
TF1 GaussDist(Double_t mu, Double_t sigma) const
Definition NcMath.cxx:1357
NcMath::NegBinomialnCDF
TF1 NegBinomialnCDF(Int_t k, Double_t p) const
Definition NcMath.cxx:1054
NcMath::NegBinomialnDist
TF1 NegBinomialnDist(Int_t k, Double_t p) const
Definition NcMath.cxx:1022
NcMath::PsiPvalue
Double_t PsiPvalue(Double_t psi0, Double_t nr, Double_t n, Int_t m, Double_t *p=0, Int_t f=0, Double_t *na=0, TH1F *psih=0, Int_t ncut=0, Double_t *nrx=0, Int_t mark=1)
Definition NcMath.cxx:3528
NcMath::FratioCDF
TF1 FratioCDF(Int_t ndf1, Int_t ndf2) const
Definition NcMath.cxx:926
NcMath::PoissonDtDist
TF1 PoissonDtDist(Double_t r, Int_t n) const
Definition NcMath.cxx:1260
NcMath::LogNfac
Double_t LogNfac(Int_t n, Int_t mode=2) const
Definition NcMath.cxx:2740
NcMath::FratioDist
TF1 FratioDist(Int_t ndf1, Int_t ndf2) const
Definition NcMath.cxx:889
NcMath::GamSer
Double_t GamSer(Double_t a, Double_t x) const
Definition NcMath.cxx:271
NcMath::StudentCDF
TF1 StudentCDF(Double_t ndf) const
Definition NcMath.cxx:857
NcMath::Chi2Value
Double_t Chi2Value(Int_t m, Int_t *n, Double_t *p=0, Int_t *ndf=0) const
Definition NcMath.cxx:3867
NcMath::GamCf
Double_t GamCf(Double_t a, Double_t x) const
Definition NcMath.cxx:316
NcMath::BesselI1
Double_t BesselI1(Double_t x) const
Definition NcMath.cxx:587
NcRandom
Generate universal random numbers and sequences on all common machines.
Definition NcRandom.h:16
NcRandom::RanBm
Double_t RanBm(Double_t nr, Double_t n, Int_t m, Double_t *p=0, Double_t *na=0, Double_t *psia=0, Double_t psi0=-1, Int_t f=0, TH1 *psih=0, Int_t ncut=0, Double_t *nrx=0)
Definition NcRandom.cxx:1143