NuStatistics Class Reference

#include <NuStatistics.h>

List of all members.

Public Member Functions

 NuStatistics ()
virtual ~NuStatistics ()
void DataSpectrum (NuMatrixSpectrum *spec)
void PredSpectrum (NuMatrixSpectrum *spec)
void FitPars (int fp)
int FitPars ()
void SetMaxEnergy (double en)
double KolmogorovShape ()
double ChisqCount ()
double KolmogorovBoth ()
double ChisqShape ()
int ChisqDOF ()
double ChisqProb ()
double GetMaxE (bool compress=false)
void ApplyMaxE (TH1 *h1, bool compress=false)
double CombineSignificance (double prb1, double prb2)

Private Member Functions

 ClassDef (NuStatistics, 0)

Private Attributes

int fitPars
double maxE
NuMatrixSpectrumdata
NuMatrixSpectrumpred
NuMatrixSpectrumdata_rb
NuMatrixSpectrumpred_rb

Detailed Description

Definition at line 7 of file NuStatistics.h.


Constructor & Destructor Documentation

NuStatistics::NuStatistics (  ) 

Definition at line 16 of file NuStatistics.cxx.

References data, data_rb, pred, and pred_rb.

00017 : fitPars(0), maxE(50.)
00018 {
00019     data = 0;
00020     data_rb = 0;
00021     pred = 0;
00022     pred_rb = 0;
00023 }

NuStatistics::~NuStatistics (  )  [virtual]

Definition at line 27 of file NuStatistics.cxx.

References data, data_rb, pred, and pred_rb.

00028 {
00029     if (data) delete data;
00030     if (data_rb) delete data_rb;
00031     if (pred) delete pred;
00032     if (pred_rb) delete pred_rb;  
00033 }


Member Function Documentation

void NuStatistics::ApplyMaxE ( TH1 *  h1,
bool  compress = false 
)

Definition at line 324 of file NuStatistics.cxx.

References GetMaxE(), and max.

Referenced by DataSpectrum(), and PredSpectrum().

00325 {
00326     double max = GetMaxE(compress);
00327     int ifirst = 0; 
00328     int ilast = h1->GetNbinsX()+1;
00329     for (int bin = ifirst; bin <= ilast; bin++) {
00330         if (h1->GetBinLowEdge(bin) >= max) {
00331             // Make sure excluded bins do not contribute
00332             h1->SetBinContent(bin, 0);
00333             h1->SetBinError(bin, 0);
00334         }
00335     }
00336 }

double NuStatistics::ChisqCount (  ) 

Definition at line 210 of file NuStatistics.cxx.

References data, pred, and NuMatrixSpectrum::Spectrum().

Referenced by KolmogorovBoth().

00211 {
00212     TH1D *h1 = data->Spectrum();
00213     TH1D *h2 = pred->Spectrum();
00214     
00215     double sum1 = 0, sum2 = 0;
00216     double ew1, ew2, w1 = 0, w2 = 0;
00217     
00218 
00219     for (int bin = 1; bin <= h1->GetNbinsX(); bin++) {
00220         sum1 += h1->GetBinContent(bin);
00221         sum2 += h2->GetBinContent(bin);
00222         ew1   = h1->GetBinError(bin);
00223         ew2   = h2->GetBinError(bin);
00224         w1   += ew1*ew1;
00225         w2   += ew2*ew2;
00226     }
00227     
00228     // My simple chisquared method
00229     double delta = sum1 - sum2;
00230     double sigmasq = w1 + w2;
00231     double chi2 = delta*delta / sigmasq;
00232     double prob = TMath::Prob(chi2, 1);
00233     
00234     printf(" Chisq Norm  h1 = %s, sum bin content =%g  sumsq error =%g\n",h1->GetName(),sum1,w1);
00235     printf(" Chisq Norm  h2 = %s, sum bin content =%g  sumsq error =%g\n",h2->GetName(),sum2,w2);
00236     printf(" Chisq Norm     = %g\n",prob);
00237     
00238     return prob;
00239 }

