#include "TH1.h" #include "TMath.h" #include "TCanvas.h" #include "TStyle.h" #include "TF1.h" #include "TRandom3.h" #include #include #include // Divide and add 1D Histograms void format_h(TH1D* h, int linecolor){ h->SetLineWidth(3); h->SetLineColor(linecolor); } Double_t crystalballfunction(Double_t *x, Double_t *par) { double A = ( (par[1]/TMath::Power(TMath::Abs(par[0]),par[1]) ) )*TMath::Exp(-(TMath::Abs(par[0])*TMath::Abs(par[0]))/2); double B = par[1]/TMath::Abs(par[0]) - TMath::Abs(par[0]); double C = ( par[1]/TMath::Abs(par[0]) )*( 1/(par[1]-1) )*TMath::Exp(-(par[0]*par[0])/2); double D = TMath::Sqrt(TMath::Pi()/2)*(1+ (TMath::Erf(par[0]/TMath::Sqrt(2) ) ) ); double N = 1/( par[3]*(C+D) ); Double_t fitval = 0.; if( (x[0]-par[2])/par[3] > -par[0]) fitval = N*TMath::Exp(-(0.5*( (x[0]-par[2])*(x[0]-par[2]) )/(par[3]*par[3]) ) ); else fitval = A*TMath::Power( B- ( (x[0]-par[2])/par[3] ), -par[1] ); return fitval; } void check_statistics() { TH1D* h1=new TH1D("Log(M_d) distribution first dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 5., // Lower X Boundary 10.); // Upper X Boundary TH1D* h2=new TH1D("Log(M_d) distribution second dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 0., // Lower X Boundary 10.); // Upper X Boundary TH1D* h3=new TH1D("Log(M_d) distribution third dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 5., // Lower X Boundary 10); // Upper X Boundary TH1D* h4=new TH1D("Log(M_d) distribution fourth dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 5., // Lower X Boundary 10); // Upper X Boundary TH1D* h5=new TH1D("Log(M_d) distribution fifth dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 5., // Lower X Boundary 10.); // Upper X Boundary TH1D* h6=new TH1D("Log(M_d) distribution sixth dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 5., // Lower X Boundary 10.); // Upper X Boundary TH1D* h7=new TH1D("Log(M_d) distribution seventh dataset", "Log(M_d);disk masses;# occurencies", 50, // Number of Bins 5., // Lower X Boundary 10.); // Upper X Boundary ifstream in1; ifstream in2; ifstream in3; ifstream in4; ifstream in5; ifstream in6; ifstream in7; Double_t x1; Double_t x2; Double_t x3; Double_t x4; Double_t x5; Double_t x6; Double_t x7; in1.open("diskmasses1.txt"); in2.open("diskmasses2.txt"); in3.open("diskmasses3.txt"); in4.open("diskmasses4.txt"); in5.open("diskmasses5.txt"); in6.open("diskmasses6.txt"); in7.open("diskmasses7.txt"); /*while (1) { in >> x; if (!in.good()) break; h1->Fill(x); } h1->Draw();*/ /*const float mean_count=3.6; TRandom3 rndgen; // simulate the measurements for (int imeas=0;imeas<400;imeas++) cnt_r_h->Fill(rndgen.Poisson(mean_count)); TCanvas* c= new TCanvas(); cnt_r_h->Draw(); TCanvas* c_norm= new TCanvas(); cnt_r_h->DrawNormalized(); // Print summary cout << "Moments of Distribution:\n" << " - Mean = " << cnt_r_h->GetMean() << " +- " << cnt_r_h->GetMeanError() << "\n" << " - RMS = " << cnt_r_h->GetRMS() << " +- " << cnt_r_h->GetRMSError() << "\n" << " - Skewness = " << cnt_r_h->GetSkewness() << "\n" << " - Kurtosis = " << cnt_r_h->GetKurtosis() << "\n";*/ /*while (1) { in2 >> x2; if (!in2.good()) break; h2->Fill(x2); } h2->Draw(); while (1) { in3 >> x3; if (!in3.good()) break; h3->Fill(x3); } h3->Draw();*/ for(Int_t i=0; i<25;i++) { in1 >> x1; h1->Fill(x1); } for(Int_t i=0; i<24;i++) { in2 >> x2; h2->Fill(x2); } for(Int_t i=0; i<53;i++) { in3 >> x3; h3->Fill(x3); } for(Int_t i=0; i<70;i++) { in4 >> x4; h4->Fill(x4); } for(Int_t i=0; i<47;i++) { in5 >> x5; h5->Fill(x5); } for(Int_t i=0; i<49;i++) { in6 >> x6; h6->Fill(x6); } for(Int_t i=0; i<12;i++) { in7 >> x7; h7->Fill(x7); } // Format Histograms /*TH1D* histos[3]={h1,h2,h3}; for (int i=0;i<3;++i){ histos[i]->Sumw2(); // *Very* Important format_h(histos[i],i+1); }*/ //TCanvas* c_sum= new TCanvas(); format_h(h1,1); format_h(h2,1); format_h(h3,1); format_h(h4,1); format_h(h5,1); format_h(h6,1); format_h(h7,1); TCanvas* c1= new TCanvas(); gStyle->SetOptStat(110000111); gStyle->SetOptFit(111); //TF1 *fit1 = (TF1*)h1->GetFunction("gaus"); // TF1* fit1 = new TF1("crystal",crystalballfunction,5,10,4); //gau->SetParameters(25,4.3); h1->Fit("crystal","EM"); h1->Draw(); fit1->Draw("same"); TCanvas* c2= new TCanvas(); gStyle->SetOptFit(111); //h2->Fit("gaus"); //TF1 *fit2 = (TF1*)h2->GetFunction("gaus"); TF1* fit2 = new TF1("fit2","gaus",5,10); h2->Draw(); h2->Fit("fit2","EM"); fit2->Draw("same"); //fit2->Draw("same"); TCanvas* c3= new TCanvas(); gStyle->SetOptFit(111); TF1* fit3 = new TF1("fit3","gaus",5,10); h3->Draw(); h3->Fit("fit3","EM"); //h3->Fit("gaus"); //TF1 *fit3 = (TF1*)h3->GetFunction("gaus"); fit3->Draw("same"); TCanvas* c4= new TCanvas(); gStyle->SetOptFit(111); TF1* fit4 = new TF1("fit4","gaus",5,10); h4->Draw(); h4->Fit("fit4","EM"); //h4->Fit("gaus"); //TF1 *fit4 = (TF1*)h4->GetFunction("gaus"); fit4->Draw("same"); TCanvas* c5= new TCanvas(); gStyle->SetOptFit(111); TF1* fit5 = new TF1("fit5","gaus",5,10); h5->Draw(); h5->Fit("fit5","EM"); //h5->Fit("gaus"); //TF1 *fit5 = (TF1*)h5->GetFunction("gaus"); //h5->Draw("hist"); fit5->Draw("same"); TCanvas* c6= new TCanvas(); gStyle->SetOptFit(111); TF1* fit6 = new TF1("fit6","gaus",5,10); h6->Draw(); h6->Fit("fit6","EM"); fit6->Draw("same"); //h6->Fit("gaus"); //TF1 *fit6 = (TF1*)h6->GetFunction("gaus"); //h6->Draw("hist"); //fit6->Draw("same"); TCanvas* c7= new TCanvas(); gStyle->SetOptFit(111); TF1* fit7 = new TF1("fit7","gaus",5,10); h7->Draw(); h7->Fit("fit7","EM"); fit7->Draw("same"); //h7->Fit("gaus"); //TF1 *fit7 = (TF1*)h7->GetFunction("gaus"); //h7->Draw("hist"); //fit7->Draw("same"); in1.close(); in2.close(); in3.close(); in4.close(); in5.close(); in6.close(); in7.close(); }