ExperimentStatisticsExample.md

post link

plot1.cpp

#include "TROOT.h"
#include "TCanvas.h"
#include "TPaveText.h"
#include "TF1.h"
#include "TRandom3.h"
#include "TH1D.h"
#include "TStyle.h"
#include "TH2D.h"
#include "TFitResultPtr.h"
#include "TLine.h"

void plot1(){
    gStyle->SetOptStat(0);
    gStyle->SetOptTitle(0);
    gStyle->SetOptDate(0);

    Double_t mean = 5.0;
    Double_t stdev = 2.0;
    
    TCanvas * c1 = new TCanvas("c1", "c1", 1200, 800);
    TF1 * f1 = new TF1("gaus", TString::Format("TMath::Gaus(x, %.2f, %.2f)", mean, stdev), -5, 15);
    f1->GetXaxis()->SetTitle("Measured Value");
    f1->GetYaxis()->SetTitle("Relative Frequency");
    f1->SetLineColor(kBlack);
    f1->SetLineWidth(6.0);
    f1->Draw();

    TRandom3 * r = new TRandom3();

    TH1D * h1 = new TH1D("h1", "h1", 150, -5, 15);
    std::cout.precision(3);
    
    for(Int_t i=0; i<5; i++){
        Double_t val = r->Gaus(mean, stdev);
        Double_t height = f1->Eval(val);

        std::cout<<val<<std::endl;

        TLine * line = new TLine(val, 0, val, height);
        line->SetLineColor(kRed);
        line->Draw();

        h1->Fill(val);
    }
    c1->SaveAs("p00_c1.root");
    c1->SaveAs("p00_c1.png");
    
    TCanvas * c2 = new TCanvas("c2", "c2", 1200, 800);
    h1->Draw();
    h1->SetLineWidth(3.0);
    h1->GetXaxis()->SetTitle("Measured Value");
    h1->GetYaxis()->SetTitle("Entries");

    c2->SaveAs("p00_c2.root");
    c2->SaveAs("p00_c2.png");


    TH1D * h2_1 = new TH1D("h2_1", "h2", 40, -15, 25);
    TH2D * h2_2 = new TH2D("h2_2", "h2", 200, 0.5, 200.5, 150, -5, 15);
    TH2D * h2_3 = new TH2D("h2_3", "h2", 200, 0.5, 200.5, 150, -5, 15);

    Double_t means[200];
    Double_t val2[200];

    for(Int_t i=0; i<20000; i++){
        Double_t val = r->Gaus(mean, stdev);
        for(Int_t j=0; j<200; j++){
            if(i==0){
                means[j]=0;
            }
            means[j] += val/Double_t(j+1);
            if(j!=0){
                val2[j]  += TMath::Power(val,2)/Double_t(j+1);
            }

            if(i%(j+1)==0 && i!=0){
                h2_2->Fill(j, means[j]);
                if(j!=0){
                    h2_3->Fill(j, TMath::Sqrt(  (val2[j]-TMath::Power(means[j],2)) * j/(j-1) )    );
                }
                means[j] = 0;
                val2[j] = 0;
            }
        }
        h2_1->Fill(val);
    }

    TCanvas *c3 = new TCanvas("c3", "c3", 1800, 1200);

    c3->Divide(1,2);
    c3->cd(1);
    h2_2->Draw("COLZ");

    c3->cd(2);
    h2_3->Draw("COLZ");

    TCanvas *c4 = new TCanvas("c4", "c4", 1800, 1000);
    c4->Divide(2,1);
    c4->Draw();
    for(int i=0; i<60; i++){
        TH1D * hh1 = h2_2->ProjectionY(TString::Format("h22_py_%0d", i), i, i);
        TH1D * hh2 = h2_3->ProjectionY(TString::Format("h23_py_%0d", i), i, i);

        TText * t1 = new TText(0.2, 0.7, TString::Format("N: %d", i+1));
        t1->SetNDC();
        
        hh1->GetXaxis()->SetTitle("Mean of Sampled Value");
        hh1->GetYaxis()->SetTitle("Entries");
        hh1->SetLineWidth(2.0);

        TF1 * f1 = new TF1(TString::Format("f1_%0d",i), "TMath::Gaus(x, [1], [2])*[0]", -5, 15);

        hh2->GetXaxis()->SetTitle("STDEV of Sampled Value");
        hh2->GetYaxis()->SetTitle("Entries");
        hh2->SetLineWidth(2.0);
        // TText * meantext = new TText(0.8, 0.8, TString::Format("ASDF"));
        TPaveText * pave1 = new TPaveText(0.6, 0.65, 0.85, 0.85, "NDC");
        pave1->AddText(TString::Format("Histogram Entity"));
        pave1->AddText(TString::Format("Entries: %.0f", hh1->GetEntries()));
        pave1->AddText(TString::Format("Integral: %.3f", hh1->Integral("width")));
        pave1->AddText(TString::Format("Mean: %.3f", hh1->GetMean()));
        pave1->AddText(TString::Format("STDEV: %.3f", hh1->GetStdDev()));
        pave1->AddText(TString::Format("STDEV x sqrt(N): %.3f", TMath::Sqrt(i+1)*hh1->GetStdDev()));

        c4->cd(1);
        hh1->Draw();
        t1->Draw();
        pave1->Draw();
        
        f1->SetParameter(0, hh1->Integral("width"));
        f1->SetParameter(1, hh1->GetMean());
        f1->SetParameter(2, hh1->GetStdDev());
        TFitResultPtr ff = hh1->Fit(f1, "S");
        
        // TPaveText * pave1_1 = new TPaveText(0.15, 0.60, 0.40, 0.8, "NDC");
        TPaveText * pave1_1 = new TPaveText(0.6, 0.43, 0.85, 0.60, "NDC");
        pave1_1->AddText(TString::Format("Gaussian Fitting"));
        pave1_1->AddText(TString::Format("Integral: %.2f", f1->Integral(-5, 15)));
        pave1_1->AddText(TString::Format("Mean: %.3f", f1->GetParameter(1)));
        pave1_1->AddText(TString::Format("STDEV: %.3f", f1->GetParameter(2)));
        pave1_1->AddText(TString::Format("Chi2/ndf: %.3f", ff->Chi2()/ff->Ndf()));
        pave1_1->Draw();

        c4->cd(2);
        TPaveText * pave2 = new TPaveText(0.6, 0.65, 0.85, 0.8, "NDC");
        pave2->AddText(TString::Format("N: %d", i+1));
        pave2->AddText(TString::Format("Mean: %.3f", hh2->GetMean()));
        pave2->AddText(TString::Format("STDEV: %.3f", hh2->GetStdDev()));
        pave2->AddText(TString::Format("STDEV x sqrt(N): %.3f", TMath::Sqrt(i+1)*hh2->GetStdDev()));

        hh2->Draw();
        pave2->Draw();

        c4->SaveAs(TString::Format("Projection_%0d.png", i));
    }
}