int NuStatistics::ChisqDOF (  ) 

Definition at line 286 of file NuStatistics.cxx.

References FitPars(), GetMaxE(), pred_rb, and NuMatrixSpectrum::Spectrum().

Referenced by ChisqProb().

00287 {
00288     TH1D *h1 = pred_rb->Spectrum();
00289     
00290     int dof = -1 * FitPars();
00291     for (int i = 1; i <= h1->GetNbinsX(); i++){
00292         if (h1->GetBinLowEdge(i) >= GetMaxE(true) ) break;
00293         ++dof;
00294     }
00295     
00296     return dof;
00297 }

double NuStatistics::ChisqProb (  ) 

Definition at line 301 of file NuStatistics.cxx.

References ChisqDOF(), and ChisqShape().

00302 {
00303     double chi2 = ChisqShape();
00304     int dof = ChisqDOF();
00305     return TMath::Prob(chi2, dof);
00306 }

double NuStatistics::ChisqShape (  ) 

Definition at line 254 of file NuStatistics.cxx.

References data_rb, pred_rb, and NuMatrixSpectrum::Spectrum().

Referenced by ChisqProb().

00255 {
00256     TH1D *h1 = data_rb->Spectrum();
00257     TH1D *h2 = pred_rb->Spectrum();
00258     
00259     double chisq = 0;
00260     double da, mc, eda, emc;
00261     double delta, sigmasq;
00262 
00263     for (int i=1; i<= h1->GetNbinsX(); i++){
00264         da = h1->GetBinContent(i);
00265         mc = h2->GetBinContent(i);
00266         
00267         eda = h1->GetBinError(i);
00268         emc = h2->GetBinError(i);
00269         
00270         delta = da - mc;
00271         sigmasq = eda*eda + emc*emc;
00272         
00273         if (sigmasq <= 0) continue;
00274         
00275         chisq += delta*delta/sigmasq;
00276     }
00277     
00278     printf(" Chisq Shape  chi2 = %g\n",chisq);
00279     
00280     
00281     return chisq;
00282 }

NuStatistics::ClassDef ( NuStatistics  ,
 
) [private]
double NuStatistics::CombineSignificance ( double  prb1,
double  prb2 
)

Definition at line 340 of file NuStatistics.cxx.

Referenced by KolmogorovBoth().

00341 {
00342     // see Eadie et al., section 11.6.2
00343     
00344     double prob;
00345     if (prb1 > 0 && prb2 > 0) prob = prb1*prb2*(1-TMath::Log(prb1*prb2));
00346     else                      prob = 0;    
00347     
00348     return prob;
00349 }

void NuStatistics::DataSpectrum ( NuMatrixSpectrum spec  ) 

Definition at line 37 of file NuStatistics.cxx.

References ApplyMaxE(), NuMatrixSpectrum::BlessedFD(), data, data_rb, NuMatrixSpectrum::RemoveErrors(), and NuMatrixSpectrum::Spectrum().

00038 {
00039     if (data) delete data;
00040     data = new NuMatrixSpectrum(*spec);
00041     data->RemoveErrors();
00042     data->Spectrum()->SetName("Data");
00043     ApplyMaxE(data->Spectrum());
00044     
00045     if (data_rb) delete data_rb;
00046     data_rb = new NuMatrixSpectrum(*spec);
00047     data_rb->BlessedFD();
00048     data_rb->RemoveErrors();
00049     data_rb->Spectrum()->SetName("DataDisplay");
00050     ApplyMaxE(data_rb->Spectrum(), true);
00051 }

int NuStatistics::FitPars (  )  [inline]

Definition at line 18 of file NuStatistics.h.

References fitPars.

Referenced by ChisqDOF().

00018 {return fitPars;};

void NuStatistics::FitPars ( int  fp  )  [inline]

Definition at line 17 of file NuStatistics.h.

References fitPars.

00017 {fitPars = fp;};

