NuStatistics.cxx

Go to the documentation of this file.

#include <iostream>

#include "TF1.h"
#include "TH1.h"
#include "TMath.h"

#include "NtupleUtils/NuStatistics.h"
#include "NtupleUtils/NuMatrixSpectrum.h"
#include "MessageService/MsgService.h"

using namespace std;

CVSID("");

//____________________________________________________________________72
NuStatistics::NuStatistics() 
: fitPars(0), maxE(50.)
{
    data = 0;
    data_rb = 0;
    pred = 0;
    pred_rb = 0;
}


//____________________________________________________________________72
NuStatistics::~NuStatistics() 
{
    if (data) delete data;
    if (data_rb) delete data_rb;
    if (pred) delete pred;
    if (pred_rb) delete pred_rb;  
}


//____________________________________________________________________72
void NuStatistics::DataSpectrum(NuMatrixSpectrum *spec)
{
    if (data) delete data;
    data = new NuMatrixSpectrum(*spec);
    data->RemoveErrors();
    data->Spectrum()->SetName("Data");
    ApplyMaxE(data->Spectrum());
    
    if (data_rb) delete data_rb;
    data_rb = new NuMatrixSpectrum(*spec);
    data_rb->BlessedFD();
    data_rb->RemoveErrors();
    data_rb->Spectrum()->SetName("DataDisplay");
    ApplyMaxE(data_rb->Spectrum(), true);
}

//____________________________________________________________________72
void NuStatistics::PredSpectrum(NuMatrixSpectrum *spec)
{
    if (pred) delete pred;
    pred = new NuMatrixSpectrum(*spec);
    pred->RemoveErrors();
    pred->GenerateSystErrors();
    pred->GenerateErrors();
    pred->Spectrum()->SetName("Pred");
    ApplyMaxE(pred->Spectrum());  

    if (pred_rb) delete pred_rb;
    pred_rb = new NuMatrixSpectrum(*spec);
    pred_rb->RebinToScheme(NuBinningScheme::kDisplayFD);
    pred_rb->RemoveErrors();
    pred_rb->GenerateSystErrors();
    pred_rb->GenerateErrors();
    pred_rb->BinWidthNormalize();
    pred_rb->NuBarCompressX();
    pred_rb->Spectrum()->SetName("PredDisplay");
    ApplyMaxE(pred_rb->Spectrum(), true);
}


//____________________________________________________________________72
double NuStatistics::KolmogorovShape()
{
    if (data->GetNbinsX() < 90) {
        cout << "Warning.  This spectrum has probably already been formatted for display." << endl
        << "The Kolmogorv test works  best with very small or no binning.  You should" << endl
        << "do the test BEFORE rebinning (i.e. BlessedND or BlessedFD)." << endl;
    }

    TH1D *h1 = data->Spectrum();
    TH1D *h2 = pred->Spectrum();
    
    // Get Kolmogorov probability
    //double prob = h1->KolmogorovTest(h2, "D");
    //double dfmax = h1->KolmogorovTest(h2, "M");
    
    Double_t prob = 0;

    TAxis *axis1 = h1->GetXaxis();
    TAxis *axis2 = h2->GetXaxis();
    Int_t ncx1   = axis1->GetNbins();
    Int_t ncx2   = axis2->GetNbins();
    
    // Check consistency of dimensions
    if (h1->GetDimension() != 1 || h2->GetDimension() != 1) {
        MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Histograms must be 1-D" << endl;
        return 0;
    }
    
    // Check consistency in number of channels
    if (ncx1 != ncx2) {
        MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Number of channels is different, "
        << ncx1 << " and " << ncx2 << endl;
        return 0;
    }
    
    // Check consistency in channel edges
    Double_t difprec = 1e-5;
    Double_t diff1 = TMath::Abs(axis1->GetXmin() - axis2->GetXmin());
    Double_t diff2 = TMath::Abs(axis1->GetXmax() - axis2->GetXmax());
    if (diff1 > difprec || diff2 > difprec) {
        MSG("NuStatistics",Msg::kError) << "KolmogorovTest: histograms with different binning" << endl;
        return 0;
    }
    
    Bool_t afunc1 = kFALSE;
    Bool_t afunc2 = kFALSE;
    Double_t sum1 = 0, sum2 = 0;
    Double_t ew1, ew2, w1 = 0, w2 = 0;
    Int_t bin;
    Int_t ifirst = 1;
    Int_t ilast = ncx1;
    // integral of all bins (use underflow/overflow if option)
    for (bin = ifirst; bin <= ilast; bin++) {
        sum1 += h1->GetBinContent(bin);
        sum2 += h2->GetBinContent(bin);
        ew1   = h1->GetBinError(bin);
        ew2   = h2->GetBinError(bin);
        w1   += ew1*ew1;
        w2   += ew2*ew2;
    }
    if (sum1 == 0) {
        MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Histogram1 " 
        << h1->GetName() << " integral is zero" << endl;
        return 0;
    }
    if (sum2 == 0) {
        MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Histogram2 " 
        << h2->GetName() << " integral is zero" << endl;
        return 0;
    }
    
    // calculate the effective entries.  
    // the case when errors are zero (w1 == 0 or w2 ==0) are equivalent to 
    // compare to a function. In that case the rescaling is done only on sqrt(esum2) or sqrt(esum1) 
    Double_t esum1 = 0, esum2 = 0; 
    if (w1 > 0) 
        esum1 = sum1 * sum1 / w1; 
    else 
        afunc1 = kTRUE;    // use later for calculating z
    
    if (w2 > 0) 
        esum2 = sum2 * sum2 / w2; 
    else 
        afunc2 = kTRUE;    // use later for calculating z
    
    if (afunc2 && afunc1) { 
        MSG("NuStatistics",Msg::kError) << "KolmogorovTest: Errors are zero for both histograms" << endl;
        return 0;
    }
    
    
    Double_t s1 = 1/sum1;
    Double_t s2 = 1/sum2;
    
    // Find largest difference for Kolmogorov Test
    Double_t dfmax =0, rsum1 = 0, rsum2 = 0;
    
    for (bin=ifirst;bin<=ilast;bin++) {
        rsum1 += s1*h1->GetBinContent(bin);
        rsum2 += s2*h2->GetBinContent(bin);
        dfmax = TMath::Max(dfmax,TMath::Abs(rsum1-rsum2));
    }
    
    // Get Kolmogorov probability
    Double_t z;
    
    // case h1 is exact (has zero errors)
    if  (afunc1)  
        z = dfmax*TMath::Sqrt(esum2);
    // case h2 has zero errors
    else if (afunc2) 
        z = dfmax*TMath::Sqrt(esum1);
    else 
        // for comparison between two data sets 
        z = dfmax*TMath::Sqrt(esum1*esum2/(esum1+esum2));
    
    prob = TMath::KolmogorovProb(z);
    
    
    // debug printout
    printf(" Kolmo Shape h1 = %s, sum bin content =%g  effective entries =%g\n",h1->GetName(),sum1,esum1);
    printf(" Kolmo Shape h2 = %s, sum bin content =%g  effective entries =%g\n",h2->GetName(),sum2,esum2);
    printf(" Kolmo Shape  p = %g, Max Dist = %g\n",prob,dfmax);

    
    
    return prob;
    
}


