// ----------------------------------------------------------------------------------- // /* ROOT macro for fitting two exponentials and considering the correlations involved. This ROOT macro illustrates the importance of correlations in fits, using two exponentials as an example. The pseudo experiments fitted are generated along, and the user has full control of what is fitted. 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 Npar = 4; const int Nexp = 1; const int Nevents = 1000; double tau0 = 2.0; double tau1 = 30.0; double frac0 = 0.5; int n0 = int(double(Nevents) * frac0); int n1 = int(double(Nevents) * (1.0-frac0)); // ------------------------------------------------------------------ // // Make Histograms: // ------------------------------------------------------------------ // TH1F* Hist_tau0 = new TH1F("Hist_tau0", "Hist_tau0", 50, 0.0, 5.0); TH1F* Hist_etau0 = new TH1F("Hist_etau0", "Hist_etau0", 50, 0.0, 1.0); TH1F* Hist_tau1 = new TH1F("Hist_tau1", "Hist_tau1", 50, 0.0, 50.0); TH1F* Hist_etau1 = new TH1F("Hist_etau1", "Hist_etau1", 50, 0.0, 10.0); TH1F* Hist_Prob = new TH1F("Hist_Prob", "Hist_Prob", 20, 0.0, 1.0); TH2F* Hist_tau0tau1 = new TH2F("Hist_tau0tau1", "Hist_tau0tau1", 50, 1.0, 3.0, 50, 20.0, 40.0); 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 < n0; i++) { double t = r.Exp(tau0); Hist_times->Fill(t); } for (int i=0; i < n1; i++) { double t = r.Exp(tau1); Hist_times->Fill(t); } // Fit data: // --------- // Original function: TF1* fit = new TF1("fit", "[0]*exp(-x/[1]) + [2]*exp(-x/[3])", 0.0, 100.0); // fit->SetParameters(double(n0)/10.0, tau0, double(n1)/10.0, tau1); // Alternative function: // TF1* fit = new TF1("fit", "[0] * ([2]/[1]*exp(-x/[1]) + (1.0-[2])/[3]*exp(-x/[3]))", 0.0, 100.0); // fit->SetParameters(double(Nevents), tau0, frac0, tau1); 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 tau1_fit = fit->GetParameter(3); double etau0_fit = fit->GetParError(1); double etau1_fit = fit->GetParError(3); double Prob = fit->GetProb(); printf(" %4d Fit result: %6.3f+-%5.3f %6.3f+-%5.3f Prob=%6.4f \n", iexp, tau0_fit, etau0_fit, tau1_fit, etau1_fit, Prob); // Fill histograms: Hist_tau0->Fill(tau0_fit); Hist_etau0->Fill(etau0_fit); Hist_tau1->Fill(tau1_fit); Hist_etau1->Fill(etau1_fit); Hist_tau0tau1->Fill(tau0_fit,tau1_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); // c0->SetLogy(); // Make the y-axis logarithmic. Hist_times_last->SetLineWidth(2); Hist_times_last->SetLineColor(4); // Original function: TF1* fit = new TF1("fit", "[0]*exp(-x/[1]) + [2]*exp(-x/[3])", 0.0, 100.0); // fit->SetParameters(double(n0)/10.0, tau0, double(n1)/10.0, tau1); // Alternative function: // TF1* fit = new TF1("fit", "[0] * ([2]/[1]*exp(-x/[1]) + (1.0-[2])/[3]*exp(-x/[3]))", 0.0, 100.0); // fit->SetParameters(double(Nevents), tau0, frac0, tau1); 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.46, 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: 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_tau1->SetLineWidth(2); Hist_tau1->SetLineColor(4); Hist_tau1->Draw(); c1->cd(4); Hist_etau1->SetLineWidth(2); Hist_etau1->SetLineColor(4); Hist_etau1->Draw(); // Make a white canvas: c2 = new TCanvas("c2","",300,120,600,450); c2->SetFillColor(0); Hist_tau0tau1->SetMarkerStyle(20); Hist_tau0tau1->SetMarkerSize(0.8); Hist_tau0tau1->SetMarkerColor(4); Hist_tau0tau1->Draw("box"); printf(" The linear correlation coefficient between tau0 and tau1 is: %5.2f \n", Hist_tau0tau1->GetCorrelationFactor()); } } //---------------------------------------------------------------------------------- /* Questions: ---------- 1) Run the macro (after reading through it, of course), and acknowledge that a fit with several correlated parameters requires (good) initial values!!! Note that the fit is defined twice, both in line 91 and line 154, as the first is the general fit (possibly used many times), while the second is for showing one case. After uncommenting to get initial values, see if the fit behaves like you expected. Are the fit results satisfactory and is everything behaving like it should? Try to plot the fit on a log-y plot. Is it now easier to see, that there are two different exponentials involved? 2) Consider the correlations. Are they negligible? And what is correlated with what? Can you think of a way of deminishing these correlations? Try the suggested but commented out fit (with other obvious modifications!!!), and see how well it performs. Compare result and correlations with before. After trying these, go back to the original fitting functions... 3) Now try to run the event generation and fit many times. A new type of plot comes up, which shows the fitted lifetimes in a 2D histogram. Any correlation to be seen? Do they match what the single fit suggested? Also, note how the results and their errors distribute themselves, as you will be asked to compare these with other situations in question 5. 4) The two initial lifetimes (2 and 30) are very different. What would you expect happens, when they are not too different? Think about this first, and then try the values 5 and 10 instead, going back to simply running one fit/experiment. How large are the correlations now? And do you think that this makes the fit reliable? NOTE: When you change the lifetimes, you should also change the histogram ranges! 5) Try to run many experiments with these more correlated lifetimes. First, consider the result of the fitting values and their errors by themselves. Any notable changes from question 3? Which and are they for the better? And what do you see in the 2D plot? What is the linear correlation coefficient, and would you say that the correlation is linear? Is this generally a desirable fitting situation? 6) Now repeat question 5, but using the less correlated functional form. Does this improve the results? Advanced questions: ------------------- 1) Try to make a 2D fit to the 2D distribution of (tau0,tau1) as it came out in question 3. What function do you use? Is it linear? */ //----------------------------------------------------------------------------------