ROOT  6.06/09
Reference Guide
testTStatistic.cxx
Go to the documentation of this file.
1 
2 #include "TRandom3.h"
3 #include "TStatistic.h"
4 #include "TStopwatch.h"
5 #include <iostream>
6 
7 bool gVerbose = false;
8 
9 void testTStatistic(Int_t n = 10000)
10 {
11 
12  // Array
13  std::vector<Double_t> xx(n);
14  std::vector<Double_t> ww(n);
15 
16  double true_mean = 10;
17  double true_sigma = 1;
18  double eps1 = 5*true_sigma/sqrt(double(n)); // 5 times error on the mean
19  double eps2 = 5*true_sigma/sqrt(double(2.*n)); // 5 times error on the RMS
20 
21  // Gaussian first
22  TRandom3 rand3;
23  for (Int_t i = 0; i < n; i++) {
24  xx[i] = rand3.Gaus(true_mean, true_sigma);
25  ww[i] = rand3.Uniform(0.01, 3.00);
26  }
27 
28 
29  TStopwatch stp;
30  bool error = false;
31 
32  printf("\nTest without using weights : ");
33 
34  TStatistic st0("st0", n, xx.data());
35  stp.Stop();
36  if (!TMath::AreEqualAbs(st0.GetMean(),true_mean, eps1) ) { Error("TestTStatistic-GetMean","Different value obtained for the unweighted data"); error = true; }
37  if (!TMath::AreEqualAbs(st0.GetRMS(),true_sigma, eps2) ) { Error("TestTStatistic-GetRMS","Different value obtained for the unweighted data"); error = true; }
38 
39  if (error) printf("Failed\n");
40  else printf("OK\n");
41  if (error || gVerbose) {
42  stp.Print();
43  st0.Print();
44  }
45 
46  // test with TMath
47  printf("\nTest using TMath: ");
48  error = false;
49  stp.Start();
50  double mean = TMath::Mean(xx.begin(), xx.end() );
51  double rms = TMath::RMS(xx.begin(), xx.end() );
52  stp.Stop();
53  if (!TMath::AreEqualAbs(mean,true_mean, eps1) ) { Error("TestTStatistic::TMath::Mean","Different value obtained for the unweighted data"); error = true; }
54  if (!TMath::AreEqualAbs(rms,true_sigma, eps2) ) { Error("TestTStatistic::TMath::RMS","Different value obtained for the unweighted data"); error = true; }
55 
56  if (error) printf("Failed\n");
57  else printf("OK\n");
58  if (error || gVerbose) {
59  stp.Print();
60  printf(" TMATH mu = %.5g +- %.4g \t RMS = %.5g \n",mean, rms/sqrt(double(xx.size()) ), rms);
61  }
62 
63 
64  printf("\nTest using Weights : ");
65  error = false;
66  stp.Start();
67  TStatistic st1("st1", n, xx.data(), ww.data());
68  stp.Stop();
69 
70  if (!TMath::AreEqualAbs(st1.GetMean(),true_mean, eps1) ) { Error("TestTStatistic-GetMean","Different value obtained for the weighted data"); error = true; }
71  if (!TMath::AreEqualAbs(st1.GetRMS(),true_sigma, eps2) ) { Error("TestTStatistic-GetRMS","Different value obtained for the weighted data"); error = true; }
72 
73  if (error) printf("Failed\n");
74  else printf("OK\n");
75  if (error || gVerbose) {
76  stp.Print();
77  st1.Print();
78  }
79 
80 
81  // Incremental test
82  printf("\nTest incremental filling : ");
83  error = false;
84  TStatistic st2("st2");
85  stp.Start();
86  for (Int_t i = 0; i < n; i++) {
87  st2.Fill(xx[i], ww[i]);
88  }
89  stp.Stop();
90 
91  if (!TMath::AreEqualRel(st1.GetMean(),st2.GetMean(), 1.E-15) ) { Error("TestTStatistic-GetMean","2 Different values obtained for the weighted data"); error = true; }
92  if (!TMath::AreEqualRel(st1.GetRMS(),st2.GetRMS(), 1.E-15) ) { Error("TestTStatistic-GetRMS","2 Different values obtained for the weighted data"); error = true; }
93 
94  if (error) printf("Failed\n");
95  else printf("OK\n");
96  if (error || gVerbose) {
97  stp.Print();
98  st2.Print();
99  }
100 
101 
102  // test merge
103  int n1 = rand3.Uniform(10,n-10);
104 
105  // sort the data to have then two biased samples
106  std::sort(xx.begin(), xx.end() );
107 
108 
109  printf("\nTest merge : ");
110  error = false;
111  TStatistic sta("sta");
112  TStatistic stb("stb");
113  stp.Start();
114  for (int i = 0; i < n ; ++i) {
115  if (i < n1) sta.Fill(xx[i],ww[i] );
116  else stb.Fill(xx[i],ww[i] );
117  }
118 
119  TList l; l.Add(&stb);
120  sta.Merge(&l);
121  stp.Stop();
122 
123  if (!TMath::AreEqualAbs(sta.GetMean(),true_mean, eps1) ) { Error("TestTStatistic-GetMean","Different value obtained for the merged data"); error = true; }
124  if (!TMath::AreEqualAbs(sta.GetRMS(),true_sigma, eps2) ) { Error("TestTStatistic-GetRMS","Different value obtained for the merged data"); error = true; }
125 
126  if (error) printf("Failed\n");
127  else printf("OK\n");
128  if (error || gVerbose) {
129  stp.Print();
130  sta.Print();
131  }
132 
133  printf("\nTest sorted data : ");
134  error = false;
135  stp.Start();
136  TStatistic st3("st3", n, xx.data(), ww.data());
137  stp.Stop();
138 
139  if (!TMath::AreEqualAbs(st3.GetMean(), sta.GetMean(), 1.E-10 ) ) { Error("TestTStatistic-GetMean","Different value obtained for the sorted data"); error = true; }
140  if (!TMath::AreEqualAbs(st3.GetRMS(), sta.GetRMS() , 1.E-10 ) ) { Error("TestTStatistic-GetRMS","Different value obtained for the sorted data"); error = true; }
141 
142  if (error) printf("Failed\n");
143  else printf("OK\n");
144  if (error || gVerbose) {
145  stp.Print();
146  st3.Print();
147  }
148 
149 
150  // test with TMath
151  printf("\nTest TMath with weights : ");
152  error = false;
153  stp.Start();
154  double meanw = TMath::Mean(xx.begin(), xx.end(), ww.begin() );
155  double rmsw = TMath::RMS(xx.begin(), xx.end(), ww.begin() );
156  double neff = st2.GetW() * st2.GetW() / st2.GetW2();
157  stp.Stop();
158 
159  if (!TMath::AreEqualAbs(meanw,true_mean, eps1) ) { Error("TestTStatistic::TMath::Mean","Different value obtained for the weighted data"); error = true; }
160  if (!TMath::AreEqualAbs(rmsw,true_sigma, eps2) ) { Error("TestTStatistic::TMath::RMS","Different value obtained for the weighted data"); error = true; }
161 
162  if (error) printf("Failed\n");
163  else printf("OK\n");
164  if (error || gVerbose) {
165  stp.Print();
166  printf(" TMATH mu = %.5g +- %.4g \t RMS = %.5g \n",meanw, rmsw/sqrt(neff ), rmsw);
167  }
168 
169 
170 }
171 
172 int main(int argc, char **argv)
173 {
174  // Parse command line arguments
175  for (Int_t i=1 ; i<argc ; i++) {
176  std::string arg = argv[i] ;
177  if (arg == "-v") {
178  gVerbose = true;
179  }
180  if (arg == "-h") {
181  std::cout << "Usage: " << argv[0] << " [-v]\n";
182  std::cout << " where:\n";
183  std::cout << " -v : verbose mode";
184  std::cout << std::endl;
185  return -1;
186  }
187  }
188 
189  testTStatistic();
190 
191  return 0;
192 
193 }
bool gVerbose
void Print(Option_t *="") const
This method must be overridden when a class wants to print itself.
Definition: TStatistic.cxx:76
void Print(Option_t *option="") const
Print the real and cpu time passed between the start and stop events.
Definition: TStopwatch.cxx:217
Random number generator class based on M.
Definition: TRandom3.h:29
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
int Int_t
Definition: RtypesCore.h:41
Double_t GetW2() const
Definition: TStatistic.h:73
Double_t RMS(Long64_t n, const T *a, const Double_t *w=0)
Definition: TMath.h:916
double sqrt(double)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
Statistical variable, defined by its mean and variance (RMS).
Definition: TStatistic.h:45
Double_t GetMean() const
Definition: TStatistic.h:68
void Error(const char *location, const char *msgfmt,...)
A doubly linked list.
Definition: TList.h:47
Bool_t AreEqualAbs(Double_t af, Double_t bf, Double_t epsilon)
Definition: TMath.h:192
Bool_t AreEqualRel(Double_t af, Double_t bf, Double_t relPrec)
Definition: TMath.h:196
Double_t Mean(Long64_t n, const T *a, const Double_t *w=0)
Definition: TMath.h:824
TLine * l
Definition: textangle.C:4
int main(int argc, char **argv)
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
void Fill(Double_t val, Double_t w=1.)
Definition: TStatistic.cxx:48
Double_t GetW() const
Definition: TStatistic.h:72
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
Int_t Merge(TCollection *in)
Definition: TStatistic.cxx:85
void testTStatistic(Int_t n=10000)
virtual void Add(TObject *obj)
Definition: TList.h:81
Double_t GetRMS() const
Definition: TStatistic.h:70
const Int_t n
Definition: legend1.C:16
Stopwatch class.
Definition: TStopwatch.h:30