//____________________________________________________________________72
double NuStatistics::ChisqCount()
{
    TH1D *h1 = data->Spectrum();
    TH1D *h2 = pred->Spectrum();
    
    double sum1 = 0, sum2 = 0;
    double ew1, ew2, w1 = 0, w2 = 0;
    

    for (int bin = 1; bin <= h1->GetNbinsX(); bin++) {
        sum1 += h1->GetBinContent(bin);
        sum2 += h2->GetBinContent(bin);
        ew1   = h1->GetBinError(bin);
        ew2   = h2->GetBinError(bin);
        w1   += ew1*ew1;
        w2   += ew2*ew2;
    }
    
    // My simple chisquared method
    double delta = sum1 - sum2;
    double sigmasq = w1 + w2;
    double chi2 = delta*delta / sigmasq;
    double prob = TMath::Prob(chi2, 1);
    
    printf(" Chisq Norm  h1 = %s, sum bin content =%g  sumsq error =%g\n",h1->GetName(),sum1,w1);
    printf(" Chisq Norm  h2 = %s, sum bin content =%g  sumsq error =%g\n",h2->GetName(),sum2,w2);
    printf(" Chisq Norm     = %g\n",prob);
    
    return prob;
}

//____________________________________________________________________72
double NuStatistics::KolmogorovBoth()
{
    double prb2 = ChisqCount();
    double prb1 = KolmogorovShape();
    double prob = CombineSignificance(prb1, prb2);
    
    printf(" Kolmo Both     = %g \n",prob);
    
    return prob;
}

//____________________________________________________________________72
double NuStatistics::ChisqShape()
{
    TH1D *h1 = data_rb->Spectrum();
    TH1D *h2 = pred_rb->Spectrum();
    
    double chisq = 0;
    double da, mc, eda, emc;
    double delta, sigmasq;

    for (int i=1; i<= h1->GetNbinsX(); i++){
        da = h1->GetBinContent(i);
        mc = h2->GetBinContent(i);
        
        eda = h1->GetBinError(i);
        emc = h2->GetBinError(i);
        
        delta = da - mc;
        sigmasq = eda*eda + emc*emc;
        
        if (sigmasq <= 0) continue;
        
        chisq += delta*delta/sigmasq;
    }
    
    printf(" Chisq Shape  chi2 = %g\n",chisq);
    
    
    return chisq;
}


//____________________________________________________________________72
int NuStatistics::ChisqDOF()
{
    TH1D *h1 = pred_rb->Spectrum();
    
    int dof = -1 * FitPars();
    for (int i = 1; i <= h1->GetNbinsX(); i++){
        if (h1->GetBinLowEdge(i) >= GetMaxE(true) ) break;
        ++dof;
    }
    
    return dof;
}


//____________________________________________________________________72
double NuStatistics::ChisqProb()
{
    double chi2 = ChisqShape();
    int dof = ChisqDOF();
    return TMath::Prob(chi2, dof);
}



//____________________________________________________________________72
double NuStatistics::GetMaxE(bool compress)
{
    if (compress) {
        TF1 *f = NuMatrixSpectrum::GetAxisFunc();
        return f->Eval(maxE);
    }
    else {
        return maxE;
    }
}


//____________________________________________________________________72
void NuStatistics::ApplyMaxE(TH1* h1, bool compress)
{
    double max = GetMaxE(compress);
    int ifirst = 0; 
    int ilast = h1->GetNbinsX()+1;
    for (int bin = ifirst; bin <= ilast; bin++) {
        if (h1->GetBinLowEdge(bin) >= max) {
            // Make sure excluded bins do not contribute
            h1->SetBinContent(bin, 0);
            h1->SetBinError(bin, 0);
        }
    }
}


//____________________________________________________________________72
double NuStatistics::CombineSignificance(double prb1, double prb2)
{
    // see Eadie et al., section 11.6.2
    
    double prob;
    if (prb1 > 0 && prb2 > 0) prob = prb1*prb2*(1-TMath::Log(prb1*prb2));
    else                      prob = 0;    
    
    return prob;
}

---