// ----------------------------------------------------------------------------------- // /* ROOT macro for fitting single exponential. This ROOT macro illustrates a fit to a single exponential, and how the results will distribute themselves - warm up problem for the two exponential case. References: Cowan: Chapter 1.5 Bevington: Chapter 11.2 Author: Troels C. Petersen (CERN) Email: Troels.Petersen@cern.ch Date: 24th of September 2012 */ // ----------------------------------------------------------------------------------- // // ----------------------------------------------------------------------------------- // { // ----------------------------------------------------------------------------------- // gROOT->Reset(); // Setting of general plotting style: gStyle->SetCanvasColor(0); gStyle->SetFillColor(0); // Setting what to be shown in statistics box: gStyle->SetOptStat("emr"); gStyle->SetOptFit(1111); // Random numbers from the Mersenne-Twister: TRandom3 r; // r.SetSeed(0); // Initializes with the current time (i.e. random!) r.SetSeed(1); // Initializes with 1, "fixing" random number array. const int Nexp = 1; const int Nevents = 1000; double tau0 = 20.0; // ------------------------------------------------------------------ // // Make Histograms: // ------------------------------------------------------------------ // TH1F* Hist_tau0 = new TH1F("Hist_tau0", "Hist_tau0", 50,10.0, 30.0); TH1F* Hist_etau0 = new TH1F("Hist_etau0", "Hist_etau0", 50, 0.0, 2.0); TH1F* Hist_Prob = new TH1F("Hist_Prob", "Hist_Prob", 20, 0.0, 1.0); const int Npar = 2; double CorrMat_last[Npar][Npar]; TH1F* Hist_times = new TH1F("Hist_times", "Hist_times", 100, 0.0, 100.0); TH1F* Hist_times_last = new TH1F("Hist_times_last", "Hist_times_last", 100, 0.0, 100.0); // ------------------------------------------------------------------ // // Loop over generating and fitting data: // ------------------------------------------------------------------ // for (int iexp=0; iexp < Nexp; iexp++) { // Generate data: // -------------- for (int i=0; i < Nevents; i++) { double t = r.Exp(tau0); Hist_times->Fill(t); } // Fit data: // --------- TF1* fit = new TF1("fit", "[0]*exp(-x/[1])", 0.0, 100.0); fit->SetParameters(1000.0, 20.0); Double_t CovMat[Npar][Npar]; TMinuit *gMinuit = new TMinuit(Npar); Hist_times->Fit("fit","rlqn"); gMinuit->mnemat(&CovMat[0][0],Npar); // Get correlation matrix (and save a copy): for (int i=0;i 0.0 && CovMat[j][j] > 0) { CorrMat_last[i][j] = CovMat[i][j] / (sqrt(CovMat[i][i])*sqrt(CovMat[j][j])); } else { printf(" ERROR: Diagonal covariance entries 0 or negative! %9.4f %9.4f \n", CovMat[i][i], CovMat[j][j]); break; } } } double tau0_fit = fit->GetParameter(1); double etau0_fit = fit->GetParError(1); double Prob = fit->GetProb(); printf(" %4d Fit result: %6.3f+-%5.3f Prob=%6.4f \n", iexp, tau0_fit, etau0_fit, Prob); // Fill histograms: Hist_tau0->Fill(tau0_fit); Hist_etau0->Fill(etau0_fit); Hist_Prob->Fill(Prob); if (Nexp > 1) Hist_times->Reset("ICE"); else Hist_times_last = Hist_times; } // ------------------------------------------------------------------ // // Fit the data and plot the result: // ------------------------------------------------------------------ // if (Nexp == 1) { // Make a white canvas and draw the example fit in: c0 = new TCanvas("c0","",200,20,600,450); c0->SetFillColor(0); Hist_times_last->SetLineWidth(2); Hist_times_last->SetLineColor(4); TF1* fit = new TF1("fit", "[0]*exp(-x/[1])", 0.0, 100.0); fit->SetParameters(500.0, 20.0); Hist_times_last->Fit("fit","r"); Hist_times_last->Draw("e"); // Plot correlation matrix in plot: TPad* pad = new TPad("pad","testpad", 0.21, 0.49, 0.61, 0.99, 0, 1, 1); pad->Draw(); pad->cd(); TH2F* Corr = new TH2F("Corr","",Npar,-0.5,double(Npar)-0.5,Npar,-0.5,double(Npar)-0.5); printf(" Correlation matrix for fitting constants: \n"); for (int i=0;iSetBinContent(i+1,j+1, int(100.0 * CorrMat_last[i][j])); printf(" %6.3f", CorrMat_last[i][j]); } printf("\n"); } Corr->Draw("colz text"); } // ------------------------------------------------------------------ // // Show the distribution of fitting results: // ------------------------------------------------------------------ // if (Nexp > 1) { // Make a white canvas and draw the example fit in: c1 = new TCanvas("c1","",250,70,600,450); c1->SetFillColor(0); c1->Divide(2,2); c1->cd(1); Hist_tau0->SetLineWidth(2); Hist_tau0->SetLineColor(4); Hist_tau0->Draw(); c1->cd(2); Hist_etau0->SetLineWidth(2); Hist_etau0->SetLineColor(4); Hist_etau0->Draw(); c1->cd(3); Hist_Prob->SetLineWidth(2); Hist_Prob->SetLineColor(4); Hist_Prob->Draw(); } } //---------------------------------------------------------------------------------- /* Questions: ---------- 1) Do a quick evaluation of the fit, and make sure you consider the following aspects: - Visual inspection of the fit - Chi2 / Ndof and probability - Fit result and errors compared to input Did you wonder, why the normalization did not really match the number of events? Did you think about the posibility that some events fell below or above the range? Did you consider what uncertainties ROOT assigned to each bin? If not, then go directly to jail! Even if you pass start, you will not cash 4000 Kr. 2) Consider the correlation matrix. How large are the correlations, and are they something to worry about? Now, try to add the lifetime in the normalization of the fitting PDF, and repeat the fit. Did this change things? And if so, do you think you understand why? 3) Now run the experiment many times and see if the results and errors are as you would expect them, both with and without the lifetime in the normalization. */ //----------------------------------------------------------------------------------