plot2.cpp

#include "TROOT.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TText.h"

void plot2(){

    gStyle->SetOptTitle(0);

    TF1 * f1_1 = new TF1("f1_1", "TMath::Erf( ((x-0.5)/0.1) )*0.5+0.5", 0, 1);
    TF1 * f1_2 = new TF1("f1_2", "0.5-TMath::Erf( ((x-0.5)/0.1) )*0.5", 0, 1);

    f1_1->SetLineColor(kBlack);
    f1_2->SetLineColor(kPink);
    f1_2->SetLineStyle(10);

    f1_1->GetXaxis()->SetTitle("True Value");
    f1_1->GetYaxis()->SetTitle("Probability of measurement");

    TF1 * f2_1 = new TF1("f2_1", "TMath::Erf( ((x-0.5)/0.01) )*0.5+0.5", 0, 1);
    TF1 * f2_2 = new TF1("f2_2", "0.5-TMath::Erf( ((x-0.5)/0.01) )*0.5", 0, 1);

    f2_1->SetLineColor(kBlack);
    f2_2->SetLineColor(kPink);
    f2_2->SetLineStyle(10);

    f2_1->GetXaxis()->SetTitle("True Value");
    f2_1->GetYaxis()->SetTitle("Probability of measurement");

    TCanvas * c1 = new TCanvas("c1", "c1", 1600, 800);
    c1->Draw();
    c1->Divide(2,1);

    c1->cd(1);

    f1_1->Draw();
    f1_2->Draw("SAME");

    TLegend* legend1 = new TLegend(0.15, 0.375, 0.45, 0.575);
    legend1->AddEntry(f1_1, "1 (Ceiling)");
    legend1->AddEntry(f1_2, "0 (Floor)");
    legend1->SetTextSize(0.044);
    legend1->Draw();
    TText * text1 = new TText(0.6, 0.5, "sigma=0.1");
    text1->Draw("SAME");

    c1->cd(2);

    f2_1->Draw();
    f2_2->Draw("SAME");
    

    TLegend* legend2 = new TLegend(0.15, 0.375, 0.45, 0.575);
    legend2->AddEntry(f2_1, "1 (Ceiling)");
    legend2->AddEntry(f2_2, "0 (Floor)");
    legend2->SetTextSize(0.044);
    // legend2->SetEntrySeparation(-2.);
    // legend2->SetMargin(0.3);
    legend2->Draw();
    TText * text2 = new TText(0.6, 0.5, "sigma=0.01");
    
    text2->Draw("SAME");

    c1->SaveAs("p01_fig1.png");
}