double NuStatistics::GetMaxE ( bool  compress = false  ) 

Definition at line 311 of file NuStatistics.cxx.

References NuMatrixSpectrum::GetAxisFunc(), and maxE.

Referenced by ApplyMaxE(), and ChisqDOF().

00312 {
00313     if (compress) {
00314         TF1 *f = NuMatrixSpectrum::GetAxisFunc();
00315         return f->Eval(maxE);
00316     }
00317     else {
00318         return maxE;
00319     }
00320 }

double NuStatistics::KolmogorovBoth (  ) 

Definition at line 242 of file NuStatistics.cxx.

References ChisqCount(), CombineSignificance(), and KolmogorovShape().

00243 {
00244     double prb2 = ChisqCount();
00245     double prb1 = KolmogorovShape();
00246     double prob = CombineSignificance(prb1, prb2);
00247     
00248     printf(" Kolmo Both     = %g \n",prob);
00249     
00250     return prob;
00251 }

double NuStatistics::KolmogorovShape (  ) 

Definition at line 78 of file NuStatistics.cxx.

References data, MuELoss::e, NuMatrixSpectrum::GetNbinsX(), Msg::kError, MSG, pred, and NuMatrixSpectrum::Spectrum().

Referenced by KolmogorovBoth().

00079 {
00080     if (data->GetNbinsX() < 90) {
00081         cout << "Warning.  This spectrum has probably already been formatted for display." << endl
00082         << "The Kolmogorv test works  best with very small or no binning.  You should" << endl
00083         << "do the test BEFORE rebinning (i.e. BlessedND or BlessedFD)." << endl;
00084     }
00085 
00086     TH1D *h1 = data->Spectrum();
00087     TH1D *h2 = pred->Spectrum();
00088     
00089     // Get Kolmogorov probability
00090     //double prob = h1->KolmogorovTest(h2, "D");
00091     //double dfmax = h1->KolmogorovTest(h2, "M");
00092     
00093     Double_t prob = 0;
00094 
00095     TAxis *axis1 = h1->GetXaxis();
00096     TAxis *axis2 = h2->GetXaxis();
00097     Int_t ncx1   = axis1->GetNbins();
00098     Int_t ncx2   = axis2->GetNbins();
00099     
00100     // Check consistency of dimensions
00101     if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
00102         MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Histograms must be 1-D" << endl;
00103         return 0;
00104     }
00105     
00106     // Check consistency in number of channels
00107     if (ncx1 != ncx2) {
00108         MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Number of channels is different, "
00109         << ncx1 << " and " << ncx2 << endl;
00110         return 0;
00111     }
00112     
00113     // Check consistency in channel edges
00114     Double_t difprec = 1e-5;
00115     Double_t diff1 = TMath::Abs(axis1->GetXmin() - axis2->GetXmin());
00116     Double_t diff2 = TMath::Abs(axis1->GetXmax() - axis2->GetXmax());
00117     if (diff1 > difprec || diff2 > difprec) {
00118         MSG("NuStatistics",Msg::kError) << "KolmogorovTest: histograms with different binning" << endl;
00119         return 0;
00120     }
00121     
00122     Bool_t afunc1 = kFALSE;
00123     Bool_t afunc2 = kFALSE;
00124     Double_t sum1 = 0, sum2 = 0;
00125     Double_t ew1, ew2, w1 = 0, w2 = 0;
00126     Int_t bin;
00127     Int_t ifirst = 1;
00128     Int_t ilast = ncx1;
00129     // integral of all bins (use underflow/overflow if option)
00130     for (bin = ifirst; bin <= ilast; bin++) {
00131         sum1 += h1->GetBinContent(bin);
00132         sum2 += h2->GetBinContent(bin);
00133         ew1   = h1->GetBinError(bin);
00134         ew2   = h2->GetBinError(bin);
00135         w1   += ew1*ew1;
00136         w2   += ew2*ew2;
00137     }
00138     if (sum1 == 0) {
00139         MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Histogram1 " 
00140         << h1->GetName() << " integral is zero" << endl;
00141         return 0;
00142     }
00143     if (sum2 == 0) {
00144         MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Histogram2 " 
00145         << h2->GetName() << " integral is zero" << endl;
00146         return 0;
00147     }
00148     
00149     // calculate the effective entries.  
00150     // the case when errors are zero (w1 == 0 or w2 ==0) are equivalent to 
00151     // compare to a function. In that case the rescaling is done only on sqrt(esum2) or sqrt(esum1) 
00152     Double_t esum1 = 0, esum2 = 0; 
00153     if (w1 > 0) 
00154         esum1 = sum1 * sum1 / w1; 
00155     else 
00156         afunc1 = kTRUE;    // use later for calculating z
00157     
00158     if (w2 > 0) 
00159         esum2 = sum2 * sum2 / w2; 
00160     else 
00161         afunc2 = kTRUE;    // use later for calculating z
00162     
00163     if (afunc2 && afunc1) { 
00164         MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Errors are zero for both histograms" << endl;
00165         return 0;
00166     }
00167     
00168     
00169     Double_t s1 = 1/sum1;
00170     Double_t s2 = 1/sum2;
00171     
00172     // Find largest difference for Kolmogorov Test
00173     Double_t dfmax =0, rsum1 = 0, rsum2 = 0;
00174     
00175     for (bin=ifirst;bin<=ilast;bin++) {
00176         rsum1 += s1*h1->GetBinContent(bin);
00177         rsum2 += s2*h2->GetBinContent(bin);
00178         dfmax = TMath::Max(dfmax,TMath::Abs(rsum1-rsum2));
00179     }
00180     
00181     // Get Kolmogorov probability
00182     Double_t z;
00183     
00184     // case h1 is exact (has zero errors)
00185     if  (afunc1)  
00186         z = dfmax*TMath::Sqrt(esum2);
00187     // case h2 has zero errors
00188     else if (afunc2) 
00189         z = dfmax*TMath::Sqrt(esum1);
00190     else 
00191         // for comparison between two data sets 
00192         z = dfmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
00193     
00194     prob = TMath::KolmogorovProb(z);
00195     
00196     
00197     // debug printout
00198     printf(" Kolmo Shape h1 = %s, sum bin content =%g  effective entries =%g\n",h1->GetName(),sum1,esum1);
00199     printf(" Kolmo Shape h2 = %s, sum bin content =%g  effective entries =%g\n",h2->GetName(),sum2,esum2);
00200     printf(" Kolmo Shape  p = %g, Max Dist = %g\n",prob,dfmax);
00201 
00202     
00203     
00204     return prob;
00205     
00206 }

void NuStatistics::PredSpectrum ( NuMatrixSpectrum spec  ) 
void NuStatistics::SetMaxEnergy ( double  en  )  [inline]

Definition at line 20 of file NuStatistics.h.

References maxE.

00020 {maxE = en;};


Member Data Documentation

Definition at line 40 of file NuStatistics.h.

Referenced by ChisqCount(), DataSpectrum(), KolmogorovShape(), NuStatistics(), and ~NuStatistics().

Definition at line 43 of file NuStatistics.h.

Referenced by ChisqShape(), DataSpectrum(), NuStatistics(), and ~NuStatistics().

int NuStatistics::fitPars [private]

Definition at line 37 of file NuStatistics.h.

Referenced by FitPars().

double NuStatistics::maxE [private]

Definition at line 38 of file NuStatistics.h.

Referenced by GetMaxE(), and SetMaxEnergy().

Definition at line 41 of file NuStatistics.h.

Referenced by ChisqCount(), KolmogorovShape(), NuStatistics(), PredSpectrum(), and ~NuStatistics().

Definition at line 44 of file NuStatistics.h.

Referenced by ChisqDOF(), ChisqShape(), NuStatistics(), PredSpectrum(), and ~NuStatistics().


The documentation for this class was generated from the following files:

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1