commit f3d3ecb7dbc25d5a6092f701f8777b10141c7286 Author: Ryan@SOLARIS_testStation Date: Fri Aug 16 16:24:21 2024 -0400 init. commit diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..a608475 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.root \ No newline at end of file diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json new file mode 100644 index 0000000..cca8d1b --- /dev/null +++ b/.vscode/c_cpp_properties.json @@ -0,0 +1,17 @@ +{ + "configurations": [ + { + "name": "Linux", + "includePath": [ + "${workspaceFolder}/**", + "/opt/root/**" + ], + "defines": [], + "compilerPath": "/usr/bin/gcc", + "cStandard": "c17", + "cppStandard": "gnu++17", + "intelliSenseMode": "linux-gcc-x64" + } + ], + "version": 4 +} \ No newline at end of file diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..4798fc2 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,13 @@ +{ + "files.associations": { + "AutoFit.C": "cpp", + "peaks.C": "cpp", + "array": "cpp", + "deque": "cpp", + "list": "cpp", + "vector": "cpp", + "string_view": "cpp", + "span": "cpp", + "SaveTH1IntoText.C": "cpp" + } +} \ No newline at end of file diff --git a/AutoFit.C b/AutoFit.C new file mode 100644 index 0000000..3624e5e --- /dev/null +++ b/AutoFit.C @@ -0,0 +1,2866 @@ +/*************************************************** + * This is a root macro for Auto fitting + * + * Created by Tsz Leung (Ryan) TANG, around 2019. + * updated 01-04-2022 + * + * contact goluckyryan@gmail.com + ***************************************************/ + +#ifndef AutoFit_C +#define AutoFit_C + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +// namespace AutoFit{ + +//Global fit paramaters +std::vector BestFitMean; +std::vector BestFitCount; +std::vector BestFitSigma; + +TString recentFitMethod; + +void ShowFitMethod(){ + printf("----------------------- Method of Fitting ---------------\n"); + printf("---------------------------------------------------------\n"); + printf(" fitAuto() - estimate BG, find peak, and fit n-Gauss \n"); + printf(" fitGaussPol() - fit 1 Gauss + pol-n BG\n"); + printf(" fit2GaussP1() - fit 2 Gauss + pol-1 BG \n"); + printf(" fitGF3Pol() - fit GF3 + pol-n BG \n"); + //printf(" fitNGF3() - fit n-GF3, estimated BG \n"); + printf(" fitNGauss() - fit n-Gauss, estimated BG, need input\n"); + printf(" fitNGaussSub() - fit estimated BG with Pol, subtract, fit n-Guass\n"); + printf(" fitNGaussPol() - fit n-Gauss + pol-n BG \n"); + printf(" fitNGaussPolSub() - subtract Pol-n BG, fit n-Gauss \n"); + printf("\n"); + printf("------- Mouse click Fit : \n"); + printf(" clickFitNGaussPol() - fit n-Gauss + pol-n BG \n"); + printf(" clickFitNGaussPolSub() - Fit Pol-n BG, subtract, fit n-Gauss\n"); + printf("\n"); + printf("------- Utility : \n"); + printf(" SaveFitPara() - Save the inital/Best fit parameters.\n"); + printf(" ShowFitMerhod() - Show this menual.\n"); + printf("---------------------------------------------------------\n"); +} + +void AutoFit(){ + ShowFitMethod(); +} + +std::vector SplitStrAF(std::string tempLine, std::string splitter, int shift = 0){ + + std::vector output; + + size_t pos; + do{ + pos = tempLine.find(splitter); // fine splitter + if( pos == 0 ){ //check if it is splitter again + tempLine = tempLine.substr(pos+1); + continue; + } + + std::string secStr; + if( pos == std::string::npos ){ + secStr = tempLine; + }else{ + secStr = tempLine.substr(0, pos+shift); + tempLine = tempLine.substr(pos+shift); + } + + //check if secStr is begin with space + if( secStr.substr(0, 1) == " "){ + secStr = secStr.substr(1); + } + + output.push_back(secStr); + //printf(" |%s---\n", secStr.c_str()); + + }while(pos != std::string::npos ); + + return output; +} + +TColor RGBWheel(double ang){ + + ang = ang * TMath::DegToRad(); + + double r = std::max(0., (1+2*cos(ang))/3.); + double g = std::max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.); + double b = std::max(0., (1 - cos(ang) - sqrt(3)* sin(ang))/3.); + + TColor col(r,g,b); + + return col; + +} + +int nPeaks = 16; +Double_t nGauss(Double_t *x, Double_t *par) { + Double_t result = 0; + for (Int_t p=0;pGetNpar(); + int count = totPar/numParPerPeak; + printf("ID "); + for( int i = 0; i < numParPerPeak; i++){ + printf("par%d ", i); + } + printf("\n"); + + for( int i = 0; i < count ; i++){ + printf("%3d ", i); + for( int j = 0; j < numParPerPeak; j++){ + int parID = numParPerPeak * i + j; + printf("%.3f(%.3f) ", fit->GetParameter(parID), fit->GetParError(parID)); + } + printf("\n"); + } +} + +void GoodnessofFit(TH1F * hist, TF1 * fit){ + + int nBin = hist->GetNbinsX(); + int effBin = 0; + double mean = 0; + double ysq = 0; + double SSR = 0; + double chisq = 0; //with estimated error be sqrt(y) + double Xsq = 0; // for Pearson's chi-sq test + for( int i = 1; i <= nBin; i++){ + + double e = hist->GetBinError(i); + if( e > 0 ) { + effBin ++; + double y = hist->GetBinContent(i); + double x = hist->GetBinCenter(i); + double ybar = fit->Eval(x); + ysq += y*y; + mean += y; + SSR += (y - ybar)*(y-ybar); + chisq += (y - ybar)*(y-ybar)/e/e; + + + if( ybar > e ) { + Xsq += (y - ybar)*(y-ybar)/ybar; + }else{ + Xsq += (y - ybar)*(y-ybar)/e; + } + //printf(" %d | x : %f, y : %f, ybar : %f , X-sq : %f\n", i, x, y, ybar, Xsq); + } + } + mean = mean / nBin; + double SSTotal = ysq + mean*mean; + + int npar = fit->GetNpar(); + int ndf = effBin - npar; + printf("#################################################\n"); + printf("## Goodness of Fit. ##\n"); + printf("#################################################\n"); + printf(" eff. Bin(%d) - numPar(%d) = ndf = %d \n", effBin, npar, ndf); + + printf("============== Regression analysis\n"); + printf("----------------- R-sq \n"); + printf(" SSTotal = %f \n", SSTotal); + printf(" SSR = %f \n", SSR); + printf(" MSR = %f <-- is it similar to sample variance?\n", SSR/ndf); + double Rsq = 1 - SSR/SSTotal; + printf(" R-sq = %f \n", Rsq ); + + printf("----------------- Chi-sq \n"); + printf(" Chi-sq = %f \n", chisq); + printf(" rd. Chi-sq = %f \n", chisq/ndf); + printf("ROOT Chi-Sq = %f , NDF = %d \n", fit->GetChisquare(), fit->GetNDF()); + //================ chi-sq test + printf("============== Hypothesis testing\n"); + printf(" Null Hypothesis : the fitting model is truth. \n"); + printf(" * p-value = prob. that Null Hypo. is truth. \n"); + printf(" * the Pearson's test in here only for background free data \n"); + printf(" Pearson's X-sq = %.2f \n", Xsq); + double p = 1 - TMath::Prob(Xsq, ndf); + printf(" Pearson's p-value = %.2f %s 0.05 | %s\n", p, p < 0.05 ? "<": ">", p < 0.05 ? "reject" : "cannot reject"); + double pchi = 1 - TMath::Prob(chisq, ndf); + printf(" Chi-sq p-value = %.2f %s 0.05 | %s\n", pchi, pchi < 0.05 ? "<": ">", pchi < 0.05 ? "reject" : "cannot reject"); + double pRoot = 1- fit->GetProb(); + printf("ROOT Chi-sq p-value = %.2f %s 0.05 | %s\n", pRoot, pRoot < 0.05 ? "<": ">", pRoot < 0.05 ? "reject" : "cannot reject"); + printf("################################################\n"); +} + +std::vector energy, height, sigma, lowE, highE ; +std::vector energyFlag, sigmaFlag; + +bool loadFitParameters(TString fitParaFile){ + + energy.clear(); energyFlag.clear(); + sigma.clear(); sigmaFlag.clear(); + + lowE.clear(); highE.clear(); + + height.clear(); + + bool paraRead = false; + + printf("====================================================================== \n"); + printf("----- loading fit parameters from : %s", fitParaFile.Data()); + std::ifstream file; + file.open(fitParaFile.Data()); + + if( !file){ + printf("\ncannot read file : %s \n", fitParaFile.Data()); + return 0; + } + + while( file.good()) { + std::string tempLine; + getline(file, tempLine); + + if( tempLine.substr(0, 1) == "#" ) continue; + if( tempLine.substr(0, 2) == "//" ) continue; + if( tempLine.size() < 5 ) continue; + + ///printf("%s\n", tempLine.c_str()); + + std::vector temp = SplitStrAF(tempLine, " "); + + if( temp.size() < 7 ) continue; + + energy.push_back( atof(temp[0].c_str())); + lowE.push_back( atof(temp[1].c_str())); + highE.push_back( atof(temp[2].c_str())); + energyFlag.push_back(atoi(temp[3].c_str())); + sigma.push_back( atof(temp[4].c_str())); + sigmaFlag.push_back( atoi(temp[5].c_str())); + height.push_back( atof(temp[6].c_str())); + + } + + printf("... done.\n"); + + int n = energy.size(); + TString limStr = "(fixed)"; + printf("%2s| %34s | %10s \n", "ID", "Peak [MeV]", "Sigma [MeV]"); + for( int j = 0; j < n; j ++){ + if( energyFlag[j] == 0 ) limStr.Form("(limited, %6.3f - %6.3f)", lowE[j], highE[j]); + printf("%2d| %7.4f %-26s | %7.4f (%5s) \n", j, energy[j], limStr.Data(), sigma[j], sigmaFlag[j] == 1 ? "fixed" : "free"); + } + + paraRead = true; + + printf("====================================================================== \n"); + + return paraRead; + +} + +TCanvas * NewCanvas(TString name, TString title, int divX, int divY, int padSizeX, int padSizeY){ + TCanvas * output = NULL; + if( gROOT->FindObjectAny(name) == NULL ){ + output = new TCanvas(name, title, divX * padSizeX, divY * padSizeY); + }else{ + output = (TCanvas *) gROOT->FindObjectAny(name) ; + output->Clear(); + } + output->Divide(divX, divY); + return output; +} + +void ScaleAndDrawHist(TH1F * hist, double xMin, double xMax){ + + if ( xMin != xMax ) hist->GetXaxis()->SetRangeUser(xMin, xMax); + int maxBin = hist->GetMaximumBin(); + double ymax = hist->GetBinContent(maxBin); + hist->GetYaxis()->SetRangeUser(0, 1.1 * ymax); + hist->Draw(); + +} + +void PlotResidual(TH1F * hist, TF1 * fit){ + + TH1F * hRes = (TH1F*) hist->Clone(); + hRes->GetListOfFunctions()->Clear(); + hRes->SetTitle("Residual"); + hRes->SetName("hRes"); + hRes->SetYTitle("Hist - fit"); + hRes->Sumw2(0); + hRes->Sumw2(1); + hRes->Add(fit, -1); + hRes->Draw(); + +} + + +//######################################## +//### Fit a Gauss + Pol-n +//######################################## +void fitGaussPol(TH1F * hist, double mean, double sigmaMax, int degPol, double xMin, double xMax, TString optStat = ""){ + + printf("=========================================================\n"); + printf("================ fit 1-Gauss + Pol-%d BG ================\n", degPol); + printf(" * mean Range +- 5 sigmaMax \n"); + printf(" * inital parameters of the polynomial is random/pow(10, i) \n"); + printf("==========================================================\n"); + + recentFitMethod = "fitGaussPol"; + + gStyle->SetOptStat(optStat); + TCanvas * cFitGaussPol = NewCanvas("cFitGaussPol", Form("fit Gauss & Pol-%d | fitGaussPol", degPol), 1, 2, 800, 350); + cFitGaussPol->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + const int nPar = 3 + degPol + 1; + + TString funcExp = "[0] * TMath::Gaus(x, [1], [2], 1)"; + for( int i = 0; i < degPol+1 ; i++){ + funcExp += Form(" + [%d]*TMath::Power(x,%d)", i+3 , i); + } + TF1 * fit = new TF1("fit", funcExp.Data(), xMin, xMax); + + double * para = new double[nPar]; + para[0] = 100 * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[1] = mean; + para[2] = sigmaMax/2.; + for( int i = 0 ; i < degPol+1; i++){ + para[3 + i ] = gRandom->Rndm()/TMath::Power(10, i); + } + + fit->SetLineWidth(2); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + fit->SetParLimits(0, 0, 1e9); + fit->SetParLimits(1, mean - 5*sigmaMax, mean + 5 * sigmaMax); + fit->SetParLimits(2, 0, sigmaMax); + + hist->Fit("fit", "Rq"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + double bw = hist->GetBinWidth(1); + + printf("histogram name : %s \n====== Gaussian:\ncount: %8.0f(%3.0f)\nmean : %8.4f(%8.4f)\nsigma: %8.4f(%8.4f) \n", + hist->GetName(), + paraA[0] / bw, paraE[0] /bw, + paraA[1], paraE[1], + paraA[2], paraE[2]); + + printf("------- the polnomail BG:\n"); + for(int i = 0 ; i < degPol+1; i++){ + printf("%2d | %8.4f(%8.4f) \n", i, paraA[3+i], paraE[3+i]); + } + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.12, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + text.DrawLatex(0.12, 0.75,Form("count: %4.0f(%3.0f)", paraA[0] / bw, paraE[0] /bw)); + text.DrawLatex(0.12, 0.70,Form("E_{x}: %6.3f(%5.3f) MeV", paraA[1], paraE[1])); + text.DrawLatex(0.12, 0.65,Form("#sigma: %3.0f(%3.0f) keV", paraA[2] * 1000., paraE[2] * 1000.)); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.60, 0.85 - 0.05*i ,Form("%3s: %8.3f(%8.3f)", Form("p%d", i), paraA[3+i], paraE[3+i])); + } + + TString expression = "[0] "; + for( int j = 1; j < degPol + 1; j++){ + expression += Form(" + [%d]*TMath::Power(x, %d)", j, j); + } + TF1 * g0 = new TF1("g0", expression.Data(), xMin, xMax); + for( int j = 0; j < degPol + 1 ; j++){ + g0->SetParameter(j, paraA[3+j]); + } + g0->SetLineColor(4); + g0->Draw("same"); + + TF1* f0 = new TF1("f0", "[0] * TMath::Gaus(x, [1], [2], 1)", xMin, xMax); + f0->SetParameter(0, paraA[0]); + f0->SetParameter(1, paraA[1]); + f0->SetParameter(2, paraA[2]); + f0->SetLineColor(2); + f0->SetNpx(1000); + f0->Draw("same"); + +// GoodnessofFit(hist, fit); + + cFitGaussPol->cd(2); + PlotResidual(hist, fit); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + BestFitCount.push_back(paraA[0]); + BestFitMean.push_back(paraA[1]); + BestFitSigma.push_back(paraA[2]); + +} + +//######################################## +//#### fit 2 gauss + pol-1 // not updated +//######################################## +std::vector fit2GaussP1(TH1F * hist, double mean1, double sigma1, + double mean2, double sigma2, + double xMin, double xMax, TString optStat = "", bool newCanvas = false){ + + + printf("=========================================================\n"); + printf("================ fit 2-Gauss + Pol-1 BG ================\n" ); + printf(" NOT updated. It works, but the code is old \n"); + printf("==========================================================\n"); + + recentFitMethod = "fit2GaussP1"; + + std::vector output; + output.clear(); + + gStyle->SetOptStat(optStat); + TCanvas * cFit2GaussP1 = NewCanvas("cFit2GaussP1", "fit Gauss & P1 | fit2GaussP1", 1, 1, 800, 350); + cFit2GaussP1->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + TF1 * fit = new TF1("fit", "[0] * TMath::Gaus(x, [1], [2], 1) + [3] * TMath::Gaus(x, [4], [5], 1) + [6] + [7]*x", xMin, xMax); + + double * para = new double[8]; + para[0] = 20 * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[1] = mean1; + para[2] = sigma1; + para[3] = 100 * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[4] = mean2; + para[5] = sigma2; + para[6] = 1; + para[7] = 0; + + fit->SetLineWidth(2); + fit->SetLineColor(2); + fit->SetNpx(1000); + fit->SetParameters(para); + + hist->Fit("fit", "Rq"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + double bw = hist->GetBinWidth(1); + + printf("%7s ====== count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + hist->GetName(), + paraA[0] / bw, paraE[0] /bw, + paraA[1], paraE[1], + paraA[2], paraE[2]); + printf("%7s ====== count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + "", + paraA[3] / bw, paraE[3] /bw, + paraA[4], paraE[4], + paraA[5], paraE[5]); + + output.push_back( paraA[0]/bw); + output.push_back( paraE[0]/bw); + output.push_back( paraA[1]); + output.push_back( paraE[1]); + output.push_back( paraA[2]); + output.push_back( paraE[2]); + + output.push_back( paraA[3]/bw); + output.push_back( paraE[3]/bw); + output.push_back( paraA[4]); + output.push_back( paraE[4]); + output.push_back( paraA[5]); + output.push_back( paraE[5]); + + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + text.DrawLatex(0.15, 0.75,Form("count: %4.0f(%3.0f), E_{x}: %6.3f(%5.3f) MeV, #sigma: %3.0f(%3.0f) keV ", + paraA[0] / bw, paraE[0] /bw, + paraA[1], paraE[1], + paraA[2] * 1000., paraE[2] * 1000.)); + text.DrawLatex(0.15, 0.7, Form("count: %4.0f(%3.0f), E_{x}: %6.3f(%5.3f) MeV, #sigma: %3.0f(%3.0f) keV ", + paraA[3] / bw, paraE[3] /bw, + paraA[4], paraE[4], + paraA[5] * 1000., paraE[5] * 1000.)); + + text.DrawLatex(0.15, 0.6, Form("Line : %6.3f(%5.3f) + %6.3f(%5.3f)x ", + paraA[6], paraE[6], + paraA[7], paraE[7])); + + GoodnessofFit(hist, fit); + + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < 2; i++){ + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } + + return output; +} + + +//######################################## +//#### fit for gamma + pol-n BG +//######################################## +void fitGF3Pol(TH1F * hist, double mean, double sigmaMax, double ratio, double beta, double step, int degPol, double xMin, double xMax, TString optStat = ""){ + + printf("=========================================================\n"); + printf("================ fit GF1 + Pol-%d BG ================\n", degPol); + printf(" * mean Range = xMin, xMax \n"); + printf(" * inital parameters of the polynomial is random/pow(10, i) \n"); + printf("==========================================================\n"); + + recentFitMethod = "fitGF3Pol"; + + gStyle->SetOptStat(optStat); + + gStyle->SetOptStat(optStat); + TCanvas * cFitGF3Pol = NewCanvas("cFitGF3Pol", Form("fit GF3 + pol-%d | fitGF3Pol", degPol), 1, 2, 800, 350); + cFitGF3Pol->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + nPeaks = 1; + nPols = degPol; + int nPar = 6*nPeaks + degPol + 1; + + TF1 * fit = new TF1("fit", nGF3Pol, xMin, xMax, nPar); + + fit->Print(); + + double * para = new double[nPar]; + para[0] = hist->GetMaximum() *4; + para[1] = mean; + para[2] = sigmaMax/2.; + para[3] = ratio ; + para[4] = beta ; + para[5] = step ; + for( int i = 0 ; i < degPol + 1; i++){ + para[6+i] = gRandom->Rndm()/TMath::Power(10, i); + } + + fit->SetLineWidth(2); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + fit->SetParLimits(0, 0, 1e9); + fit->SetParLimits(1, xMin, xMax); + fit->SetParLimits(2, 0.00000001, sigmaMax); + fit->SetParLimits(3, 0, 0.5); + fit->SetParLimits(4, 1, 400); + fit->SetParLimits(5, 0, 0.5); + + hist->Fit("fit", "Rq"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + double chisquare = fit->GetChisquare(); + int ndf = fit->GetNDF(); + double bw = hist->GetBinWidth(1); + + printf("histogram : %s \n", hist->GetName()); + printf("========= The Gaussian \n"); + printf("count: %8.0f(%3.0f)\n", paraA[0] / bw, paraE[0] /bw); + printf("mean : %8.4f(%8.4f)\n", paraA[1], paraE[1]); + printf("sigma: %8.4f(%8.4f)\n", paraA[2], paraE[2]); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + text.DrawLatex(0.12, 0.65, Form("count : %5.0f(%5.0f)", paraA[0]/bw, paraE[0]/bw)); + text.DrawLatex(0.12, 0.60, Form(" mean : %5.3f(%5.3f) keV", paraA[1], paraE[1])); + text.DrawLatex(0.12, 0.55, Form("sigma : %5.3f(%5.3f) keV", paraA[2], paraE[2])); + text.DrawLatex(0.12, 0.50, Form(" FWHM : %5.3f(%5.3f) keV", paraA[2] *2.355, paraE[2]*2.355)); + + text.DrawLatex(0.12, 0.40, Form("#chi^2/ndf : %5.3f", chisquare/ndf)); + + //GoodnessofFit(hist, fit); + + /// 0 1 2 3 4 5 + std::string label[8] = {"Area", "mean", "sigma", "ratio", "beta", "step"}; + printf("---------- The detail\n"); + for(int i = 0 ; i < 6 ; i++){ + printf("%d | %8s | %f (%f) \n", i, label[i].c_str(), paraA[i], paraE[i]); + text.DrawLatex(0.65, 0.85-0.05*i, Form("%6s: %5.3f(%5.3f)", label[i].c_str(), paraA[i], paraE[i])); + } + for(int i = 6 ; i < nPar; i++){ + printf("%d | %8s | %f (%f) \n", i, Form("p%d", i-6), paraA[i], paraE[i]); + text.DrawLatex(0.65, 0.85-0.05*i, Form("%6s: %5.3f (%5.3f) \n", Form("p%d", i-6), paraA[i], paraE[i])); + } + + /// norm * (1-ratio)* TMath::Gaus(x[0], mean, sigma, 1) + TF1 * g0 = new TF1("g0", "[0] * (1.0-[3]) * TMath::Gaus(x, [1], [2], 1)", xMin, xMax); + g0->SetParameter(0, paraA[0]); + g0->SetParameter(1, paraA[1]); + g0->SetParameter(2, paraA[2]); + g0->SetParameter(3, paraA[3]); + g0->SetNpx(1000); + g0->SetLineColor(kRed); + + /// norm * ratio * exp( sigma * sigma/2/beta/beta)* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ; + TF1 * g1 = new TF1("g1", "[0] * [3] * exp( [2] * [2]/2/[4]/[4]) / (2* [4])* exp((x-[1])/[4]) * TMath::Erfc( (x-[1])/([2] * sqrt(2)) + [2]/[4]/sqrt(2)) ", xMin, xMax); + g1->SetParameter(0, paraA[0]); + g1->SetParameter(1, paraA[1]); + g1->SetParameter(2, paraA[2]); + g1->SetParameter(3, paraA[3]); + g1->SetParameter(4, paraA[4]); + g1->SetNpx(1000); + g1->SetLineColor(kGreen +3); + + /// norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) ); + TF1 * g2 = new TF1("g2", "[0] * [3] * TMath::Erfc( (x-[1])/([2] * sqrt(2)) );", xMin, xMax); + g2->SetParameter(0, paraA[0]); + g2->SetParameter(1, paraA[1]); + g2->SetParameter(2, paraA[2]); + g2->SetParameter(3, paraA[5]); + g2->SetNpx(1000); + g2->SetLineColor(kViolet); + + + TString expression = "[0] "; + for( int j = 1; j < degPol + 1; j++){ + expression += Form(" + [%d]*TMath::Power(x, %d)", j, j); + } + TF1 * g3 = new TF1("g3", expression.Data(), xMin, xMax); + for( int j = 0; j < degPol + 1 ; j++){ + g3->SetParameter(j, paraA[6+j]); + } + g3->SetLineColor(kBlue); + g3->Draw("same"); + + + g0->Draw("same"); + g1->Draw("same"); + g2->Draw("same"); + g3->Draw("same"); + + cFitGF3Pol->cd(2); + PlotResidual(hist, fit); + +} + +//############################################## +//##### Auto Fit n-Gauss with estimated BG +//############################################## +std::vector fitAuto(TH1F * hist, int bgEst = 10, + double peakThreshold = 0.05, + double sigmaMax = 0, + int peakDensity = 4, + TString optStat = ""){ + + printf("================================================================\n"); + printf("========== Auto Fit n-Gauss with estimated BG ==================\n"); + printf(" * bgEst = parameter of BG estimation, larger BG, more linear \n"); + printf(" * peakThreshold = precentage of the highest peak that will count \n"); + printf(" * sigmaMax = maximum sigma, if -1, program try to find the sigma \n"); + printf(" * peakDensity = peak will closer when the number is larger "); + printf(" \n"); + printf(" after peaks found, the i-th peaks will be limited by the mid-point\n"); + printf(" by the (i-1)-th peak and the i-th peak, and the mid-point of the\n"); + printf(" i-th peak and (i+1)-th peak \n"); + printf(" i.e. [peak(i-1)+peak(i)]/2 < limit of peak(i) < [peak(i)+peak(i+1)]/2 \n"); + printf("================================================================\n"); + + recentFitMethod = "fitAuto"; + + gStyle->SetOptStat(optStat); + TCanvas *cFitAuto = NewCanvas("cFitAuto","Auto Fitting | fitAuto", 1, 4, 800, 300); + cFitAuto->cd(1); + + ScaleAndDrawHist(hist, 0, 0); + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitted spectrum (BG=%d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + ///specS->GetXaxis()->SetTitleSize(0.06); + ///specS->GetYaxis()->SetTitleSize(0.06); + ///specS->GetXaxis()->SetTitleOffset(0.7); + ///specS->GetYaxis()->SetTitleOffset(0.6); + + //=================== find peak and fit + gStyle->SetOptFit(0); + TSpectrum * peak = new TSpectrum(50); + nPeaks = peak->Search(hist, peakDensity, "", peakThreshold); + + if( bgEst > 0 ) { + printf("============= estimating background...\n"); + TH1 * h1 = peak->Background(hist, bgEst); + h1->Draw("same"); + printf("============= substracting the linear background...\n"); + specS->Sumw2(); + specS->Add(h1, -1.); + } + + cFitAuto->cd(2)->SetGrid(); + cFitAuto->cd(2); + specS->Draw("hist"); + + //========== Fitting + printf("============= Fitting....."); + printf(" found %d peaks \n", nPeaks); + + double * xpos = peak->GetPositionX(); + double * ypos = peak->GetPositionY(); + + int * inX = new int[nPeaks]; + TMath::Sort(nPeaks, xpos, inX, 0 ); + std::vector energy, height; + for( int j = 0; j < nPeaks; j++){ + energy.push_back(xpos[inX[j]]); + height.push_back(ypos[inX[j]]); + } + for( int j = 0; j < nPeaks; j++){ + printf(" energy : %f , %f \n", energy[j], height[j]); + } + + + if( sigmaMax == 0 ){ + printf("------------- Estimate sigma for each peak \n"); + sigma.clear(); + int binMin = hist->FindBin(xMin); + int binMax = hist->FindBin(xMax); + for( int i = 0; i < nPeaks ; i++){ + int b0 = hist->FindBin(energy[i]); + double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.; + //---- backward search, stop when + for( int b = b0-1 ; b > binMin ; b-- ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y < (height[i])/2. ) { + sMin = energy[i] - hist->GetBinCenter(b); + break; + } + } + //---- forward search, stop when + for( int b = b0+1 ; b < binMax ; b++ ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y < (height[i])/2. ) { + sMax = hist->GetBinCenter(b) - energy[i]; + break; + } + } + + double temp = TMath::Min(sMin, sMax); + /// When there are multiple peaks closely packed : + if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back(); + sigma.push_back(temp); + + printf("%2d | x : %8.2f | sigma(est) %f \n", i, energy[i], sigma[i]); + } + }else if( sigmaMax < 0 ){ + printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax)); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax)); + }else if( sigmaMax > 0 ){ + printf("========== use user input sigma : %f/2. \n", sigmaMax/2.); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.); + } + + + + int numParPerPeak = 3; + const int n = numParPerPeak * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[numParPerPeak*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[numParPerPeak*i+1] = energy[i]; + if( sigmaMax == 0 ){ + para[numParPerPeak*i+2] = sigma[i]; + }else if(sigmaMax < 0 ){ + para[numParPerPeak*i+2] = abs(sigmaMax); + }else if(sigmaMax > 0 ){ + para[numParPerPeak*i+2] = sigmaMax/2.; + } + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3 * nPeaks ); + fit->SetLineWidth(2); + fit->SetLineColor(2); + fit->SetNpx(1000); + fit->SetParameters(para); + + if( nPeaks > 1 ){ + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(numParPerPeak*i+0, 0, 1e+9); + double de1 = 1, de2 = 1; + if( i == 0 ){ + de2 = (energy[i+1] - energy[i])/2.; + de1 = de2; + }else if( i < nPeaks -1 ){ + de1 = (energy[i] - energy[i-1])/2.; + de2 = (energy[i+1] - energy[i])/2.; + }else{ + de1 = (energy[i] - energy[i-1])/2.; + de2 = de1; + } + + fit->SetParLimits(numParPerPeak*i+1, energy[i] - de1 , energy[i] + de2); + if( sigmaMax== 0 ) fit->SetParLimits(numParPerPeak*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma + if( sigmaMax < 0 ) fit->FixParameter(numParPerPeak*i+2, abs(sigmaMax)); + if( sigmaMax > 0 ) fit->SetParLimits(numParPerPeak*i+2, 0, sigmaMax); + } + }else{ + fit->SetParLimits(0, 0, 1e+9); + fit->SetParLimits(2, 0, sigmaMax); + } + + specS->Fit("fit", "q"); + + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + //======== calculate reduce chi-squared + //GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + std::vector exPos; + for(int i = 0; i < nPeaks ; i++){ + exPos.push_back(paraA[numParPerPeak*i+1]); + printf("%2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[numParPerPeak*i] / bw, paraE[numParPerPeak*i] /bw, + paraA[numParPerPeak*i+1], paraE[numParPerPeak*i+1], + paraA[numParPerPeak*i+2], paraE[numParPerPeak*i+2]); + + } + cFitAuto->Update(); + + //draw the indivual fit + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[numParPerPeak*i]); + gFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + gFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + + //======== print text on plot + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + cFitAuto->cd(3); + PlotResidual(specS, fit); + + cFitAuto->cd(4); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } + + return exPos; + +} + + +//######################################## +//###### NOT tested +//######################################## +std::vector fitNGF3(TH1 * hist, int bgEst = 10, + double peakThreshold = 0.1, + double sigmaMax = 20, + int peakDensity = 4, + TString optStat = "", bool newPlot = true){ + + TCanvas *cFitAuto = NULL; + if( newPlot ){ + cFitAuto = new TCanvas("cFitAuto","Auto Fitting", 100, 100, 800,800); + cFitAuto->Divide(1,2); + + gStyle->SetOptStat(optStat); + cFitAuto->cd(1); + hist->Draw(); + } + + recentFitMethod = "fitNGF3"; + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitted spectrum (BG=%d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + ///specS->GetXaxis()->SetTitleSize(0.06); + ///specS->GetYaxis()->SetTitleSize(0.06); + ///specS->GetXaxis()->SetTitleOffset(0.7); + ///specS->GetYaxis()->SetTitleOffset(0.6); + + //=================== find peak and fit + gStyle->SetOptFit(0); + TSpectrum * peak = new TSpectrum(50); + nPeaks = peak->Search(hist, peakDensity, "", peakThreshold); + + if( bgEst > 0 ) { + printf("============= estimating background...\n"); + TH1 * h1 = peak->Background(hist, bgEst); + h1->Draw("same"); + printf("============= substracting the linear background...\n"); + specS->Sumw2(); + specS->Add(h1, -1.); + } + + if( newPlot ){ + cFitAuto->cd(2)->SetGrid(); + cFitAuto->cd(2); + } + specS->Draw("hist"); + + + //========== Fitting + if( newPlot ){ + printf("============= Fitting....."); + printf(" found %d peaks \n", nPeaks); + } + double * xpos = peak->GetPositionX(); + double * ypos = peak->GetPositionY(); + + int * inX = new int[nPeaks]; + TMath::Sort(nPeaks, xpos, inX, 0 ); + std::vector energy, height; + for( int j = 0; j < nPeaks; j++){ + energy.push_back(xpos[inX[j]]); + height.push_back(ypos[inX[j]]); + } + if( newPlot ){ + for( int j = 0; j < nPeaks; j++){ + printf(" energy : %f , %f \n", energy[j], height[j]); + } + } + + int numParPerPeak = 6; + const int n = numParPerPeak * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[numParPerPeak*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[numParPerPeak*i+1] = energy[i]; + para[numParPerPeak*i+2] = sigmaMax; + para[numParPerPeak*i+3] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()) * 0.1; + para[numParPerPeak*i+4] = sigmaMax; + para[numParPerPeak*i+5] = -4; + } + + TF1 * fit = new TF1("fit", nGF3, xMin, xMax, numParPerPeak * nPeaks ); + fit->SetLineWidth(2); + fit->SetLineColor(2); + fit->SetNpx(1000); + fit->SetParameters(para); + + if( nPeaks > 1 ){ + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(numParPerPeak*i+0, 0, 1e+9); + double de1 = 1, de2 = 1; + if( i == 0 ){ + de2 = (energy[i+1] - energy[i])/2.; + de1 = de2; + }else if( i < nPeaks -1 ){ + de1 = (energy[i] - energy[i-1])/2.; + de2 = (energy[i+1] - energy[i])/2.; + }else{ + de1 = (energy[i] - energy[i-1])/2.; + de2 = de1; + } + + fit->SetParLimits(numParPerPeak*i+1, energy[i] - de1 , energy[i] + de2); + fit->SetParLimits(numParPerPeak*i+2, 0, sigmaMax * 5); + + fit->SetParLimits(numParPerPeak*i+3, 0, 1e+9); + fit->SetParLimits(numParPerPeak*i+4, 0, sigmaMax); + fit->SetParLimits(numParPerPeak*i+5, -10, -2); + + } + }else{ + fit->SetParLimits(0, 0, 1e+9); + fit->SetParLimits(2, 0, sigmaMax); + fit->SetParLimits(3, 0, 1e+9); + fit->SetParLimits(4, 0, sigmaMax); + fit->SetParLimits(5, -10, -2); + } + specS->Fit("fit", "q"); + + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + //======== calculate reduce chi-squared + if( newPlot ) GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + std::vector exPos; + for(int i = 0; i < nPeaks ; i++){ + exPos.push_back(paraA[numParPerPeak*i+1]); + double totCount = paraA[numParPerPeak*i] + paraA[numParPerPeak*i+3]; + double totCountErr = sqrt(paraE[numParPerPeak*i]*paraE[numParPerPeak*i] + paraE[numParPerPeak*i+3]*paraE[numParPerPeak*i+3]); + printf("%2d , count: %8.0f(%3.0f)+%8.0f(%3.0f)=%8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f), skewneww: %4.1f(%4.1f) \n", + i, + paraA[numParPerPeak*i] / bw, paraE[numParPerPeak*i] /bw, + paraA[numParPerPeak*i+3] / bw, paraE[numParPerPeak*i+3] /bw, + totCount / bw, totCountErr /bw, + paraA[numParPerPeak*i+1], paraE[numParPerPeak*i+1], + paraA[numParPerPeak*i+2], paraE[numParPerPeak*i+2], + paraA[numParPerPeak*i+5], paraE[numParPerPeak*i+5]); + + //PrintPar(fit, numParPerPeak); + } + if( newPlot ) cFitAuto->Update(); + + //draw the indivual fit + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + TF1 ** kFit = new TF1 *[nn]; + TF1 ** zFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1) + [3] * TMath::Gaus(x,[1],[4], 1) * ( 1 + TMath::Erf( [5]*(x-[1])/sqrt(2)/[4] ))", xMin, xMax); + gFit[i]->SetParameter(0, paraA[numParPerPeak*i]); + gFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + gFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]); + gFit[i]->SetParameter(3, paraA[numParPerPeak*i+3]); + gFit[i]->SetParameter(4, paraA[numParPerPeak*i+4]); + gFit[i]->SetParameter(5, paraA[numParPerPeak*i+5]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + + kFit[i] = new TF1(Form("kFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1) * ( 1 + TMath::Erf( [3]*(x-[1])/sqrt(2)/[2] ))", xMin, xMax); + kFit[i]->SetParameter(0, paraA[numParPerPeak*i+3]); + kFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + kFit[i]->SetParameter(2, paraA[numParPerPeak*i+4]); + kFit[i]->SetParameter(3, paraA[numParPerPeak*i+5]); + kFit[i]->SetLineColor(i+1); + kFit[i]->SetNpx(1000); + kFit[i]->SetLineWidth(1); + kFit[i]->Draw("same"); + + zFit[i] = new TF1(Form("zFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + zFit[i]->SetParameter(0, paraA[numParPerPeak*i]); + zFit[i]->SetParameter(1, paraA[numParPerPeak*i+1]); + zFit[i]->SetParameter(2, paraA[numParPerPeak*i+2]); + zFit[i]->SetLineColor(i+1); + zFit[i]->SetNpx(1000); + zFit[i]->SetLineWidth(1); + zFit[i]->Draw("same"); + + } + + specS->Draw("hist same"); + + return exPos; + +} + +//######################################## +//###### fit N Gauss with estimated BG +//######################################## +void fitNGauss(TH1F * hist, int bgEst = -1, TString fitFile = "AutoFit_para.txt", TString optStat = ""){ + + printf("================================================================\n"); + printf("================ fit N-Gauss with estimated BG ================\n"); + printf(" * bgEst = larger of bgEst, more linear the estimated BG \n"); + printf(" * need the file input \n"); + printf(" \n"); + printf(" 1) The histogram will be subtracted by the estimated BG. \n"); + printf(" 2) n-Gauss will then be fitted the BG subtracted histogram \n"); + printf("================================================================\n"); + + recentFitMethod = "fitNGauss"; + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + nPeaks = energy.size(); + + gStyle->SetOptStat(optStat); + TCanvas *cFitNGauss = NewCanvas("cFitNGauss","Fit n-Gauss | fitNGauss", 1,4, 800, 300);; + cFitNGauss->cd(1); + + ScaleAndDrawHist(hist, 0, 0); + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitNGauss (BG = %2d); Ex [MeV]; Count / %4.0f keV", bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + + //=================== find peak and fi + + gStyle->SetOptFit(0); + //cFitNGauss->cd(2)->SetGrid(); + cFitNGauss->cd(2); + + if( bgEst > 0 ) { + printf("============= estimating background...\n"); + TSpectrum * peak = new TSpectrum(50); + TH1 * h1 = peak->Background(hist, bgEst); + cFitNGauss->cd(1); + h1->Draw("same"); + cFitNGauss->cd(2); + printf("============= substracting the estimated background...\n"); + specS->Sumw2(); + specS->Add(h1, -1.); + } + + specS->Draw("hist"); + + //========== Fitting + printf("============= Fitting %d-Gauss..... \n", nPeaks); + + const int n = 3 * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3* nPeaks ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + specS->Fit("fit", "q"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + //======== calculate reduce chi-squared + //GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + //draw the indivual fit + specS->Draw("hist"); + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + //specS->Draw("E same"); + + + //======== print text on plot + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + cFitNGauss->cd(3); + PlotResidual(specS, fit); + + cFitNGauss->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } + + cFitNGauss->Update(); +} + + +//######################################## +//#### not updated +//######################################## +void fitNGaussSub(TH1F * hist, int bgEst = 10, int degPol = 1, TString fitFile = "AutoFit_para.txt", TString optStat = ""){ + + + printf("==================================================================\n"); + printf("======== fit N-Gauss with estimated BG (method-2) ================\n"); + printf(" * bgEst = larger of bgEst, more linear the estimated BG \n"); + printf(" * degPol = degree of polynomial \n"); + printf(" * need the file input \n"); + printf(" \n"); + printf(" 1) A BG is estimated, and then the BG is fitted by a polynomial. \n"); + printf(" 2) The histogram will be subtracted by the polynomial. \n"); + printf(" 3) n-Gauss will then be fitted the subtracted histogram \n"); + printf("================================================================\n"); + + recentFitMethod = "fitNGaussSub"; + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + nPeaks = energy.size(); + nPols = degPol; + + gStyle->SetOptStat(optStat); + TCanvas *cFitNGaussSub = NewCanvas("cFitNGaussSub","Fit n-Gauss, replace estimated BG with Pol-n | fitNGauss2", 1, 4, 800, 300 ); + //if(! cFitNGaussSub->GetShowEventStatus() ) cFitNGaussSub->ToggleEventStatus(); + + cFitNGaussSub->cd(1); + ScaleAndDrawHist(hist, 0, 0); + + TH1F * specS = (TH1F*) hist->Clone(); + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + TString titleH; + titleH.Form("fitNGauss2 (replace Est. BG with Pol-%d) (BG = %2d); Ex [MeV]; Count / %4.0f keV", degPol, bgEst, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + + printf("============= estimating background...\n"); + TSpectrum * peak = new TSpectrum(50); + TH1 * h1 = peak->Background(hist, bgEst); + + printf("============= fit the est-background with a polynomial function...\n"); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + + bg->SetParameter(0, 50); + bg->SetParameter(0, 0); + bg->SetLineColor(2); + bg->SetNpx(1000); + h1->Fit("bg", "R", ""); + + hist->Draw(); + bg->Draw("same"); + + //======== print text on plot + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + + const Double_t * paraAt = bg->GetParameters(); + const Double_t * paraEt = bg->GetParErrors(); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[i], paraEt[i])); + } + + gStyle->SetOptFit(0); +// cFitNGaussSub->cd(2)->SetGrid(); + cFitNGaussSub->cd(2); + + printf("============= substracting the polynomial background...\n"); + specS->Sumw2(); + specS->Add(bg, -1.); + specS->Draw("hist"); + + //========== Fitting + printf("============= Fitting..... \n"); + + const int n = 3 * nPeaks; + double * para = new double[n]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, 3* nPeaks ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + specS->Fit("fit", "q"); + + const Double_t* paraE = fit->GetParErrors(); + const Double_t* paraA = fit->GetParameters(); + + GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + //draw the indivual fit + specS->Draw("hist"); + fit->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + //specS->Draw("E same"); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + cFitNGaussSub->cd(3); + PlotResidual(specS, fit); + + cFitNGaussSub->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } + +} + + +//######################################## +//#### fit N-Gauss with pol-n BG +//######################################## +void fitNGaussPol(TH1F * hist, int degPol, TString fitFile = "AutoFit_para.txt",double xMin = 0, double xMax = 0, TString optStat = ""){ + + printf("================================================================\n"); + printf("================ fit N-Gauss with Pol-%1d BG ==================\n", degPol); + printf(" * degPol = degree of polynomial \n"); + printf(" * need the file input \n"); + printf(" * xMin, xMax = if left empty, full range will be used\n"); + printf(" \n"); + printf(" 1) The histogram will be fitted by n-Gauss + Pol \n"); + printf("================================================================\n"); + + recentFitMethod = "fitNGaussPol"; + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + gStyle->SetOptStat(optStat); + nPeaks = energy.size(); + nPols = degPol; + + TCanvas * cFitNGaussPol = NewCanvas("cFitNGaussPol", Form("Fitting with n-Gauss + pol-%d | fitNGaussPol", degPol), 1, 3, 800, 300); + //if(! cFitNGaussPol->GetShowEventStatus() ) cFitNGaussPol->ToggleEventStatus(); + cFitNGaussPol->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + if( xMin == xMax){ + xMin = hist->GetXaxis()->GetXmin(); + xMax = hist->GetXaxis()->GetXmax(); + } + int xBin = hist->GetXaxis()->GetNbins(); + + + printf("============= find the polynomial background ..... \n"); + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + for(int i = 0; i < nPols + 1; i++){ + //if( ggg == NULL ){ + para[3*nPeaks+i] = gRandom->Rndm()/TMath::Power(10, 3*i); + //}else{ + // para[3*nPeaks+i] = gPara[i]; + //} + } + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + hist->Fit("fit", "Rq"); + + //=========== get the polynomial part + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + for( int i = 0; i < degPol + 1; i++){ + bg->SetParameter(i, paraA[3*nPeaks+i]); + } + bg->SetNpx(1000); + + for( int i = 0; i < degPol + 1; i++){ + printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]); + } + printf("====================================\n"); + + cFitNGaussPol->cd(1); + bg->Draw("same"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + //GoodnessofFit(specS, fit); + + double bw = hist->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + cFitNGaussPol->Update(); + + cFitNGaussPol->cd(2); + PlotResidual(hist, fit); + + cFitNGaussPol->cd(3); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } +} + + + + +//######################################## +//#### fit N-Gauss with pol-n BG +//######################################## +void fitNGaussPolSub(TH1F * hist, int degPol, TString fitFile = "AutoFit_para.txt",double xMin = 0, double xMax = 0, TString optStat = ""){ + + printf("================================================================\n"); + printf("========= fit N-Gauss with Pol-%1d BG / 2nd method ============\n", degPol); + printf(" * degPol = degree of polynomial \n"); + printf(" * need the file input \n"); + printf(" * xMin, xMax = if left empty, full range will be used\n"); + printf(" \n"); + printf(" 1) The histogram will be fitted by n-Gauss + Pol -> to get the estimated BG \n"); + printf(" 2) The histogram will be subtracted by the Pol BG. \n"); + printf(" 3) n-Gauss will then be fitted the BG subtracted histogram \n"); + printf("================================================================\n"); + + recentFitMethod = "fitNGaussPolSub"; + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + gStyle->SetOptStat(optStat); + nPeaks = energy.size(); + nPols = degPol; + + TCanvas * cFitNGaussPolSub = NewCanvas("cFitNGaussPolSub", Form("Fitting with n-Gauss + pol-%d (method-2) | fitGaussPol2", degPol), 1, 4, 800, 300); + //if(! cFitNGaussPol->GetShowEventStatus() ) cFitNGaussPol->ToggleEventStatus(); + cFitNGaussPolSub->cd(1); + + ScaleAndDrawHist(hist, xMin, xMax); + + if( xMin == xMax){ + xMin = hist->GetXaxis()->GetXmin(); + xMax = hist->GetXaxis()->GetXmax(); + } + int xBin = hist->GetXaxis()->GetNbins(); + + printf("============= find the polynomial background ..... \n"); + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + for(int i = 0; i < nPols + 1; i++){ + para[3*nPeaks+i] = gRandom->Rndm()/TMath::Power(10, 3*i); + } + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + hist->Fit("fit", "nq"); + + //=========== get the polynomial part and substract + const Double_t* paraAt = fit->GetParameters(); + const Double_t* paraEt = fit->GetParErrors(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + for( int i = 0; i < degPol + 1; i++){ + bg->SetParameter(i, paraAt[3*nPeaks+i]); + } + bg->SetNpx(1000); + + for( int i = 0; i < degPol + 1; i++){ + printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[3*nPeaks+i], paraEt[3*nPeaks+i]); + } + printf("====================================\n"); + + cFitNGaussPolSub->cd(1); + bg->Draw("same"); + + TH1F * specS = (TH1F*) hist->Clone(); + TString titleH; + titleH.Form("pol-%d BG Subtracted spectrum (fitNGaussPol-2); Ex [MeV]; Count / %4.0f keV", degPol, (xMax-xMin)*1000./xBin ); + specS->SetTitle(titleH); + specS->SetName("specS"); + + //=================== find peak and fit + + gStyle->SetOptFit(0); + ///cFitNGaussPol->cd(2)->SetGrid(); + cFitNGaussPolSub->cd(2); + + printf("============= substracting the polynomial background...\n"); + specS->Sumw2(); + specS->Add(bg, -1.); + specS->Draw("hist"); + + //======= fit again + printf("============= fitting the subtracted spectrum.... \n"); + + nPar = 3* nPeaks; + TF1 * fita = new TF1("fita", nGauss, xMin, xMax, nPar ); + fita->SetLineWidth(3); + fita->SetLineColor(1); + fita->SetNpx(1000); + fita->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fita->SetParLimits(3*i , 0, 1e9); + + if( energyFlag[i] == 1 ) { + fita->FixParameter(3*i+1, energy[i]); + }else{ + fita->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fita->FixParameter(3*i+2, sigma[i]); + }else{ + fita->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + specS->Fit("fita", "q"); + + const Double_t* paraE = fita->GetParErrors(); + const Double_t* paraA = fita->GetParameters(); + + //GoodnessofFit(specS, fit); + + double bw = specS->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + //draw the indivual fit + specS->Draw("hist"); + fita->Draw("same"); + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + specS->Draw("hist same"); + //specS->Draw("E same"); + + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraAt[3*nPeaks+i], paraEt[3*nPeaks+i])); + } + + double chi2 = fita->GetChisquare(); + int ndf = fita->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + cFitNGaussPolSub->cd(3); + PlotResidual(specS, fita); + + cFitNGaussPolSub->cd(4); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } + + +} + +//################################################# +//#### fit N-Gauss with pol-n BG using mouse click +//################################################# + +int nClick = 0; +bool peakFlag = 1; +std::vector xPeakList; +std::vector yPeakList; +std::vector xBGList; +std::vector yBGList; + +TH1F * tempHist; +int markerStyle = 23; +int markerColor = 2; +int markerBGColor = 4; + +void Clicked() { + int event = gPad->GetEvent(); + if (event != 11) return; + TObject *select = gPad->GetSelected(); + if (!select) return; + + TH1F *h = (TH1F*)select; + int px = gPad->GetEventX(); + double xd = gPad->AbsPixeltoX(px); + float x = gPad->PadtoX(xd); + + if( peakFlag ) { + xPeakList.push_back(x); + }else{ + xBGList.push_back(x); + } + int b = tempHist->FindBin(x); + double y = tempHist->GetBinContent(b); + if( peakFlag ){ + yPeakList.push_back(y); + }else{ + yBGList.push_back(y); + } + // add marker in the histogram + TMarker * mark = new TMarker(x,y, markerStyle); + mark->SetMarkerColor(markerColor); + tempHist->GetListOfFunctions()->Add(mark); + + printf("%2d | x : %8.2f , y : %.0f \n", nClick, x, y); + + nClick ++; + +} + + +void SaveFitPara(bool isBestFit = true, TString fileName = "AutoFit_para.txt"){ + + if( xPeakList.size() == 0 && BestFitMean.size() == 0 ){ + printf(" no fit paramaters availible. \n"); + return; + } + + if( recentFitMethod == "fitGF3Pol" || recentFitMethod == "fitGF3" ){ + printf(" Not support for GF3 fitting. \n"); + return; + } + + printf("Save to : %s \n", fileName.Data()); + FILE * file_out; + file_out = fopen (fileName.Data(), "w+"); + + fprintf(file_out, "# for n-Gauss fit, can use \"#\", or \"//\" to comment out whole line\n"); + fprintf(file_out, "# peak low high fixed? sigma_Max fixed? hight\n"); + + + if ( xPeakList.size() == 0 || isBestFit ){ + for( int i = 0 ; i < BestFitMean.size() ; i++){ + fprintf(file_out, "%.3f %.3f %.3f 0 %.3f 0 %.0f\n", + BestFitMean[i], + BestFitMean[i] - 5*BestFitSigma[i], + BestFitMean[i] + 5*BestFitSigma[i], + BestFitSigma[i], + BestFitCount[i]); + } + }else{ + for( int i = 0 ; i < xPeakList.size() ; i++){ + fprintf(file_out, "%.3f %.3f %.3f 0 %.3f 0 %.0f\n", + xPeakList[i], + xPeakList[i] - 5*sigma[i], + xPeakList[i] + 5*sigma[i], + sigma[i], + yPeakList[i]); + } + } + fclose(file_out); +} + + +void clickFitNGaussPol(TH1F * hist, int degPol, double sigmaMax = 0, double meanRange = 0){ + + printf("=========================================================\n"); + printf("======== fit n-Gauss + Pol-%d BG using mouse click =====\n", degPol ); + printf("==========================================================\n"); + + recentFitMethod = "clickFitNGaussPol"; + + gStyle->SetOptStat(""); + gStyle->SetOptFit(0); + + TCanvas * cClickFitNGaussPol = NULL; + if( gROOT->FindObjectAny("cClickFitGaussPol") == NULL ){ + cClickFitNGaussPol = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click | clickFitGaussPol", degPol), 1400, 1200); + }else{ + delete gROOT->FindObjectAny("cClickFitNGaussPol") ; + cClickFitNGaussPol = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click | clickFitGaussPol", degPol), 1400, 1200); + } + cClickFitNGaussPol->Divide(1, 4); + for(int i = 1; i <= 4 ; i++){ + cClickFitNGaussPol->cd(i)->SetGrid(0,0); + } + + if(! cClickFitNGaussPol->GetShowEventStatus() ) cClickFitNGaussPol->ToggleEventStatus(); + cClickFitNGaussPol->cd(1); + + + hist->GetListOfFunctions()->Clear(); + ScaleAndDrawHist(hist, 0, 0); + TH1F* hspec = (TH1F*) hist->Clone(); + hspec->Sumw2(); + cClickFitNGaussPol->Update(); + cClickFitNGaussPol->Draw(); + + TLatex helpMsg; + helpMsg.SetNDC(); + helpMsg.SetTextFont(82); + helpMsg.SetTextSize(0.06); + helpMsg.SetTextColor(kRed); + helpMsg.DrawLatex(0.15, 0.8, "Click for peaks: (Double-click / x / Ctrl to end) "); + + printf("--------- Click the histogram for peaks: (Double-click / x / Ctrl to end) \n"); + nClick = 0; + xPeakList.clear(); + yPeakList.clear(); + markerColor = 2; + markerStyle = 23; + peakFlag = 1; + cClickFitNGaussPol->cd(1)->SetCrosshair(1); + cClickFitNGaussPol->cd(1)->AddExec("ex", "Clicked()"); + tempHist = hist; + TObject * obj ; + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + + if( degPol >= 0 ){ + printf("--------- Click the histogram for Background: (Double-click / x / Ctrl to end) \n"); + printf(" * when no input, program will estimate \n"); + + cClickFitNGaussPol->cd(1)->Clear(); + hist->Draw(); + helpMsg.SetTextColor(markerBGColor); + helpMsg.DrawLatex(0.15, 0.8, "Click for BG: (Double-click / x / Ctrl to end) "); + helpMsg.DrawLatex(0.15, 0.75, "* when no input, program will estimate"); + cClickFitNGaussPol->Update(); + cClickFitNGaussPol->Draw(); + + nClick = 0; + xBGList.clear(); + yBGList.clear(); + markerColor = markerBGColor; + markerStyle = 33; + peakFlag = 0; + cClickFitNGaussPol->cd(1)->AddExec("ex", "Clicked()"); + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + + } + + cClickFitNGaussPol->cd(1)->DeleteExec("ex"); + cClickFitNGaussPol->cd(1)->SetCrosshair(0); + cClickFitNGaussPol->cd(1)->Clear(); + hist->Draw(); + + tempHist = NULL; + + cClickFitNGaussPol->Update(); + cClickFitNGaussPol->Draw(); + + nPols = degPol; + double xMin = hspec->GetXaxis()->GetXmin(); + double xMax = hspec->GetXaxis()->GetXmax(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++) polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + TF1 *bg = new TF1("bg", polExp.Data(), xMin, xMax); + bg->SetNpx(1000); + bg->SetLineColor(markerBGColor); + bg->SetLineWidth(1); + if( xBGList.size() > 0 ) { + printf("---------------- fit the BG with Pol-%d \n", nPols); + TGraph * gBG = new TGraph((int) xBGList.size(), &xBGList[0], &yBGList[0]); + for( int i = 0; i < degPol + 1; i++) bg->SetParameter(i, gRandom->Rndm()/TMath::Power(10, 3*i)); + gBG->Fit("bg", "Rq"); + bg->Draw("same"); + //printf("--------------- Subtracting the BG \n"); + //hspec->Add(bg, -1); + }else{ + for( int i = 0; i < nPols+1; i++) bg->SetParameter(i, 0.); + } + + nPeaks = (int) xPeakList.size(); + if( sigmaMax == 0 ){ + printf("------------- Estimate sigma for each peak \n"); + sigma.clear(); + int binMin = hist->FindBin(xMin); + int binMax = hist->FindBin(xMax); + for( int i = 0; i < nPeaks ; i++){ + int b0 = hist->FindBin(xPeakList[i]); + double estBG = bg->Eval(xPeakList[i]); + double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.; + //---- backward search, stop when + for( int b = b0-1 ; b > binMin ; b-- ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMin = xPeakList[i] - hist->GetBinCenter(b); + break; + } + } + //---- forward search, stop when + for( int b = b0+1 ; b < binMax ; b++ ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMax = hist->GetBinCenter(b) - xPeakList[i]; + break; + } + } + + double temp = TMath::Min(sMin, sMax); + /// When there are multiple peaks closely packed : + if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back(); + sigma.push_back(temp); + + printf("%2d | x : %8.2f | sigma(est) %f \n", i, xPeakList[i], sigma[i]); + + } + + //---- use the mean of the sigma + double sigma0 = 0; + for( int i = 0; i < nPeaks ; i++) sigma0 += sigma[i]; + sigma0 = sigma0/(nPeaks+1); + for( int i = 0; i < nPeaks ; i++) sigma[i] = sigma0; + printf("========== use the mean sigma : %f \n", sigma0); + }else if( sigmaMax < 0 ){ + printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax)); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax)); + }else if( sigmaMax > 0 ){ + printf("========== use user input sigma : %f/2. \n", sigmaMax/2.); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.); + } + + printf("-------------- Fit the spectrum with %d-Gauss + Pol-%d\n", nPeaks, nPols ); + cClickFitNGaussPol->cd(2); + hspec->Draw("hist"); + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = yPeakList[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = xPeakList[i]; + if( sigmaMax == 0){ + para[3*i+2] = sigma[i]; + }else if(sigmaMax < 0 ){ + para[3*i+2] = abs(sigmaMax); + }else if(sigmaMax > 0 ){ + para[3*i+2] = sigmaMax/2.; + } + } + for(int i = 0; i < nPols+1 ; i++){ + para[3*nPeaks+i] = bg->GetParameter(i); + } + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //limit parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + if( meanRange <= 0 ) { + double dE1, dE2; + if( i == 0 ){ + dE2 = xPeakList[i+1] - xPeakList[i]; + dE1 = dE2; + }else if ( i == nPeaks-1 ){ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = dE1; + }else{ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = xPeakList[i+1] - xPeakList[i]; + } + fit->SetParLimits(3*i+1, xPeakList[i] - dE1 , xPeakList[i] + dE2 ); + }else{ + fit->SetParLimits(3*i+1, xPeakList[i] - meanRange/2. , xPeakList[i] + meanRange/2. ); + } + if( sigmaMax== 0 ) fit->SetParLimits(3*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma + if( sigmaMax < 0 ) fit->FixParameter(3*i+2, abs(sigmaMax)); + if( sigmaMax > 0 ) fit->SetParLimits(3*i+2, 0, sigmaMax); + } + + hspec->Fit("fit", "Rq"); + fit->Draw("same"); + + //=========== get the polynomial part + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + double bw = hspec->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + TF1 ** gFit = new TF1 *[nPeaks]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + if( degPol >= 0 ){ + TF1 *bg2 = new TF1("bg", polExp.Data(), xMin, xMax); + bg2->SetNpx(1000); + bg2->SetLineColor(markerBGColor); + bg2->SetLineWidth(1); + for( int i = 0; i < nPols + 1; i++){ + bg2->SetParameter(i, paraA[3*nPeaks+i]); + } + bg2->Draw("same"); + } + + cClickFitNGaussPol->cd(3); + PlotResidual(hspec, fit); + + cClickFitNGaussPol->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } + +} + + + +void clickFitNGaussPolSub(TH1F * hist, int degPol, double sigmaMax = 0, double meanRange = 0){ + + printf("=========================================================\n"); + printf("= fit n-Gauss + Pol-%d BG using mouse click (method-2) =\n", degPol ); + printf("==========================================================\n"); + + recentFitMethod = "clickFitNGaussPolSub"; + + gStyle->SetOptStat(""); + gStyle->SetOptFit(0); + + TCanvas * cClickFitNGaussPolsub = NULL; + if( gROOT->FindObjectAny("cClickFitGaussPolsub") == NULL ){ + cClickFitNGaussPolsub = new TCanvas("cClickFitNGaussPol", Form("fit Gauss & Pol-%d by mouse click (sub) | clickFitGaussPolsub", degPol), 1400, 1200); + }else{ + delete gROOT->FindObjectAny("cClickFitNGaussPolsub") ; + cClickFitNGaussPolsub = new TCanvas("cClickFitNGaussPolsub", Form("fit Gauss & Pol-%d by mouse click (sub) | clickFitGaussPolsub", degPol), 1400, 1200); + } + cClickFitNGaussPolsub->Divide(1, 4); + for(int i = 1; i <= 4 ; i++){ + cClickFitNGaussPolsub->cd(i)->SetGrid(0,0); + } + + if(! cClickFitNGaussPolsub->GetShowEventStatus() ) cClickFitNGaussPolsub->ToggleEventStatus(); + cClickFitNGaussPolsub->cd(1); + + + hist->GetListOfFunctions()->Clear(); + ScaleAndDrawHist(hist, 0, 0); + TH1F* hspec = (TH1F*) hist->Clone(); + hspec->Sumw2(); + hspec->SetName("subtraced"); + cClickFitNGaussPolsub->Update(); + cClickFitNGaussPolsub->Draw(); + + TLatex helpMsg; + helpMsg.SetNDC(); + helpMsg.SetTextFont(82); + helpMsg.SetTextSize(0.06); + helpMsg.SetTextColor(kRed); + helpMsg.DrawLatex(0.15, 0.8, "Click for peaks: (Double-click / x / Ctrl to end) "); + + + printf("--------- Click the histogram for peaks: (Double-click / x / Ctrl to end) \n"); + nClick = 0; + xPeakList.clear(); + yPeakList.clear(); + markerColor = 2; + markerStyle = 23; + peakFlag = 1; + cClickFitNGaussPolsub->cd(1)->SetCrosshair(1); + cClickFitNGaussPolsub->cd(1)->AddExec("ex", "Clicked()"); + tempHist = hist; + TObject * obj ; + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + + if( degPol >= 0 ){ + printf("--------- Click the histogram for Background: (Double-click / x / Ctrl to end) \n"); + + cClickFitNGaussPolsub->cd(1)->Clear(); + hist->Draw(); + helpMsg.SetTextColor(markerBGColor); + helpMsg.DrawLatex(0.15, 0.8, "Click for BG: (Double-click / x / Ctrl to end) "); + cClickFitNGaussPolsub->Update(); + cClickFitNGaussPolsub->Draw(); + + nClick = 0; + xBGList.clear(); + yBGList.clear(); + markerColor = markerBGColor; + markerStyle = 33; + peakFlag = 0; + do{ + obj = gPad->WaitPrimitive(); + if( obj == NULL ) break; + }while( obj != NULL); + } + cClickFitNGaussPolsub->cd(1)->DeleteExec("ex"); + cClickFitNGaussPolsub->cd(1)->SetCrosshair(0); + cClickFitNGaussPolsub->cd(1)->Clear(); + hist->Draw(); + + tempHist = NULL; + + if( xBGList.size() == 0 ) helpMsg.DrawLatex(0.15, 0.75, " No BG defined ! fitting could be problematics. "); + + cClickFitNGaussPolsub->Update(); + cClickFitNGaussPolsub->Draw(); + + nPols = degPol; + double xMin = hspec->GetXaxis()->GetXmin(); + double xMax = hspec->GetXaxis()->GetXmax(); + + TString polExp = "[0]"; + for( int i = 1; i < degPol + 1; i++) polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + TF1 *bg = new TF1("bg", polExp.Data(), xMin, xMax); + bg->SetNpx(1000); + bg->SetLineColor(markerBGColor); + bg->SetLineWidth(1); + if( xBGList.size() > 0 ) { + printf("---------------- fit the BG with Pol-%d \n", nPols); + + TGraph * gBG = new TGraph((int) xBGList.size(), &xBGList[0], &yBGList[0]); + + for( int i = 0; i < degPol + 1; i++) bg->SetParameter(i, gRandom->Rndm()/TMath::Power(10, 3*i)); + + gBG->Fit("bg", "Rq"); + bg->Draw("same"); + + printf(">>>> fitted parameters of pol.\n"); + for( int i = 0; i< degPol + 1; i++){ + printf("%2i | %.f\n", i, bg->GetParameter(i)); + } + + printf("--------------- Subtracting the BG \n"); + hspec->Add(bg, -1); + }else{ + for( int i = 0; i < nPols+1; i++) bg->SetParameter(i, 0.); + } + + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.06); + + double chi2BG = bg->GetChisquare(); + int ndfBG = bg->GetNDF(); + + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2BG/ndfBG)); + + //=========== get the polynomial BG + const Double_t* paraAt = bg->GetParameters(); + const Double_t* paraEt = bg->GetParErrors(); + + for( int i = 0; i < degPol + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.2f(%8.2f)\n", Form("p%d", i), paraAt[i], paraEt[i])); + } + + nPeaks = (int) xPeakList.size(); + if( sigmaMax == 0 ){ + printf("------------- Estimate sigma for each peak \n"); + sigma.clear(); + int binMin = hist->FindBin(xMin); + int binMax = hist->FindBin(xMax); + + for( int i = 0; i < nPeaks ; i++){ + int b0 = hist->FindBin(xPeakList[i]); + double estBG = bg->Eval(xPeakList[i]); + double sMin = (xMax-xMin)/5., sMax = (xMax-xMin)/5.; + //---- backward search, stop when + for( int b = b0-1 ; b > binMin ; b-- ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMin = xPeakList[i] - hist->GetBinCenter(b); + break; + } + } + //---- forward search, stop when + for( int b = b0+1 ; b < binMax ; b++ ){ + double y = hist->GetBinContent(b); + double x = hist->GetBinCenter(b); + if( y - (bg->Eval(x)) < (yPeakList[i]-estBG)/2. ) { + sMax = hist->GetBinCenter(b) - xPeakList[i]; + break; + } + } + + double temp = TMath::Min(sMin, sMax); + /// When there are multiple peaks closely packed : + if( i > 0 && temp > 2.5 * sigma.back() ) temp = sigma.back(); + sigma.push_back(temp); + + printf("%2d | x : %8.2f | sigma(est) %f \n", i, xPeakList[i], sigma[i]); + + } + }else if( sigmaMax < 0 ){ + printf("========== use user input sigma : %f (fixed)\n", abs(sigmaMax)); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(abs(sigmaMax)); + }else if( sigmaMax > 0 ){ + printf("========== use user input sigma : %f/2. \n", sigmaMax/2.); + sigma.clear(); + for( int i = 0; i < nPeaks ; i++) sigma.push_back(sigmaMax/2.); + } + + printf("-------------- Fit the subtracted spectrum with %d-Gauss\n", nPeaks ); + cClickFitNGaussPolsub->cd(2); + hspec->Draw("hist"); + + int nPar = 3 * nPeaks; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = yPeakList[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = xPeakList[i]; + if( sigmaMax == 0){ + para[3*i+2] = sigma[i]; + }else if(sigmaMax < 0 ){ + para[3*i+2] = abs(sigmaMax); + }else if(sigmaMax > 0 ){ + para[3*i+2] = sigmaMax/2.; + } + } + + TF1 * fit = new TF1("fit", nGauss, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //limit parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + + if( meanRange <= 0 ) { + double dE1, dE2; + if( i == 0 ){ + dE2 = xPeakList[i+1] - xPeakList[i]; + dE1 = dE2; + }else if ( i == nPeaks-1 ){ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = dE1; + }else{ + dE1 = xPeakList[i] - xPeakList[i-1]; + dE2 = xPeakList[i+1] - xPeakList[i]; + } + fit->SetParLimits(3*i+1, xPeakList[i] - dE1 , xPeakList[i] + dE2 ); + }else{ + fit->SetParLimits(3*i+1, xPeakList[i] - meanRange/2. , xPeakList[i] + meanRange/2. ); + } + if( sigmaMax== 0 ) fit->SetParLimits(3*i+2, 0, 1.5*sigma[i]); // add 50% margin of sigma + if( sigmaMax < 0 ) fit->FixParameter(3*i+2, abs(sigmaMax)); + if( sigmaMax > 0 ) fit->SetParLimits(3*i+2, 0, sigmaMax); + } + + hspec->Fit("fit", "Rq"); + fit->Draw("same"); + + //=========== get the fit parameters + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + double bw = hspec->GetBinWidth(1); + + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + TF1 ** gFit = new TF1 *[nPeaks]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + cClickFitNGaussPolsub->cd(3); + PlotResidual(hspec, fit); + + cClickFitNGaussPolsub->cd(4); + + text.SetTextSize(0.05); + text.SetTextColor(2); + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + BestFitCount.clear(); + BestFitMean.clear(); + BestFitSigma.clear(); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + + BestFitCount.push_back(paraA[3*i]/ bw); + BestFitMean.push_back(paraA[3*i+1]); + BestFitSigma.push_back(paraA[3*i+2]); + } +} + + +//######################################## +//#### fit N-Gauss with pol-1 BG at fixed for < 0 +//######################################## +void fitSpecial(TH1F * hist, TString fitFile = "AutoFit_para.txt"){ + + printf("================================================================\n"); + printf("================ Special fit for h074_82Kr =====================\n"); + printf(" * need the file input \n"); + printf("================================================================\n"); + + + bool isParaRead = loadFitParameters(fitFile); + if( !isParaRead ) { + printf("Please provide a valid input file\n"); + return; + } + + gStyle->SetOptStat(""); + gStyle->SetOptFit(0); + nPeaks = energy.size(); + nPols = 1; + + TCanvas * cFitSpecial = NewCanvas("cFitSpecial", "Fitting for h074_82Kr", 1, 3, 800, 300); + //if(! cFitSpecial->GetShowEventStatus() ) cFitSpecial->ToggleEventStatus(); + cFitSpecial->cd(1); + + ScaleAndDrawHist(hist, 0, 0); + + double xMin = hist->GetXaxis()->GetXmin(); + double xMax = hist->GetXaxis()->GetXmax(); + int xBin = hist->GetXaxis()->GetNbins(); + + printf("============= find the pol-1 background ..... \n"); + + TF1 * fitpol1 = new TF1("fitpol1", "pol1", xMin, -0.3); + fitpol1->SetParameter(0, gRandom->Rndm()); + fitpol1->SetParameter(1, gRandom->Rndm()); + + hist->Fit("fitpol1", "Rq"); + + double x0 = fitpol1->GetParameter(0); + double x1 = fitpol1->GetParameter(1); + + + int nPar = 3 * nPeaks + nPols + 1; + double * para = new double[nPar]; + for(int i = 0; i < nPeaks ; i++){ + para[3*i+0] = height[i] * 0.05 * TMath::Sqrt(TMath::TwoPi()); + para[3*i+1] = energy[i]; + para[3*i+2] = sigma[i]/2.; + } + + para[3*nPeaks+0] = x0; + para[3*nPeaks+1] = x1; + + + TF1 * fit = new TF1("fit", nGaussPol, xMin, xMax, nPar ); + fit->SetLineWidth(3); + fit->SetLineColor(1); + fit->SetNpx(1000); + fit->SetParameters(para); + + //fixing parameters + for( int i = 0; i < nPeaks; i++){ + fit->SetParLimits(3*i , 0, 1e9); + if( energyFlag[i] == 1 ) { + fit->FixParameter(3*i+1, energy[i]); + }else{ + fit->SetParLimits(3*i+1, lowE[i], highE[i]); + } + if( sigmaFlag[i] == 1 ) { + fit->FixParameter(3*i+2, sigma[i]); + }else{ + fit->SetParLimits(3*i+2, 0, sigma[i]); + } + } + + fit->FixParameter(3*nPeaks + 0, x0); + fit->FixParameter(3*nPeaks + 1, x1); + + hist->Fit("fit", "Rq"); + + //=========== get the polynomial part + const Double_t* paraA = fit->GetParameters(); + const Double_t* paraE = fit->GetParErrors(); + + TString polExp = "[0]"; + for( int i = 1; i < nPols + 1; i++){ + polExp += Form("+[%d]*TMath::Power(x,%d)", i, i ); + } + + TF1 * bg = new TF1("bg", polExp.Data(), xMin, xMax); + for( int i = 0; i < nPols + 1; i++){ + bg->SetParameter(i, paraA[3*nPeaks+i]); + } + bg->SetNpx(1000); + + for( int i = 0; i < nPols + 1; i++){ + printf("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i]); + } + printf("====================================\n"); + + cFitSpecial->cd(1); + bg->Draw("same"); + + TLatex text; + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.04); + text.SetTextColor(1); + + for( int i = 0; i < nPols + 1; i++){ + text.DrawLatex(0.6, 0.85 - 0.05*i, Form("%4s : %8.4e(%8.4e)\n", Form("p%d", i), paraA[3*nPeaks+i], paraE[3*nPeaks+i])); + } + + double chi2 = fit->GetChisquare(); + int ndf = fit->GetNDF(); + text.SetTextSize(0.06); + text.DrawLatex(0.15, 0.8, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf)); + + + //GoodnessofFit(specS, fit); + + double bw = hist->GetBinWidth(1); + + for(int i = 0; i < nPeaks ; i++){ + printf(" %2d , count: %8.0f(%3.0f), mean: %8.4f(%8.4f), sigma: %8.4f(%8.4f) \n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2]); + } + printf("\n"); + + + const int nn = nPeaks; + TF1 ** gFit = new TF1 *[nn]; + for( int i = 0; i < nPeaks; i++){ + gFit[i] = new TF1(Form("gFit%d", i), "[0] * TMath::Gaus(x,[1],[2], 1)", xMin, xMax); + gFit[i]->SetParameter(0, paraA[3*i]); + gFit[i]->SetParameter(1, paraA[3*i+1]); + gFit[i]->SetParameter(2, paraA[3*i+2]); + gFit[i]->SetLineColor(i+1); + gFit[i]->SetNpx(1000); + gFit[i]->SetLineWidth(1); + gFit[i]->Draw("same"); + } + + cFitSpecial->Update(); + + cFitSpecial->cd(2); + PlotResidual(hist, fit); + + cFitSpecial->cd(3); + + text.SetNDC(); + text.SetTextFont(82); + text.SetTextSize(0.05); + text.SetTextColor(2); + + text.DrawLatex(0.1, 0.9, Form(" %13s, %18s, %18s", "count", "mean", "sigma")); + + for( int i = 0; i < nPeaks; i++){ + text.DrawLatex(0.1, 0.8-0.05*i, Form(" %2d, %8.0f(%3.0f), %8.4f(%8.4f), %8.4f(%8.4f)\n", + i, + paraA[3*i] / bw, paraE[3*i] /bw, + paraA[3*i+1], paraE[3*i+1], + paraA[3*i+2], paraE[3*i+2])); + } + +} + +// } + +#endif diff --git a/AutoFit_para.txt b/AutoFit_para.txt new file mode 100644 index 0000000..5528e65 --- /dev/null +++ b/AutoFit_para.txt @@ -0,0 +1,8 @@ +# for n-Gauss fit, can use "#", or "//" to comment out whole line +# peak low high fixed? sigma_Max fixed? hight +0.000 -0.1 0.1 1 0.08 0 20 +0.854 0.8 0.9 1 0.08 0 20 +1.363 1.3 1.4 1 0.08 0 10 +1.561 1.5 1.6 1 0.09 0 10 +2.005 1.9 2.1 1 0.09 0 20 +2.809 2.7 2.85 0 0.09 0 10 \ No newline at end of file diff --git a/SaveTH1IntoText.C b/SaveTH1IntoText.C new file mode 100644 index 0000000..1494e06 --- /dev/null +++ b/SaveTH1IntoText.C @@ -0,0 +1,24 @@ +#include "TFile.h" +#include "TTree.h" +#include "TH1F.h" + +void SaveTH1IntoText(TH1F * hist, TString fileName){ + + int nBin = hist->GetNbinsX(); + + + FILE * paraOut; + paraOut = fopen (fileName.Data(), "w+"); + printf("=========== save histogram to %s \n", fileName.Data()); + + for( int i = 1; i <= nBin ; i++){ + + double x = hist->GetBinCenter(i); + double y = hist->GetBinContent(i); + + fprintf(paraOut, "%14.8f\t%9.6f\n", x, y); + } + + fflush(paraOut); + fclose(paraOut); +} diff --git a/caliPeaks.cpp b/caliPeaks.cpp new file mode 100644 index 0000000..5c78fd5 --- /dev/null +++ b/caliPeaks.cpp @@ -0,0 +1,220 @@ +#include "TFile.h" +#include "TTree.h" +#include "TH1F.h" +#include "TStyle.h" +#include "TH2F.h" +#include "TCanvas.h" +#include "TTreeReader.h" + +#include "AutoFit.C" +#include "SaveTH1IntoText.C" + +// mod = 0, 1, 2 +// row = 1, 2, 3 +// pstrip = 0 - 127 + +#define MAXMOD 3 +#define MAXROW 3 +#define MAXPStrip 128 + +const double ExRange[2] = {-0.5, 3.5}; +const double ExRel = 0.05; //MeV + +#define SAVETREE false + +//^==================================== +TH1F * haha[MAXMOD][MAXROW][MAXPStrip]; + +TH2F * dudu[MAXMOD][MAXROW]; + +TH1F * hEx; + +//^==================================== +void caliPeaks(){ + + + /// Load peaks + std::vector peaks[MAXMOD][MAXROW][MAXPStrip]; + + FILE * inList = fopen("fit_result.txt", "r"); + + char buffer[256]; + + while( fgets(buffer, sizeof(buffer), inList) ){ + + std::vector haha = SplitStrAF(buffer, " "); + + // printf("%zu| %s", haha.size(), buffer); + int m = atoi(haha[0].c_str()); + int r = atoi(haha[1].c_str()) - 1; + int p = atoi(haha[2].c_str()); + + for( int i = 3; i < (int) haha.size(); i++ ) peaks[m][r][p].push_back(atof(haha[i].c_str())); + } + + fclose(inList); + + // for( int i = 0; i < MAXMOD; i++ ){ + // for( int j = 0 ; j < MAXROW; j++ ){ + // for( int p = 0 ; p < MAXPStrip; p++){ + // if( peaks[i][j][p].size() == 0 ) continue; + // printf("%d %d %3d | ", i, j+1, p); + // for( size_t h = 0; h < peaks[i][j][p].size(); h++) printf(" %6.3f", peaks[i][j][p][h]); + // printf("\n"); + // } + // } + // } + + TFile * file = new TFile("R9-172_summed_tree.root"); + + TTree * tree = (TTree *) file->Get("rxtree"); + + const int nBin = (ExRange[1] - ExRange[0])/ExRel; + + for( int i = 0; i < MAXMOD; i++ ){ + for( int j = 0; j < MAXROW; j++ ){ + + dudu[i][j] = new TH2F(Form("dudu_%d%d", i, j+1), Form("%d-%d; pStrip; Ex", i, j+1), 128, 0, 128, nBin, ExRange[0], ExRange[1]); + + for( int p = 0; p < MAXPStrip; p ++ ){ + haha[i][j][p] = new TH1F(Form("h_%d%d_%d", i, j+1, p), Form("%d-%d-%d", i, j+1, p), nBin, ExRange[0], ExRange[1]); + } + } + } + + hEx = new TH1F("hEx", "Ex", 10*nBin, ExRange[0], ExRange[1]); + + TTreeReader reader(tree); + + TTreeReaderValue Ex = {reader, "Ex"}; + TTreeReaderValue mod = {reader, "mod"}; + TTreeReaderValue row = {reader, "row"}; + TTreeReaderValue pStrip = {reader, "pstrip"}; + TTreeReaderValue ebis = {reader, "ebis_on"}; + TTreeReaderValue ThetaCM = {reader, "ThetaCM"}; + TTreeReaderValue Zmeasured = {reader, "Zmeasured"}; + + TFile * newFile = nullptr; + TTree * newTree = nullptr; + + Double_t Ex_org; + Double_t Ex_C; + Int_t mod_C; + Int_t row_C; + Int_t pStrip_C; + Double_t z; + Double_t thetaCM; + + if( SAVETREE ){ + TFile * newFile = new TFile("Cali_132Sn.root", "recreate"); + TTree * newTree = new TTree ("tree", "Calibrated Ex"); + + newTree->Branch("Ex_org", &Ex_org, "Ex_org/D"); + newTree->Branch("mod", &mod_C, "mod/I"); + newTree->Branch("row", &row_C, "row/I"); + newTree->Branch("pstrip", &pStrip_C, "pstrip/I"); + newTree->Branch("Ex", &Ex_C, "Ex/D"); + newTree->Branch("z", &z, "z/D"); + newTree->Branch("thetaCM", &thetaCM, "thetaCM/D"); + } + + const double y1 = 0.000; + const double ya = 0.854; + const double y2 = 2.005; + + while( reader.Next() ){ + + if( *ebis == 0 ) continue; + if( *row == 0 ) continue; + + int m = *mod; + int r = *row -1; + int p = *pStrip; + + // printf("%d %d %d | %f\n", m, r, p, *Ex/1000.); + + if( peaks[m][r][p].size() == 0 ) continue; + + double slope = 1; + double offset = 0; + + if( peaks[m][r][p].size() == 1 ) { + if( peaks[m][r][p][0] < 1 ){ + offset = -peaks[m][r][p][0]; + }else{ + offset = y2 - peaks[m][r][p][0]; + } + } + if( peaks[m][r][p].size() >= 2 ){ + double x1 = peaks[m][r][p][0]; + double x2 = peaks[m][r][p].back(); + + if( x1 < 0.2 ) { + if( x2 < 1.8 ) { + slope = (ya-y1)/(x2-x1); + }else{ + slope = (y2-y1)/(x2-x1); + } + offset = y1 - slope * x1; + } + if( x1 >= 0.2 ) { + // printf("%d %d %d| x1 : %f, ya : %f \n", m, r, p, x1, ya); + slope = (y2-ya)/(x2-x1); + offset = ya - slope * x1; + } + } + + Ex_org = *Ex/1000.; + mod_C = m; + row_C = r; + pStrip_C = p; + z = *Zmeasured; + thetaCM = *ThetaCM; + + Ex_C = slope * (*Ex/1000.) + offset; + + haha[m][r][p]->Fill(Ex_C); + + dudu[m][r]->Fill(p, Ex_C); + + hEx->Fill( *Ex/1000); + + if( SAVETREE ) { + newFile->cd(); + newTree->Fill(); + } + + } + + if( SAVETREE ) { + newFile->cd(); + newTree->Write(); + printf("============== %llu\n", newTree->GetEntries()); + newFile->Close(); + } + + gStyle->SetOptStat(""); + + // TCanvas * canvas = new TCanvas("canvas", "canvas", 2100, 1500); + // // haha[0][2][10]->Draw(); + // canvas->Divide(3,3); + + // canvas->cd(1); dudu[0][0]->Draw("colz"); + // canvas->cd(2); dudu[0][1]->Draw("colz"); + // canvas->cd(3); dudu[0][2]->Draw("colz"); + + // canvas->cd(4); dudu[1][0]->Draw("colz"); + // canvas->cd(5); dudu[1][1]->Draw("colz"); + // canvas->cd(6); dudu[1][2]->Draw("colz"); + + // canvas->cd(7); dudu[2][0]->Draw("colz"); + // canvas->cd(8); dudu[2][1]->Draw("colz"); + // canvas->cd(9); dudu[2][2]->Draw("colz"); + + TCanvas * canvasEx = new TCanvas("canvasEx", "canvasEx", 800, 400); + canvasEx->cd(); + hEx->Draw(); + + + +} \ No newline at end of file diff --git a/ex_sub_50keV.txt b/ex_sub_50keV.txt new file mode 100644 index 0000000..fe0c260 --- /dev/null +++ b/ex_sub_50keV.txt @@ -0,0 +1,800 @@ + -0.49750000 11.810265 + -0.49250000 20.502277 + -0.48750000 10.194973 + -0.48250000 14.888348 + -0.47750000 0.582401 + -0.47250000 7.277138 + -0.46750000 18.972553 + -0.46250000 20.668644 + -0.45750000 0.365417 + -0.45250000 11.062874 + -0.44750000 8.761005 + -0.44250000 12.459820 + -0.43750000 -6.840687 + -0.43250000 -0.140518 + -0.42750000 8.560337 + -0.42250000 17.261868 + -0.41750000 5.964081 + -0.41250000 -1.333031 + -0.40750000 7.370541 + -0.40250000 15.074791 + -0.39750000 -4.220276 + -0.39250000 5.485332 + -0.38750000 18.191624 + -0.38250000 8.898594 + -0.37750000 6.606243 + -0.37250000 24.314571 + -0.36750000 4.023582 + -0.36250000 7.733269 + -0.35750000 -4.556362 + -0.35250000 0.154686 + -0.34750000 17.866413 + -0.34250000 -5.421177 + -0.33750000 1.291908 + -0.33250000 7.005676 + -0.32750000 19.720123 + -0.32250000 0.435249 + -0.31750000 18.151058 + -0.31250000 11.867542 + -0.30750000 2.584709 + -0.30250000 9.302555 + -0.29750000 16.021080 + -0.29250000 13.740284 + -0.28750000 4.460171 + -0.28250000 14.180733 + -0.27750000 -7.098022 + -0.27250000 2.623901 + -0.26750000 4.346508 + -0.26250000 0.069790 + -0.25750000 5.793755 + -0.25250000 1.518394 + -0.24750000 8.243717 + -0.24250000 12.969723 + -0.23750000 31.696404 + -0.23250000 -0.576233 + -0.22750000 3.151806 + -0.22250000 6.880527 + -0.21750000 14.609932 + -0.21250000 -1.659988 + -0.20750000 12.070770 + -0.20250000 26.802212 + -0.19750000 22.534332 + -0.19250000 21.267132 + -0.18750000 42.000610 + -0.18250000 32.734772 + -0.17750000 44.469612 + -0.17250000 62.205132 + -0.16750000 46.941330 + -0.16250000 60.678207 + -0.15750000 69.415764 + -0.15250000 77.154007 + -0.14750000 67.892921 + -0.14250000 108.632523 + -0.13750000 113.372795 + -0.13250000 132.113754 + -0.12750000 185.855392 + -0.12250000 169.597702 + -0.11750000 189.340698 + -0.11250000 238.084381 + -0.10750000 267.828735 + -0.10250000 292.573761 + -0.09750000 358.319489 + -0.09250000 408.065887 + -0.08750000 396.812958 + -0.08250000 467.560699 + -0.07750000 515.309143 + -0.07250000 494.058258 + -0.06750000 580.808044 + -0.06250000 609.558533 + -0.05750000 688.309692 + -0.05250000 684.061523 + -0.04750000 704.814026 + -0.04250000 823.567200 + -0.03750000 812.321106 + -0.03250000 822.075623 + -0.02750000 849.830872 + -0.02250000 850.586792 + -0.01750000 868.343384 + -0.01250000 951.100647 + -0.00750000 941.858582 + -0.00250000 975.617188 + 0.00250000 956.376526 + 0.00750000 909.136536 + 0.01250000 894.897217 + 0.01750000 936.658569 + 0.02250000 894.420532 + 0.02750000 873.183289 + 0.03250000 847.946655 + 0.03750000 821.710754 + 0.04250000 774.475464 + 0.04750000 739.240906 + 0.05250000 718.007019 + 0.05750000 747.773804 + 0.06250000 701.541260 + 0.06750000 677.309387 + 0.07250000 625.078247 + 0.07750000 586.847778 + 0.08250000 567.617920 + 0.08750000 532.388794 + 0.09250000 470.160339 + 0.09750000 480.932556 + 0.10250000 460.705475 + 0.10750000 428.479065 + 0.11250000 363.253326 + 0.11750000 345.028259 + 0.12250000 322.803894 + 0.12750000 345.580200 + 0.13250000 269.357178 + 0.13750000 268.134827 + 0.14250000 226.913177 + 0.14750000 219.692200 + 0.15250000 193.471893 + 0.15750000 165.252289 + 0.16250000 169.033356 + 0.16750000 136.815094 + 0.17250000 132.597504 + 0.17750000 118.380608 + 0.18250000 109.164391 + 0.18750000 98.948853 + 0.19250000 74.733994 + 0.19750000 70.519814 + 0.20250000 55.306313 + 0.20750000 46.093491 + 0.21250000 39.881348 + 0.21750000 65.669891 + 0.22250000 42.459106 + 0.22750000 52.249008 + 0.23250000 18.039581 + 0.23750000 31.830841 + 0.24250000 19.622780 + 0.24750000 32.415398 + 0.25250000 10.208694 + 0.25750000 4.002670 + 0.26250000 29.797325 + 0.26750000 25.592659 + 0.27250000 7.388680 + 0.27750000 43.185371 + 0.28250000 31.982750 + 0.28750000 14.780800 + 0.29250000 0.579536 + 0.29750000 24.378952 + 0.30250000 -3.820953 + 0.30750000 9.979820 + 0.31250000 9.781273 + 0.31750000 -6.416588 + 0.32250000 17.386223 + 0.32750000 -1.810287 + 0.33250000 7.993889 + 0.33750000 6.798744 + 0.34250000 8.604271 + 0.34750000 14.410484 + 0.35250000 1.217377 + 0.35750000 2.024948 + 0.36250000 8.833199 + 0.36750000 14.642128 + 0.37250000 9.451744 + 0.37750000 5.262032 + 0.38250000 4.072998 + 0.38750000 -14.115349 + 0.39250000 -5.303017 + 0.39750000 0.509995 + 0.40250000 7.323677 + 0.40750000 -6.861954 + 0.41250000 -7.046906 + 0.41750000 8.768829 + 0.42250000 -27.414764 + 0.42750000 2.402321 + 0.43250000 6.220093 + 0.43750000 -1.961464 + 0.44250000 5.857666 + 0.44750000 -2.322533 + 0.45250000 10.497955 + 0.45750000 -14.680878 + 0.46250000 18.140968 + 0.46750000 -3.036507 + 0.47250000 -7.213295 + 0.47750000 6.610588 + 0.48250000 -14.564850 + 0.48750000 -9.739601 + 0.49250000 9.086319 + 0.49750000 -18.087074 + 0.50250000 5.740211 + 0.50750000 4.568176 + 0.51250000 2.396820 + 0.51750000 -9.773857 + 0.52250000 -7.943855 + 0.52750000 -2.113174 + 0.53250000 8.718185 + 0.53750000 4.550232 + 0.54250000 -2.617050 + 0.54750000 12.216354 + 0.55250000 9.050438 + 0.55750000 -9.114807 + 0.56250000 3.720634 + 0.56750000 21.556755 + 0.57250000 -15.606445 + 0.57750000 2.231041 + 0.58250000 -13.930801 + 0.58750000 -14.091965 + 0.59250000 -2.252441 + 0.59750000 0.587753 + 0.60250000 4.428635 + 0.60750000 7.270195 + 0.61250000 23.112434 + 0.61750000 -0.044647 + 0.62250000 5.798950 + 0.62750000 0.643227 + 0.63250000 15.488182 + 0.63750000 9.333817 + 0.64250000 2.180138 + 0.64750000 8.027130 + 0.65250000 22.874809 + 0.65750000 21.723167 + 0.66250000 26.572205 + 0.66750000 40.421913 + 0.67250000 42.272308 + 0.67750000 40.123390 + 0.68250000 38.975143 + 0.68750000 67.827576 + 0.69250000 96.680687 + 0.69750000 80.534485 + 0.70250000 108.388954 + 0.70750000 132.244110 + 0.71250000 166.099945 + 0.71750000 201.956451 + 0.72250000 194.813660 + 0.72750000 251.671524 + 0.73250000 295.530090 + 0.73750000 281.389313 + 0.74250000 334.249207 + 0.74750000 366.109802 + 0.75250000 424.971069 + 0.75750000 416.833038 + 0.76250000 480.695679 + 0.76750000 502.558960 + 0.77250000 529.422974 + 0.77750000 561.287598 + 0.78250000 603.152954 + 0.78750000 650.019043 + 0.79250000 681.885742 + 0.79750000 666.753113 + 0.80250000 721.621155 + 0.80750000 724.489929 + 0.81250000 685.359375 + 0.81750000 745.229492 + 0.82250000 728.100220 + 0.82750000 751.971741 + 0.83250000 736.843872 + 0.83750000 813.716675 + 0.84250000 703.590210 + 0.84750000 786.464355 + 0.85250000 776.339233 + 0.85750000 738.214783 + 0.86250000 786.091003 + 0.86750000 757.967896 + 0.87250000 716.845520 + 0.87750000 702.723755 + 0.88250000 686.602722 + 0.88750000 694.482361 + 0.89250000 664.362671 + 0.89750000 612.243652 + 0.90250000 584.125305 + 0.90750000 618.007629 + 0.91250000 522.890625 + 0.91750000 568.774353 + 0.92250000 534.658752 + 0.92750000 482.543823 + 0.93250000 504.429565 + 0.93750000 417.315979 + 0.94250000 429.203064 + 0.94750000 392.090851 + 0.95250000 352.979309 + 0.95750000 356.868469 + 0.96250000 344.758270 + 0.96750000 277.648773 + 0.97250000 238.539948 + 0.97750000 250.431808 + 0.98250000 246.324341 + 0.98750000 211.217560 + 0.99250000 231.111465 + 0.99750000 190.006042 + 1.00250000 180.901306 + 1.00750000 157.797241 + 1.01250000 124.693863 + 1.01750000 130.591156 + 1.02250000 119.489136 + 1.02750000 114.387794 + 1.03250000 84.287125 + 1.03750000 86.187141 + 1.04250000 63.087837 + 1.04750000 68.989220 + 1.05250000 52.891273 + 1.05750000 60.794006 + 1.06250000 44.697426 + 1.06750000 45.601517 + 1.07250000 25.506294 + 1.07750000 27.411743 + 1.08250000 34.317879 + 1.08750000 31.224693 + 1.09250000 26.132187 + 1.09750000 26.040359 + 1.10250000 12.949211 + 1.10750000 -0.141251 + 1.11250000 -6.231041 + 1.11750000 17.679848 + 1.12250000 8.591423 + 1.12750000 11.503677 + 1.13250000 26.416603 + 1.13750000 -3.669785 + 1.14250000 -1.755493 + 1.14750000 4.159477 + 1.15250000 2.075127 + 1.15750000 0.991463 + 1.16250000 23.908470 + 1.16750000 2.826157 + 1.17250000 7.744530 + 1.17750000 19.663574 + 1.18250000 9.583305 + 1.18750000 30.503716 + 1.19250000 25.424805 + 1.19750000 31.346573 + 1.20250000 46.269020 + 1.20750000 34.192146 + 1.21250000 47.115952 + 1.21750000 81.040443 + 1.22250000 91.965607 + 1.22750000 83.891457 + 1.23250000 95.817978 + 1.23750000 115.745186 + 1.24250000 98.673073 + 1.24750000 147.601639 + 1.25250000 179.530884 + 1.25750000 161.460815 + 1.26250000 161.391418 + 1.26750000 199.322693 + 1.27250000 252.254669 + 1.27750000 246.187302 + 1.28250000 256.120636 + 1.28750000 284.054626 + 1.29250000 242.989319 + 1.29750000 306.924683 + 1.30250000 304.860718 + 1.30750000 309.797424 + 1.31250000 323.734863 + 1.31750000 308.672913 + 1.32250000 299.611694 + 1.32750000 313.551147 + 1.33250000 354.491241 + 1.33750000 330.432068 + 1.34250000 295.373535 + 1.34750000 309.315674 + 1.35250000 319.258545 + 1.35750000 350.202057 + 1.36250000 322.146240 + 1.36750000 348.091125 + 1.37250000 328.036682 + 1.37750000 328.982941 + 1.38250000 312.929840 + 1.38750000 301.877441 + 1.39250000 282.825714 + 1.39750000 294.774689 + 1.40250000 275.724304 + 1.40750000 287.674622 + 1.41250000 274.625641 + 1.41750000 283.577301 + 1.42250000 259.529663 + 1.42750000 250.482697 + 1.43250000 242.436417 + 1.43750000 249.390808 + 1.44250000 250.345886 + 1.44750000 248.301636 + 1.45250000 235.258087 + 1.45750000 233.215195 + 1.46250000 246.172989 + 1.46750000 206.131470 + 1.47250000 235.090607 + 1.47750000 203.050446 + 1.48250000 252.010956 + 1.48750000 162.972168 + 1.49250000 199.934021 + 1.49750000 213.896576 + 1.50250000 175.859818 + 1.50750000 193.823730 + 1.51250000 172.788315 + 1.51750000 162.753601 + 1.52250000 131.719543 + 1.52750000 149.686172 + 1.53250000 142.653473 + 1.53750000 137.621475 + 1.54250000 145.590149 + 1.54750000 124.559494 + 1.55250000 138.529526 + 1.55750000 117.500237 + 1.56250000 105.471619 + 1.56750000 121.443687 + 1.57250000 102.416443 + 1.57750000 115.389870 + 1.58250000 91.363976 + 1.58750000 83.338760 + 1.59250000 86.314232 + 1.59750000 88.290382 + 1.60250000 67.267204 + 1.60750000 70.244713 + 1.61250000 57.222900 + 1.61750000 51.201767 + 1.62250000 46.181313 + 1.62750000 80.161537 + 1.63250000 58.142441 + 1.63750000 59.124031 + 1.64250000 27.106293 + 1.64750000 51.089241 + 1.65250000 52.072861 + 1.65750000 62.057167 + 1.66250000 47.042152 + 1.66750000 33.027817 + 1.67250000 27.014160 + 1.67750000 9.001183 + 1.68250000 15.988884 + 1.68750000 18.977264 + 1.69250000 3.966331 + 1.69750000 21.956070 + 1.70250000 11.946495 + 1.70750000 21.937599 + 1.71250000 2.929375 + 1.71750000 0.921837 + 1.72250000 -0.085022 + 1.72750000 10.908798 + 1.73250000 9.903297 + 1.73750000 -0.101517 + 1.74250000 15.894341 + 1.74750000 9.890884 + 1.75250000 0.888100 + 1.75750000 5.886002 + 1.76250000 0.884575 + 1.76750000 -3.116165 + 1.77250000 0.883774 + 1.77750000 28.884392 + 1.78250000 7.885689 + 1.78750000 -9.112335 + 1.79250000 10.890327 + 1.79750000 32.893661 + 1.80250000 41.897682 + 1.80750000 25.902374 + 1.81250000 30.907753 + 1.81750000 51.913811 + 1.82250000 58.920547 + 1.82750000 70.927963 + 1.83250000 108.936058 + 1.83750000 129.944824 + 1.84250000 135.954285 + 1.84750000 176.964417 + 1.85250000 170.975235 + 1.85750000 182.986725 + 1.86250000 244.998901 + 1.86750000 291.011749 + 1.87250000 359.025269 + 1.87750000 370.039490 + 1.88250000 405.054382 + 1.88750000 480.069977 + 1.89250000 518.086243 + 1.89750000 551.103149 + 1.90250000 613.120728 + 1.90750000 691.139038 + 1.91250000 702.158020 + 1.91750000 767.177673 + 1.92250000 738.197998 + 1.92750000 777.218994 + 1.93250000 815.240723 + 1.93750000 810.263062 + 1.94250000 772.286133 + 1.94750000 796.309875 + 1.95250000 804.334290 + 1.95750000 879.359375 + 1.96250000 899.385132 + 1.96750000 867.411560 + 1.97250000 887.438721 + 1.97750000 896.466492 + 1.98250000 880.494995 + 1.98750000 938.524170 + 1.99250000 873.554016 + 1.99750000 948.584534 + 2.00250000 897.615723 + 2.00750000 874.647644 + 2.01250000 877.680176 + 2.01750000 814.713440 + 2.02250000 925.747375 + 2.02750000 789.781982 + 2.03250000 863.817261 + 2.03750000 751.853210 + 2.04250000 782.889893 + 2.04750000 734.927185 + 2.05250000 727.965210 + 2.05750000 747.003906 + 2.06250000 707.043213 + 2.06750000 599.083252 + 2.07250000 629.124023 + 2.07750000 597.165405 + 2.08250000 568.207520 + 2.08750000 564.250244 + 2.09250000 547.293701 + 2.09750000 495.337830 + 2.10250000 452.382629 + 2.10750000 422.428101 + 2.11250000 430.474243 + 2.11750000 429.521088 + 2.12250000 389.568604 + 2.12750000 380.616821 + 2.13250000 328.665710 + 2.13750000 294.715271 + 2.14250000 294.765503 + 2.14750000 275.816406 + 2.15250000 244.868011 + 2.15750000 238.920288 + 2.16250000 223.973251 + 2.16750000 214.026886 + 2.17250000 185.081207 + 2.17750000 162.136200 + 2.18250000 158.191879 + 2.18750000 152.248230 + 2.19250000 163.305267 + 2.19750000 136.362991 + 2.20250000 131.421387 + 2.20750000 97.480461 + 2.21250000 66.540215 + 2.21750000 78.600655 + 2.22250000 62.661766 + 2.22750000 63.723564 + 2.23250000 65.786034 + 2.23750000 42.849190 + 2.24250000 48.913025 + 2.24750000 54.977539 + 2.25250000 49.042732 + 2.25750000 35.108604 + 2.26250000 54.175156 + 2.26750000 19.242393 + 2.27250000 33.310303 + 2.27750000 22.378899 + 2.28250000 1.448166 + 2.28750000 38.518120 + 2.29250000 8.588753 + 2.29750000 6.660065 + 2.30250000 12.732056 + 2.30750000 -2.195274 + 2.31250000 -11.121925 + 2.31750000 9.952103 + 2.32250000 17.026817 + 2.32750000 -1.897797 + 2.33250000 -5.821724 + 2.33750000 11.255020 + 2.34250000 3.332451 + 2.34750000 -2.589439 + 2.35250000 -8.510651 + 2.35750000 -3.431183 + 2.36250000 -24.351036 + 2.36750000 0.729790 + 2.37250000 7.811302 + 2.37750000 -36.106514 + 2.38250000 -1.023643 + 2.38750000 -8.940102 + 2.39250000 7.144127 + 2.39750000 -7.770966 + 2.40250000 -10.685379 + 2.40750000 -3.599113 + 2.41250000 -5.512169 + 2.41750000 -16.424545 + 2.42250000 -7.336243 + 2.42750000 4.752739 + 2.43250000 -12.157593 + 2.43750000 -2.067253 + 2.44250000 -14.976227 + 2.44750000 -14.884521 + 2.45250000 -22.792137 + 2.45750000 -15.699074 + 2.46250000 -14.605331 + 2.46750000 -2.510910 + 2.47250000 -18.415810 + 2.47750000 -10.320030 + 2.48250000 -16.223564 + 2.48750000 -26.126427 + 2.49250000 -28.028603 + 2.49750000 -12.930099 + 2.50250000 -17.830925 + 2.50750000 -30.731064 + 2.51250000 -16.630524 + 2.51750000 -17.529305 + 2.52250000 -14.427399 + 2.52750000 -15.324821 + 2.53250000 -12.221565 + 2.53750000 -8.117622 + 2.54250000 -4.013008 + 2.54750000 -6.907707 + 2.55250000 -17.801727 + 2.55750000 -31.695068 + 2.56250000 -7.587738 + 2.56750000 -21.479713 + 2.57250000 -0.371017 + 2.57750000 -5.261642 + 2.58250000 -16.151588 + 2.58750000 -15.040848 + 2.59250000 6.070564 + 2.59750000 -8.817337 + 2.60250000 -21.704567 + 2.60750000 -7.591110 + 2.61250000 -13.476974 + 2.61750000 -23.362160 + 2.62250000 -14.246666 + 2.62750000 -15.130493 + 2.63250000 -1.013634 + 2.63750000 -19.896103 + 2.64250000 -5.777893 + 2.64750000 -1.658997 + 2.65250000 -14.539421 + 2.65750000 12.580826 + 2.66250000 9.701759 + 2.66750000 -12.176628 + 2.67250000 3.945663 + 2.67750000 20.068634 + 2.68250000 30.192291 + 2.68750000 10.316620 + 2.69250000 6.441628 + 2.69750000 10.567322 + 2.70250000 12.693687 + 2.70750000 9.820740 + 2.71250000 26.948471 + 2.71750000 29.076881 + 2.72250000 27.205971 + 2.72750000 21.335739 + 2.73250000 36.466187 + 2.73750000 58.597313 + 2.74250000 14.729126 + 2.74750000 35.861610 + 2.75250000 47.994781 + 2.75750000 40.128624 + 2.76250000 40.263153 + 2.76750000 56.398361 + 2.77250000 52.534248 + 2.77750000 39.670815 + 2.78250000 38.808060 + 2.78750000 38.945984 + 2.79250000 29.084595 + 2.79750000 66.223877 + 2.80250000 40.363838 + 2.80750000 19.504486 + 2.81250000 57.645813 + 2.81750000 35.787811 + 2.82250000 46.930496 + 2.82750000 61.073860 + 2.83250000 26.217903 + 2.83750000 20.362633 + 2.84250000 39.508034 + 2.84750000 44.654114 + 2.85250000 34.800880 + 2.85750000 42.948318 + 2.86250000 48.096443 + 2.86750000 39.245247 + 2.87250000 32.394722 + 2.87750000 31.544884 + 2.88250000 42.695724 + 2.88750000 55.847252 + 2.89250000 33.999451 + 2.89750000 58.152328 + 2.90250000 18.305885 + 2.90750000 28.460129 + 2.91250000 24.615051 + 2.91750000 29.770645 + 2.92250000 16.926926 + 2.92750000 26.083885 + 2.93250000 23.241524 + 2.93750000 14.399841 + 2.94250000 18.558838 + 2.94750000 32.718513 + 2.95250000 12.878876 + 2.95750000 23.039909 + 2.96250000 1.201630 + 2.96750000 4.364021 + 2.97250000 7.527100 + 2.97750000 4.690857 + 2.98250000 10.855293 + 2.98750000 1.020409 + 2.99250000 -2.813797 + 2.99750000 4.352676 + 3.00250000 12.519836 + 3.00750000 20.687668 + 3.01250000 -8.143822 + 3.01750000 -5.974625 + 3.02250000 -11.804749 + 3.02750000 -1.634201 + 3.03250000 8.537033 + 3.03750000 -6.291054 + 3.04250000 -21.118462 + 3.04750000 -5.945183 + 3.05250000 -6.771233 + 3.05750000 -5.596603 + 3.06250000 -11.421288 + 3.06750000 -9.245300 + 3.07250000 5.931374 + 3.07750000 -6.891273 + 3.08250000 -9.713249 + 3.08750000 -6.534538 + 3.09250000 -16.355148 + 3.09750000 -1.175079 + 3.10250000 1.005676 + 3.10750000 -8.812897 + 3.11250000 0.369209 + 3.11750000 -10.447998 + 3.12250000 -7.264534 + 3.12750000 5.919617 + 3.13250000 -2.895554 + 3.13750000 2.289955 + 3.14250000 0.476143 + 3.14750000 -10.336990 + 3.15250000 -2.149445 + 3.15750000 -6.961220 + 3.16250000 -1.772308 + 3.16750000 -14.582726 + 3.17250000 3.607544 + 3.17750000 14.798485 + 3.18250000 -1.009888 + 3.18750000 -9.817581 + 3.19250000 -18.624596 + 3.19750000 8.569069 + 3.20250000 -15.236588 + 3.20750000 -5.041565 + 3.21250000 5.154137 + 3.21750000 -11.649475 + 3.22250000 -6.452415 + 3.22750000 -7.254669 + 3.23250000 -12.056252 + 3.23750000 -15.857147 + 3.24250000 -20.657364 + 3.24750000 -3.456902 + 3.25250000 1.744240 + 3.25750000 -6.053940 + 3.26250000 -11.851440 + 3.26750000 10.351746 + 3.27250000 -2.444397 + 3.27750000 0.760147 + 3.28250000 -10.034637 + 3.28750000 -2.828735 + 3.29250000 -8.622154 + 3.29750000 6.585106 + 3.30250000 -13.206955 + 3.30750000 1.001663 + 3.31250000 -3.789040 + 3.31750000 -0.579056 + 3.32250000 -3.368401 + 3.32750000 -2.157066 + 3.33250000 0.054955 + 3.33750000 12.267654 + 3.34250000 5.481026 + 3.34750000 14.695084 + 3.35250000 3.909821 + 3.35750000 4.125237 + 3.36250000 0.341331 + 3.36750000 7.558113 + 3.37250000 -6.224434 + 3.37750000 -8.006302 + 3.38250000 2.212517 + 3.38750000 5.432014 + 3.39250000 -5.347816 + 3.39750000 12.873039 + 3.40250000 -1.905426 + 3.40750000 2.316788 + 3.41250000 -1.460320 + 3.41750000 -6.236748 + 3.42250000 -18.012489 + 3.42750000 5.212440 + 3.43250000 13.438057 + 3.43750000 0.664345 + 3.44250000 -5.108681 + 3.44750000 0.118973 + 3.45250000 9.347298 + 3.45750000 -3.423691 + 3.46250000 2.806000 + 3.46750000 10.036377 + 3.47250000 1.267426 + 3.47750000 11.499153 + 3.48250000 10.731567 + 3.48750000 11.964653 + 3.49250000 -0.801575 + 3.49750000 2.432877 diff --git a/ex_sub_50keV_2.txt b/ex_sub_50keV_2.txt new file mode 100644 index 0000000..f867c81 --- /dev/null +++ b/ex_sub_50keV_2.txt @@ -0,0 +1,800 @@ + -0.49750000 3.216999 + -0.49250000 11.968987 + -0.48750000 1.721531 + -0.48250000 6.474644 + -0.47750000 -7.771683 + -0.47250000 -1.017445 + -0.46750000 10.737358 + -0.46250000 12.492722 + -0.45750000 -7.751350 + -0.45250000 3.005142 + -0.44750000 0.762196 + -0.44250000 4.519814 + -0.43750000 -14.722004 + -0.43250000 -7.963257 + -0.42750000 0.796051 + -0.42250000 9.555923 + -0.41750000 -1.683643 + -0.41250000 -8.922646 + -0.40750000 -0.161083 + -0.40250000 7.601044 + -0.39750000 -11.636269 + -0.39250000 -1.873016 + -0.38750000 10.890800 + -0.38250000 1.655178 + -0.37750000 -0.579880 + -0.37250000 17.185627 + -0.36750000 -3.048306 + -0.36250000 0.718327 + -0.35750000 -11.514477 + -0.35250000 -6.746716 + -0.34750000 11.021606 + -0.34250000 -12.209507 + -0.33750000 -5.440060 + -0.33250000 0.329952 + -0.32750000 13.100529 + -0.32250000 -6.128330 + -0.31750000 11.643372 + -0.31250000 5.415638 + -0.30750000 -3.811531 + -0.30250000 2.961861 + -0.29750000 9.735817 + -0.29250000 7.510338 + -0.28750000 -1.714581 + -0.28250000 8.061066 + -0.27750000 -13.162724 + -0.27250000 -3.385948 + -0.26750000 -1.608612 + -0.26250000 -5.830711 + -0.25750000 -0.052250 + -0.25250000 -4.273224 + -0.24750000 2.506367 + -0.24250000 7.286522 + -0.23750000 26.067238 + -0.23250000 -6.151482 + -0.22750000 -2.369637 + -0.22250000 1.412769 + -0.21750000 9.195740 + -0.21250000 -7.020725 + -0.20750000 6.763371 + -0.20250000 21.548031 + -0.19750000 17.333256 + -0.19250000 16.119045 + -0.18750000 36.905396 + -0.18250000 27.692310 + -0.17750000 39.479786 + -0.17250000 57.267826 + -0.16750000 42.056431 + -0.16250000 55.845596 + -0.15750000 64.635330 + -0.15250000 72.425629 + -0.14750000 63.216484 + -0.14250000 104.007904 + -0.13750000 108.799889 + -0.13250000 127.592438 + -0.12750000 181.385559 + -0.12250000 165.179230 + -0.11750000 184.973465 + -0.11250000 233.768265 + -0.10750000 263.563629 + -0.10250000 288.359558 + -0.09750000 354.156067 + -0.09250000 403.953094 + -0.08750000 392.750732 + -0.08250000 463.548889 + -0.07750000 511.347656 + -0.07250000 490.146973 + -0.06750000 576.946838 + -0.06250000 605.747253 + -0.05750000 684.548279 + -0.05250000 680.349854 + -0.04750000 701.151978 + -0.04250000 819.954651 + -0.03750000 808.757874 + -0.03250000 818.561707 + -0.02750000 846.366089 + -0.02250000 847.171021 + -0.01750000 864.976562 + -0.01250000 947.782593 + -0.00750000 938.589233 + -0.00250000 972.396423 + 0.00250000 953.204224 + 0.00750000 906.012512 + 0.01250000 891.821411 + 0.01750000 933.630859 + 0.02250000 891.440857 + 0.02750000 870.251465 + 0.03250000 845.062561 + 0.03750000 818.874268 + 0.04250000 771.686523 + 0.04750000 736.499390 + 0.05250000 715.312744 + 0.05750000 745.126709 + 0.06250000 698.941284 + 0.06750000 674.756348 + 0.07250000 622.571960 + 0.07750000 584.388184 + 0.08250000 565.204956 + 0.08750000 530.022278 + 0.09250000 467.840210 + 0.09750000 478.658661 + 0.10250000 458.477692 + 0.10750000 426.297302 + 0.11250000 361.117432 + 0.11750000 342.938171 + 0.12250000 320.759430 + 0.12750000 343.581299 + 0.13250000 267.403687 + 0.13750000 266.226654 + 0.14250000 225.050201 + 0.14750000 217.874298 + 0.15250000 191.698944 + 0.15750000 163.524170 + 0.16250000 167.349976 + 0.16750000 135.176315 + 0.17250000 131.003235 + 0.17750000 116.830711 + 0.18250000 107.658752 + 0.18750000 97.487358 + 0.19250000 73.316521 + 0.19750000 69.146255 + 0.20250000 53.976547 + 0.20750000 44.807411 + 0.21250000 38.638832 + 0.21750000 64.470818 + 0.22250000 41.303360 + 0.22750000 51.136475 + 0.23250000 16.970146 + 0.23750000 30.804390 + 0.24250000 18.639191 + 0.24750000 31.474556 + 0.25250000 9.310486 + 0.25750000 3.146973 + 0.26250000 28.984032 + 0.26750000 24.821648 + 0.27250000 6.659836 + 0.27750000 42.498581 + 0.28250000 31.337891 + 0.28750000 14.177757 + 0.29250000 0.018196 + 0.29750000 23.859192 + 0.30250000 -4.299240 + 0.30750000 9.542885 + 0.31250000 9.385574 + 0.31750000 -6.771172 + 0.32250000 17.072639 + 0.32750000 -2.082977 + 0.33250000 7.761963 + 0.33750000 6.607468 + 0.34250000 8.453537 + 0.34750000 14.300171 + 0.35250000 1.147369 + 0.35750000 1.995125 + 0.36250000 8.843452 + 0.36750000 14.692337 + 0.37250000 9.541786 + 0.37750000 5.391800 + 0.38250000 4.242378 + 0.38750000 -13.906479 + 0.39250000 -5.054779 + 0.39750000 0.797485 + 0.40250000 7.650314 + 0.40750000 -6.496292 + 0.41250000 -6.642334 + 0.41750000 9.212189 + 0.42250000 -26.932724 + 0.42750000 2.922920 + 0.43250000 6.779129 + 0.43750000 -1.364098 + 0.44250000 6.493240 + 0.44750000 -1.648857 + 0.45250000 11.209610 + 0.45750000 -13.931366 + 0.46250000 18.928223 + 0.46750000 -2.211624 + 0.47250000 -6.350906 + 0.47750000 7.510376 + 0.48250000 -13.627777 + 0.48750000 -8.765373 + 0.49250000 10.097603 + 0.49750000 -17.038864 + 0.50250000 6.825233 + 0.50750000 5.689896 + 0.51250000 3.555122 + 0.51750000 -8.579094 + 0.52250000 -6.712738 + 0.52750000 -0.845825 + 0.53250000 10.021652 + 0.53750000 5.889694 + 0.54250000 -1.241699 + 0.54750000 13.627464 + 0.55250000 10.497200 + 0.55750000 -7.632507 + 0.56250000 5.238350 + 0.56750000 23.109772 + 0.57250000 -14.018242 + 0.57750000 3.854309 + 0.58250000 -12.272583 + 0.58750000 -12.398903 + 0.59250000 -0.524666 + 0.59750000 2.350136 + 0.60250000 6.225502 + 0.60750000 9.101433 + 0.61250000 24.977921 + 0.61750000 1.854980 + 0.62250000 7.732597 + 0.62750000 2.610779 + 0.63250000 17.489525 + 0.63750000 11.368835 + 0.64250000 4.248711 + 0.64750000 10.129143 + 0.65250000 25.010147 + 0.65750000 23.891708 + 0.66250000 28.773834 + 0.66750000 42.656525 + 0.67250000 44.539772 + 0.67750000 42.423592 + 0.68250000 41.307968 + 0.68750000 70.192917 + 0.69250000 99.078423 + 0.69750000 82.964493 + 0.70250000 110.851120 + 0.70750000 134.738312 + 0.71250000 168.626083 + 0.71750000 204.514404 + 0.72250000 197.403290 + 0.72750000 254.292740 + 0.73250000 298.182739 + 0.73750000 284.073334 + 0.74250000 336.964478 + 0.74750000 368.856171 + 0.75250000 427.748444 + 0.75750000 419.641266 + 0.76250000 483.534668 + 0.76750000 505.428619 + 0.77250000 532.323120 + 0.77750000 564.218262 + 0.78250000 606.113892 + 0.78750000 653.010071 + 0.79250000 684.906860 + 0.79750000 669.804199 + 0.80250000 724.702087 + 0.80750000 727.600586 + 0.81250000 688.499573 + 0.81750000 748.399170 + 0.82250000 731.299316 + 0.82750000 755.200073 + 0.83250000 740.101318 + 0.83750000 817.003174 + 0.84250000 706.905579 + 0.84750000 789.808594 + 0.85250000 779.712097 + 0.85750000 741.616211 + 0.86250000 789.520874 + 0.86750000 761.426086 + 0.87250000 720.331848 + 0.87750000 706.238220 + 0.88250000 690.145142 + 0.88750000 698.052612 + 0.89250000 667.960632 + 0.89750000 615.869263 + 0.90250000 587.778442 + 0.90750000 621.688171 + 0.91250000 526.598450 + 0.91750000 572.509338 + 0.92250000 538.420715 + 0.92750000 486.332733 + 0.93250000 508.245270 + 0.93750000 421.158386 + 0.94250000 433.072052 + 0.94750000 395.986298 + 0.95250000 356.901093 + 0.95750000 360.816467 + 0.96250000 348.732391 + 0.96750000 281.648895 + 0.97250000 242.565948 + 0.97750000 254.483566 + 0.98250000 250.401749 + 0.98750000 215.320496 + 0.99250000 235.239807 + 0.99750000 194.159668 + 1.00250000 185.080109 + 1.00750000 162.001114 + 1.01250000 128.922668 + 1.01750000 134.844803 + 1.02250000 123.767494 + 1.02750000 118.690742 + 1.03250000 88.614563 + 1.03750000 90.538940 + 1.04250000 67.463882 + 1.04750000 73.389389 + 1.05250000 57.315460 + 1.05750000 65.242096 + 1.06250000 49.169289 + 1.06750000 50.097054 + 1.07250000 30.025375 + 1.07750000 31.954262 + 1.08250000 38.883713 + 1.08750000 35.813728 + 1.09250000 30.744308 + 1.09750000 30.675446 + 1.10250000 17.607147 + 1.10750000 4.539413 + 1.11250000 -1.527756 + 1.11750000 22.405640 + 1.12250000 13.339600 + 1.12750000 16.274117 + 1.13250000 31.209206 + 1.13750000 1.144852 + 1.14250000 3.081062 + 1.14750000 9.017838 + 1.15250000 6.955177 + 1.15750000 5.893074 + 1.16250000 28.831543 + 1.16750000 7.770569 + 1.17250000 12.710159 + 1.17750000 24.650314 + 1.18250000 14.591034 + 1.18750000 35.532318 + 1.19250000 30.474159 + 1.19750000 36.416565 + 1.20250000 51.359543 + 1.20750000 39.303078 + 1.21250000 52.247177 + 1.21750000 86.191833 + 1.22250000 97.137062 + 1.22750000 89.082848 + 1.23250000 101.029198 + 1.23750000 120.976120 + 1.24250000 103.923599 + 1.24750000 152.871643 + 1.25250000 184.820251 + 1.25750000 166.769409 + 1.26250000 166.719147 + 1.26750000 204.669434 + 1.27250000 257.620300 + 1.27750000 251.571716 + 1.28250000 261.523712 + 1.28750000 289.476257 + 1.29250000 248.429367 + 1.29750000 312.383057 + 1.30250000 310.337280 + 1.30750000 315.292084 + 1.31250000 329.247437 + 1.31750000 314.203369 + 1.32250000 305.159882 + 1.32750000 319.116943 + 1.33250000 360.074554 + 1.33750000 336.032745 + 1.34250000 300.991486 + 1.34750000 314.950806 + 1.35250000 324.910675 + 1.35750000 355.871094 + 1.36250000 327.832092 + 1.36750000 353.793671 + 1.37250000 333.755798 + 1.37750000 334.718475 + 1.38250000 318.681732 + 1.38750000 307.645569 + 1.39250000 288.609924 + 1.39750000 300.574890 + 1.40250000 281.540375 + 1.40750000 293.506470 + 1.41250000 280.473083 + 1.41750000 289.440277 + 1.42250000 265.408051 + 1.42750000 256.376373 + 1.43250000 248.345245 + 1.43750000 255.314697 + 1.44250000 256.284729 + 1.44750000 254.255295 + 1.45250000 241.226440 + 1.45750000 239.198151 + 1.46250000 252.170410 + 1.46750000 212.143250 + 1.47250000 241.116638 + 1.47750000 209.090591 + 1.48250000 258.065125 + 1.48750000 169.040192 + 1.49250000 206.015839 + 1.49750000 219.992065 + 1.50250000 181.968842 + 1.50750000 199.946167 + 1.51250000 178.924072 + 1.51750000 168.902542 + 1.52250000 137.881561 + 1.52750000 155.861160 + 1.53250000 148.841309 + 1.53750000 143.822021 + 1.54250000 151.803314 + 1.54750000 130.785156 + 1.55250000 144.767578 + 1.55750000 123.750542 + 1.56250000 111.734077 + 1.56750000 127.718178 + 1.57250000 108.702835 + 1.57750000 121.688065 + 1.58250000 97.673851 + 1.58750000 89.660202 + 1.59250000 92.647118 + 1.59750000 94.634598 + 1.60250000 73.622643 + 1.60750000 76.611244 + 1.61250000 63.600418 + 1.61750000 57.590149 + 1.62250000 52.580444 + 1.62750000 86.571304 + 1.63250000 64.562729 + 1.63750000 65.554710 + 1.64250000 33.547264 + 1.64750000 57.540375 + 1.65250000 58.534050 + 1.65750000 68.528290 + 1.66250000 53.523094 + 1.66750000 39.518463 + 1.67250000 33.514389 + 1.67750000 15.510880 + 1.68250000 22.507942 + 1.68750000 25.505562 + 1.69250000 10.503746 + 1.69750000 28.502487 + 1.70250000 18.501801 + 1.70750000 28.501671 + 1.71250000 9.502106 + 1.71750000 7.503113 + 1.72250000 6.504677 + 1.72750000 17.506798 + 1.73250000 16.509491 + 1.73750000 6.512741 + 1.74250000 22.516563 + 1.74750000 16.520943 + 1.75250000 7.525887 + 1.75750000 12.531395 + 1.76250000 7.537468 + 1.76750000 3.544098 + 1.77250000 7.551292 + 1.77750000 35.559059 + 1.78250000 14.567383 + 1.78750000 -2.423729 + 1.79250000 17.585716 + 1.79750000 39.595734 + 1.80250000 48.606316 + 1.80750000 32.617455 + 1.81250000 37.629158 + 1.81750000 58.641426 + 1.82250000 65.654259 + 1.82750000 77.667648 + 1.83250000 115.681610 + 1.83750000 136.696136 + 1.84250000 142.711212 + 1.84750000 183.726868 + 1.85250000 177.743073 + 1.85750000 189.759857 + 1.86250000 251.777191 + 1.86750000 297.795105 + 1.87250000 365.813538 + 1.87750000 376.832581 + 1.88250000 411.852173 + 1.88750000 486.872314 + 1.89250000 524.893066 + 1.89750000 557.914307 + 1.90250000 619.936157 + 1.90750000 697.958618 + 1.91250000 708.981567 + 1.91750000 774.005127 + 1.92250000 745.029175 + 1.92750000 784.053833 + 1.93250000 822.079102 + 1.93750000 817.104858 + 1.94250000 779.131226 + 1.94750000 803.158142 + 1.95250000 811.185608 + 1.95750000 886.213684 + 1.96250000 906.242249 + 1.96750000 874.271423 + 1.97250000 894.301147 + 1.97750000 903.331482 + 1.98250000 887.362305 + 1.98750000 945.393738 + 1.99250000 880.425720 + 1.99750000 955.458252 + 2.00250000 904.491394 + 2.00750000 881.525085 + 2.01250000 884.559326 + 2.01750000 821.594116 + 2.02250000 932.629517 + 2.02750000 796.665405 + 2.03250000 870.701904 + 2.03750000 758.738953 + 2.04250000 789.776611 + 2.04750000 741.814758 + 2.05250000 734.853516 + 2.05750000 753.892822 + 2.06250000 713.932678 + 2.06750000 605.973145 + 2.07250000 636.014160 + 2.07750000 604.055725 + 2.08250000 575.097839 + 2.08750000 571.140503 + 2.09250000 554.183777 + 2.09750000 502.227600 + 2.10250000 459.271973 + 2.10750000 429.316925 + 2.11250000 437.362427 + 2.11750000 436.408508 + 2.12250000 396.455139 + 2.12750000 387.502319 + 2.13250000 335.550079 + 2.13750000 301.598419 + 2.14250000 301.647308 + 2.14750000 282.696747 + 2.15250000 251.746765 + 2.15750000 245.797348 + 2.16250000 230.848480 + 2.16750000 220.900192 + 2.17250000 191.952454 + 2.17750000 169.005295 + 2.18250000 165.058685 + 2.18750000 159.112640 + 2.19250000 170.167175 + 2.19750000 143.222260 + 2.20250000 138.277893 + 2.20750000 104.334114 + 2.21250000 73.390884 + 2.21750000 85.448227 + 2.22250000 69.506126 + 2.22750000 70.564590 + 2.23250000 72.623619 + 2.23750000 49.683212 + 2.24250000 55.743370 + 2.24750000 61.804092 + 2.25250000 55.865372 + 2.25750000 41.927216 + 2.26250000 60.989624 + 2.26750000 26.052597 + 2.27250000 40.116135 + 2.27750000 29.180237 + 2.28250000 8.244896 + 2.28750000 45.310120 + 2.29250000 15.375916 + 2.29750000 13.442268 + 2.30250000 19.509178 + 2.30750000 4.576660 + 2.31250000 -4.355293 + 2.31750000 16.713310 + 2.32250000 23.782478 + 2.32750000 4.852211 + 2.33250000 0.922508 + 2.33750000 17.993370 + 2.34250000 10.064796 + 2.34750000 4.136780 + 2.35250000 -1.790672 + 2.35750000 3.282440 + 2.36250000 -17.643883 + 2.36750000 7.430359 + 2.37250000 14.505165 + 2.37750000 -29.419472 + 2.38250000 5.656464 + 2.38750000 -2.267044 + 2.39250000 13.810013 + 2.39750000 -1.112366 + 2.40250000 -4.034180 + 2.40750000 3.044563 + 2.41250000 1.123871 + 2.41750000 -9.796249 + 2.42250000 -0.715813 + 2.42750000 11.365189 + 2.43250000 -5.553246 + 2.43750000 4.528877 + 2.44250000 -8.388428 + 2.44750000 -8.305176 + 2.45250000 -16.221359 + 2.45750000 -9.136978 + 2.46250000 -8.052032 + 2.46750000 4.033478 + 2.47250000 -11.880455 + 2.47750000 -3.793816 + 2.48250000 -9.706619 + 2.48750000 -19.618858 + 2.49250000 -21.530533 + 2.49750000 -6.441643 + 2.50250000 -11.352196 + 2.50750000 -24.262177 + 2.51250000 -10.171600 + 2.51750000 -11.080460 + 2.52250000 -7.988754 + 2.52750000 -8.896484 + 2.53250000 -5.803658 + 2.53750000 -1.710258 + 2.54250000 2.383698 + 2.54750000 -0.521774 + 2.55250000 -11.426689 + 2.55750000 -25.331047 + 2.56250000 -1.234833 + 2.56750000 -15.138054 + 2.57250000 5.959282 + 2.57750000 1.057182 + 2.58250000 -9.844353 + 2.58750000 -8.745323 + 2.59250000 12.354271 + 2.59750000 -2.545570 + 2.60250000 -15.444855 + 2.60750000 -1.343567 + 2.61250000 -7.241722 + 2.61750000 -17.139313 + 2.62250000 -8.036339 + 2.62750000 -8.932808 + 2.63250000 5.171295 + 2.63750000 -13.724045 + 2.64250000 0.381187 + 2.64750000 4.486977 + 2.65250000 -8.406677 + 2.65750000 18.700241 + 2.66250000 15.807724 + 2.66750000 -6.084236 + 2.67250000 10.024376 + 2.67750000 26.133545 + 2.68250000 36.243279 + 2.68750000 16.353577 + 2.69250000 12.464432 + 2.69750000 16.575859 + 2.70250000 18.687843 + 2.70750000 15.800392 + 2.71250000 32.913506 + 2.71750000 35.027184 + 2.72250000 33.141426 + 2.72750000 27.256226 + 2.73250000 42.371597 + 2.73750000 64.487526 + 2.74250000 20.604019 + 2.74750000 41.721077 + 2.75250000 53.838699 + 2.75750000 45.956886 + 2.76250000 46.075630 + 2.76750000 62.194939 + 2.77250000 58.314819 + 2.77750000 45.435257 + 2.78250000 44.556252 + 2.78750000 44.677818 + 2.79250000 34.799950 + 2.79750000 71.922638 + 2.80250000 46.045891 + 2.80750000 25.169708 + 2.81250000 63.294090 + 2.81750000 41.419037 + 2.82250000 52.544548 + 2.82750000 66.670616 + 2.83250000 31.797249 + 2.83750000 25.924454 + 2.84250000 45.052216 + 2.84750000 50.180534 + 2.85250000 40.309425 + 2.85750000 48.438881 + 2.86250000 53.568893 + 2.86750000 44.699471 + 2.87250000 37.830612 + 2.87750000 36.962318 + 2.88250000 48.094589 + 2.88750000 61.227417 + 2.89250000 39.360817 + 2.89750000 63.494774 + 2.90250000 23.629295 + 2.90750000 33.764381 + 2.91250000 29.900032 + 2.91750000 35.036247 + 2.92250000 22.173019 + 2.92750000 31.310356 + 2.93250000 28.448265 + 2.93750000 19.586731 + 2.94250000 23.725754 + 2.94750000 37.865349 + 2.95250000 18.005508 + 2.95750000 28.146225 + 2.96250000 6.287506 + 2.96750000 9.429352 + 2.97250000 12.571762 + 2.97750000 9.714737 + 2.98250000 15.858276 + 2.98750000 6.002373 + 2.99250000 2.147041 + 2.99750000 9.292267 + 3.00250000 17.438057 + 3.00750000 25.584404 + 3.01250000 -3.268677 + 3.01750000 -1.121193 + 3.02250000 -6.973152 + 3.02750000 3.175453 + 3.03250000 13.324623 + 3.03750000 -1.525642 + 3.04250000 -16.375343 + 3.04750000 -1.224480 + 3.05250000 -2.073059 + 3.05750000 -0.921074 + 3.06250000 -6.768524 + 3.06750000 -4.615410 + 3.07250000 10.538269 + 3.07750000 -2.307487 + 3.08250000 -5.152687 + 3.08750000 -1.997314 + 3.09250000 -11.841385 + 3.09750000 3.315109 + 3.10250000 5.472168 + 3.10750000 -4.370216 + 3.11250000 4.787971 + 3.11750000 -6.053284 + 3.12250000 -2.893967 + 3.12750000 10.265907 + 3.13250000 1.426346 + 3.13750000 6.587341 + 3.14250000 4.748909 + 3.14750000 -6.088966 + 3.15250000 2.073730 + 3.15750000 -2.763016 + 3.16250000 2.400803 + 3.16750000 -10.434814 + 3.17250000 7.730125 + 3.17750000 18.895638 + 3.18250000 3.061707 + 3.18750000 -5.771660 + 3.19250000 -14.604454 + 3.19750000 12.563301 + 3.20250000 -11.268372 + 3.20750000 -1.099480 + 3.21250000 9.069969 + 3.21750000 -7.760017 + 3.22250000 -2.589432 + 3.22750000 -3.418289 + 3.23250000 -8.246590 + 3.23750000 -12.074318 + 3.24250000 -16.901489 + 3.24750000 0.271912 + 3.25250000 5.445869 + 3.25750000 -2.379608 + 3.26250000 -8.204521 + 3.26750000 13.971130 + 3.27250000 1.147339 + 3.27750000 4.324120 + 3.28250000 -6.498543 + 3.28750000 0.679359 + 3.29250000 -5.142174 + 3.29750000 10.036858 + 3.30250000 -9.783554 + 3.30750000 4.396606 + 3.31250000 -0.422676 + 3.31750000 2.758606 + 3.32250000 -0.059547 + 3.32750000 1.122864 + 3.33250000 3.305840 + 3.33750000 15.489372 + 3.34250000 8.673477 + 3.34750000 17.858139 + 3.35250000 7.043365 + 3.35750000 7.229156 + 3.36250000 3.415512 + 3.36750000 10.602425 + 3.37250000 -3.210091 + 3.37750000 -5.022049 + 3.38250000 5.166557 + 3.38750000 8.355728 + 3.39250000 -2.454536 + 3.39750000 15.735756 + 3.40250000 0.926620 + 3.40750000 5.118042 + 3.41250000 1.310028 + 3.41750000 -3.497421 + 3.42250000 -15.304306 + 3.42750000 7.889374 + 3.43250000 16.083618 + 3.43750000 3.278419 + 3.44250000 -2.526215 + 3.44750000 2.669716 + 3.45250000 11.866211 + 3.45750000 -0.936729 + 3.46250000 5.260895 + 3.46750000 12.459076 + 3.47250000 3.657829 + 3.47750000 13.857140 + 3.48250000 13.057014 + 3.48750000 14.257454 + 3.49250000 1.458450 + 3.49750000 4.660019 diff --git a/fit_result.txt b/fit_result.txt new file mode 100644 index 0000000..1b9233a --- /dev/null +++ b/fit_result.txt @@ -0,0 +1,1052 @@ +0 1 30 -0.040 +0 1 31 0.006 +0 1 32 -0.025 +0 1 33 -0.054 +0 1 34 0.021 +0 1 35 0.067 +0 1 36 0.043 +0 1 37 0.022 +0 1 38 -0.041 +0 1 39 0.051 +0 1 40 0.039 +0 1 41 0.058 +0 1 42 0.031 +0 1 43 0.048 +0 1 44 0.040 +0 1 45 0.037 +0 1 46 0.068 +0 1 47 0.043 +0 1 48 0.046 +0 1 49 0.024 +0 1 50 0.068 +0 1 51 0.026 +0 1 52 0.076 +0 1 53 0.024 +0 1 54 0.039 +0 1 55 0.041 +0 1 56 0.037 +0 1 57 0.084 +0 1 58 0.049 +0 1 59 0.063 +0 1 60 0.050 +0 1 61 0.092 +0 1 62 0.030 +0 1 63 0.050 +0 1 64 0.062 +0 1 65 0.027 +0 1 66 0.025 +0 1 67 0.075 +0 1 68 0.035 +0 1 69 0.037 +0 1 70 0.055 +0 1 71 0.057 +0 1 72 0.051 +0 1 73 0.058 +0 1 74 0.060 +0 1 75 0.042 +0 1 76 0.054 +0 1 77 -0.007 +0 1 78 0.030 +0 1 79 0.051 +0 1 80 0.067 +0 1 81 0.063 +0 1 82 0.072 +0 1 83 0.046 +0 1 84 0.052 +0 1 85 0.070 +0 1 86 0.034 0.849 +0 1 87 0.038 +0 1 88 0.042 0.811 +0 1 89 0.051 0.819 +0 1 90 0.041 0.835 +0 1 91 0.049 0.856 +0 1 92 0.028 0.856 +0 1 93 0.071 0.860 +0 1 94 0.022 0.882 +0 1 95 0.056 0.874 +0 1 96 0.044 0.830 +0 1 97 0.023 0.845 +0 1 98 0.056 0.841 +0 1 99 0.048 0.864 +0 1 100 0.038 0.877 +0 1 101 0.079 0.871 +0 1 102 0.019 0.897 +0 1 103 0.069 0.911 +0 1 104 0.030 0.877 +0 1 105 0.047 0.878 +0 1 106 0.093 0.887 +0 1 107 0.024 0.897 +0 1 108 0.025 0.891 +0 1 109 0.076 0.881 +0 1 110 0.071 0.941 +0 1 111 0.040 0.914 +0 1 112 0.036 0.893 +0 1 113 0.049 0.872 +0 1 114 0.080 0.866 +0 1 115 0.082 0.898 +0 1 116 0.063 0.923 +0 1 117 0.074 0.899 +0 1 118 0.011 0.896 +0 1 119 0.077 0.888 +0 1 120 0.057 0.900 +0 1 121 0.043 0.929 +0 1 122 0.020 0.884 +0 1 123 0.050 0.934 +0 1 124 0.064 0.929 +0 1 125 0.057 0.914 +0 1 126 0.023 0.911 +0 1 127 0.075 0.934 +0 2 0 -0.013 0.833 1.959 +0 2 1 -0.054 0.800 1.906 +0 2 2 -0.019 0.791 1.870 +0 2 3 -0.030 0.759 1.902 +0 2 4 -0.056 0.827 1.906 +0 2 5 0.011 0.820 1.988 +0 2 6 -0.023 0.759 1.946 +0 2 7 -0.077 0.747 1.955 +0 2 8 -0.039 0.809 1.988 +0 2 9 -0.008 0.831 1.975 +0 2 10 -0.021 0.846 1.938 +0 2 11 -0.042 0.772 1.957 +0 2 12 -0.041 0.829 1.976 +0 2 13 -0.080 0.829 1.955 +0 2 14 -0.041 0.838 1.987 +0 2 15 0.007 0.851 1.990 +0 2 16 -0.045 0.816 2.002 +0 2 17 -0.050 0.829 1.974 +0 2 18 -0.015 0.791 1.960 +0 2 19 0.004 0.788 1.961 +0 2 20 0.012 0.804 1.998 +0 2 21 -0.102 0.826 1.995 +0 2 22 -0.042 0.794 1.984 +0 2 23 -0.044 0.820 1.964 +0 2 24 0.011 0.811 1.998 +0 2 25 -0.039 0.843 1.969 +0 2 26 -0.050 0.835 1.972 +0 2 27 -0.026 0.816 1.978 +0 2 28 -0.038 0.793 1.972 +0 2 29 -0.046 0.790 2.002 +0 2 30 -0.044 0.814 2.006 +0 2 31 -0.034 0.759 2.042 +0 2 32 -0.043 0.466 1.964 +0 2 33 -0.031 0.789 1.943 +0 2 34 -0.055 0.786 1.991 +0 2 35 -0.033 0.824 2.015 +0 2 36 -0.030 0.805 1.984 +0 2 37 -0.025 0.845 1.968 +0 2 38 -0.009 0.810 1.979 +0 2 39 -0.018 0.774 1.948 +0 2 40 -0.078 0.788 2.001 +0 2 41 -0.028 0.832 2.003 +0 2 42 -0.012 0.819 1.959 +0 2 43 -0.027 0.786 1.980 +0 2 44 -0.003 0.833 1.957 +0 2 45 0.001 0.834 2.005 +0 2 46 -0.028 0.832 1.959 +0 2 47 -0.035 0.828 1.971 +0 2 48 -0.032 0.806 1.959 +0 2 49 -0.019 0.831 1.972 +0 2 50 -0.030 0.813 1.988 +0 2 51 -0.005 0.741 1.998 +0 2 52 -0.148 0.651 1.828 +0 2 53 -0.030 0.777 1.982 +0 2 54 -0.125 0.700 1.847 +0 2 55 -0.041 0.835 1.957 +0 2 56 -0.012 0.804 1.987 +0 2 57 0.040 0.849 1.961 +0 2 58 -0.038 0.815 1.960 +0 2 59 -0.084 0.801 2.010 +0 2 60 -0.031 0.818 1.962 +0 2 61 0.011 0.799 1.966 +0 2 62 -0.003 0.839 1.970 +0 2 63 -0.026 0.762 1.962 +0 2 64 0.093 0.838 1.979 +0 2 65 -0.028 0.844 1.964 +0 2 66 -0.018 0.812 1.984 +0 2 67 -0.021 0.786 1.982 +0 2 68 -0.029 0.824 1.964 +0 2 69 -0.024 0.791 1.978 +0 2 70 0.002 0.840 1.969 +0 2 71 -0.027 0.790 1.994 +0 2 72 -0.025 0.743 1.947 +0 2 73 0.018 0.806 1.990 +0 2 74 -0.044 0.793 1.953 +0 2 75 -0.012 0.820 1.941 +0 2 76 -0.022 0.779 1.965 +0 2 77 -0.019 0.827 1.963 +0 2 78 -0.014 0.831 1.982 +0 2 79 -0.014 0.800 1.982 +0 2 80 -0.023 0.816 1.989 +0 2 81 -0.003 0.767 1.941 +0 2 82 -0.026 0.815 1.981 +0 2 83 -0.047 0.835 1.963 +0 2 84 -0.063 0.786 1.941 +0 2 85 -0.047 0.779 1.966 +0 2 86 -0.021 0.786 1.966 +0 2 87 -0.043 0.824 1.962 +0 2 88 -0.004 0.811 2.005 +0 2 89 0.020 0.796 1.964 +0 2 90 -0.014 0.799 2.003 +0 2 91 -0.055 0.787 1.953 +0 2 92 -0.022 0.861 1.989 +0 2 93 0.001 0.816 1.997 +0 2 94 -0.035 0.843 1.983 +0 2 95 0.001 0.800 1.957 +0 2 96 0.007 0.788 1.973 +0 2 97 -0.051 0.814 1.969 +0 2 98 -0.023 0.812 1.977 +0 2 99 0.039 0.765 1.979 +0 2 100 0.011 0.814 1.992 +0 2 101 -0.027 0.861 1.944 +0 2 102 -0.133 0.841 1.977 +0 2 103 -0.000 0.832 1.939 +0 2 104 -0.019 0.839 1.973 +0 2 105 -0.013 0.819 2.020 +0 2 106 -0.055 0.836 1.939 +0 2 107 -0.006 0.822 1.965 +0 2 108 -0.105 0.804 1.953 +0 2 109 -0.029 0.839 1.979 +0 2 110 0.013 0.809 1.977 +0 2 111 -0.047 0.810 1.954 +0 2 112 -0.002 0.791 1.978 +0 2 113 0.037 0.853 1.964 +0 2 114 -0.009 0.801 2.043 +0 2 115 -0.039 0.825 1.976 +0 2 116 -0.007 0.815 1.934 +0 2 117 -0.032 0.857 1.953 +0 2 118 -0.030 0.787 1.965 +0 2 119 -0.046 0.825 1.966 +0 2 120 -0.065 0.807 1.969 +0 2 121 0.003 0.804 1.985 +0 2 122 0.009 0.831 1.963 +0 2 123 -0.010 0.859 1.980 +0 2 124 -0.041 0.732 1.970 +0 2 125 -0.030 0.827 1.938 +0 2 126 -0.039 0.797 1.917 +0 2 127 -0.046 0.854 1.948 +0 3 0 0.013 0.844 2.028 +0 3 1 -0.015 0.880 1.990 +0 3 2 0.003 0.833 2.000 +0 3 3 -0.021 0.839 2.026 +0 3 4 -0.017 0.861 2.000 +0 3 5 -0.019 0.873 2.031 +0 3 6 0.010 0.863 1.965 +0 3 7 0.008 0.869 1.986 +0 3 8 0.006 0.894 2.019 +0 3 9 0.005 0.875 2.032 +0 3 10 0.003 0.874 1.973 +0 3 11 0.006 0.855 2.027 +0 3 12 0.007 0.835 1.985 +0 3 13 0.011 0.876 1.979 +0 3 14 0.006 0.862 2.009 +0 3 15 -0.006 0.852 2.016 +0 3 16 0.013 0.870 2.006 +0 3 17 -0.011 0.844 2.026 +0 3 18 -0.009 0.854 1.984 +0 3 19 -0.130 0.857 2.034 +0 3 20 0.018 0.863 1.983 +0 3 21 -0.035 0.888 2.017 +0 3 22 -0.056 0.862 2.020 +0 3 23 -0.004 0.887 1.983 +0 3 24 0.011 0.856 2.008 +0 3 25 0.024 0.875 1.998 +0 3 26 0.001 0.810 2.023 +0 3 27 -0.076 0.861 2.004 +0 3 28 0.010 0.854 1.982 +0 3 29 -0.011 0.803 2.003 +0 3 30 -0.016 0.873 2.000 +0 3 31 0.003 0.852 2.025 +0 3 32 -0.000 0.877 1.825 +0 3 33 -0.021 0.746 1.986 +0 3 34 0.009 0.862 1.974 +0 3 35 0.009 0.849 1.975 +0 3 36 -0.081 0.879 1.977 +0 3 37 -0.034 0.836 1.981 +0 3 38 -0.019 0.847 2.018 +0 3 39 0.039 0.852 2.011 +0 3 40 -0.037 0.886 2.000 +0 3 41 -0.061 0.828 2.003 +0 3 42 -0.052 0.819 2.001 +0 3 43 0.011 0.869 2.035 +0 3 44 0.017 0.810 2.027 +0 3 45 0.010 0.857 2.012 +0 3 46 -0.009 0.400 2.008 +0 3 47 -0.029 0.826 1.996 +0 3 48 -0.023 0.842 2.025 +0 3 49 0.003 0.858 1.992 +0 3 50 0.037 0.854 2.004 +0 3 51 -0.038 0.845 2.004 +0 3 52 -0.123 0.887 1.997 +0 3 53 -0.029 0.842 1.986 +0 3 54 -0.017 0.873 1.990 +0 3 55 -0.002 0.843 2.038 +0 3 56 0.013 0.748 2.021 +0 3 57 -0.174 0.837 2.010 +0 3 58 0.015 0.841 2.003 +0 3 59 -0.023 0.870 1.987 +0 3 60 0.007 0.867 2.029 +0 3 61 -0.021 0.831 2.004 +0 3 62 -0.028 0.837 2.033 +0 3 63 0.051 0.860 1.950 +0 3 64 -0.083 0.834 1.926 +0 3 65 0.016 0.845 2.002 +0 3 66 -0.030 0.841 1.999 +0 3 67 0.017 0.800 2.005 +0 3 68 0.035 0.837 1.981 +0 3 69 -0.016 0.820 1.989 +0 3 70 -0.029 0.848 2.005 +0 3 71 -0.049 0.846 2.009 +0 3 72 0.019 0.876 2.005 +0 3 73 0.002 0.842 2.018 +0 3 74 -0.005 0.854 1.975 +0 3 75 -0.001 0.824 2.022 +0 3 76 -0.008 0.848 2.011 +0 3 77 -0.053 0.826 1.946 +0 3 78 0.004 0.793 1.978 +0 3 79 0.004 0.829 1.966 +0 3 80 -0.039 0.825 1.980 +0 3 81 -0.081 0.839 2.010 +0 3 82 0.003 0.843 1.990 +0 3 83 -0.020 0.806 1.961 +0 3 84 -0.014 0.813 2.001 +0 3 85 -0.034 0.861 2.014 +0 3 86 -0.027 0.815 1.941 +0 3 87 0.022 0.826 1.976 +0 3 88 -0.075 0.789 1.960 +0 3 89 -0.075 0.781 1.980 +0 3 90 -0.046 0.835 1.954 +0 3 91 -0.015 0.796 1.988 +0 3 92 -0.026 0.842 1.993 +0 3 93 -0.076 0.811 1.972 +0 3 94 0.040 0.866 2.000 +0 3 95 -0.066 0.844 1.958 +0 3 96 -0.082 0.810 2.006 +0 3 97 -0.032 0.837 2.011 +0 3 98 0.010 0.795 1.966 +0 3 99 -0.037 0.812 1.990 +0 3 100 -0.060 0.790 1.977 +0 3 101 -0.083 0.826 1.986 +0 3 102 -0.038 0.822 1.965 +0 3 103 -0.008 0.760 1.947 +0 3 104 -0.062 0.817 1.982 +0 3 105 -0.075 0.775 1.925 +0 3 106 -0.005 0.797 1.927 +0 3 107 0.075 0.825 1.975 +0 3 108 -0.019 0.795 1.910 +0 3 109 0.000 0.871 1.967 +0 3 110 0.014 0.741 1.984 +0 3 111 0.776 1.959 +0 3 112 0.794 1.966 +0 3 113 0.769 1.980 +0 3 114 0.775 1.925 +0 3 115 -0.075 1.932 +0 3 116 0.785 1.943 +0 3 117 0.778 1.961 +0 3 118 0.825 1.925 +0 3 119 0.756 1.948 +0 3 120 0.749 1.975 +0 3 121 -0.075 1.958 +0 3 123 -0.076 0.973 1.977 +0 3 124 0.025 0.725 1.925 +0 3 125 -0.125 0.825 1.925 +0 3 126 0.764 1.952 +1 1 30 0.030 +1 1 31 0.045 +1 1 32 0.068 +1 1 33 0.051 +1 1 34 0.066 +1 1 35 0.054 +1 1 36 0.054 +1 1 37 0.066 +1 1 38 0.091 +1 1 39 0.077 +1 1 40 0.082 +1 1 41 -0.031 +1 1 42 0.043 +1 1 43 0.083 +1 1 44 0.102 +1 1 45 0.051 +1 1 46 0.047 +1 1 47 0.054 +1 1 48 0.049 +1 1 49 0.054 +1 1 50 0.046 +1 1 51 0.060 +1 1 52 0.050 +1 1 53 0.035 +1 1 54 0.036 +1 1 55 0.048 +1 1 56 0.090 +1 1 57 0.050 +1 1 58 0.033 +1 1 59 0.079 +1 1 60 0.068 +1 1 61 0.048 +1 1 62 0.023 +1 1 63 0.043 +1 1 64 0.035 +1 1 65 0.051 +1 1 66 0.047 +1 1 67 0.034 +1 1 68 0.030 +1 1 69 0.075 0.875 +1 1 70 0.075 0.875 +1 1 71 0.044 +1 1 72 0.049 0.874 +1 1 73 0.038 0.912 +1 1 74 0.019 0.915 +1 1 75 0.025 0.925 +1 1 76 0.035 0.937 +1 1 77 0.056 0.926 +1 1 78 0.025 0.975 +1 1 79 0.059 0.942 +1 1 80 0.060 0.958 +1 1 81 0.036 0.937 +1 1 82 0.037 0.945 +1 1 83 0.029 0.946 +1 1 84 0.033 0.941 +1 1 85 0.041 0.952 +1 1 86 0.033 0.899 +1 1 87 0.056 0.868 +1 1 88 0.019 0.902 +1 1 89 0.031 0.905 +1 1 90 0.041 0.939 +1 1 91 0.027 0.913 +1 1 92 0.036 0.937 +1 1 93 0.035 0.899 +1 1 94 0.011 0.906 +1 1 95 0.051 0.903 +1 1 96 0.022 0.875 +1 1 97 0.022 0.922 +1 1 98 0.018 0.938 +1 1 99 0.047 0.917 +1 1 100 0.039 0.887 +1 1 101 0.052 0.870 +1 1 102 0.038 0.934 +1 1 103 0.028 0.893 +1 1 104 0.045 0.889 +1 1 105 0.034 0.861 +1 1 106 0.018 0.898 +1 1 107 0.022 0.926 +1 1 108 0.025 0.915 +1 1 109 0.075 0.924 +1 1 110 0.023 0.899 +1 1 111 0.056 0.892 +1 1 112 0.029 0.900 +1 1 113 0.038 0.871 +1 1 114 0.027 0.888 +1 1 115 0.026 0.857 +1 1 116 0.064 0.893 +1 1 117 -0.018 0.865 +1 1 118 0.025 0.848 +1 1 119 0.002 0.888 +1 1 120 0.044 0.930 +1 1 121 0.013 0.913 +1 1 122 0.059 0.920 +1 1 123 0.066 0.880 +1 1 124 0.051 0.900 +1 1 125 0.033 0.888 +1 1 126 0.032 0.886 +1 1 127 0.021 0.904 +1 2 0 0.076 0.929 2.077 +1 2 1 0.053 0.958 2.098 +1 2 2 0.068 0.915 2.095 +1 2 3 0.066 0.957 2.101 +1 2 4 0.065 0.932 2.101 +1 2 5 0.069 0.874 2.120 +1 2 6 0.066 0.941 2.161 +1 2 7 0.044 0.899 2.065 +1 2 8 0.058 0.935 2.127 +1 2 9 0.066 0.922 2.121 +1 2 10 0.065 0.935 2.102 +1 2 11 0.072 0.936 2.108 +1 2 12 0.047 0.914 2.127 +1 2 13 0.051 0.934 2.110 +1 2 14 0.072 0.870 2.126 +1 2 15 0.086 0.926 2.105 +1 2 16 0.052 0.921 2.114 +1 2 17 0.055 0.912 2.252 +1 2 18 0.066 0.926 2.097 +1 2 19 0.057 0.919 2.087 +1 2 20 0.040 0.922 2.092 +1 2 21 0.042 0.913 2.103 +1 2 22 0.067 0.956 2.109 +1 2 23 0.056 0.907 2.089 +1 2 24 0.085 0.864 2.111 +1 2 25 0.059 0.908 2.095 +1 2 26 0.022 0.889 2.093 +1 2 27 0.058 0.904 2.081 +1 2 28 0.052 0.927 2.079 +1 2 29 0.010 0.912 2.097 +1 2 30 0.042 0.936 2.104 +1 2 31 0.051 0.732 2.098 +1 2 32 0.054 0.898 2.065 +1 2 33 0.086 0.874 2.108 +1 2 34 0.064 0.889 2.112 +1 2 35 0.064 0.882 2.085 +1 2 36 0.053 0.912 2.103 +1 2 37 0.052 0.905 2.079 +1 2 38 0.026 0.907 2.072 +1 2 39 0.068 0.888 2.109 +1 2 40 0.071 0.908 2.074 +1 2 41 0.081 0.936 2.083 +1 2 42 0.068 0.891 2.089 +1 2 43 0.020 0.886 2.086 +1 2 44 0.035 0.880 2.083 +1 2 45 0.040 0.901 2.088 +1 2 46 0.075 0.925 2.063 +1 2 47 0.012 0.906 2.109 +1 2 48 0.029 0.913 2.098 +1 2 49 0.048 0.865 2.065 +1 2 50 0.032 0.919 2.048 +1 2 51 0.037 0.900 2.062 +1 2 52 0.043 0.931 2.110 +1 2 53 0.047 0.923 2.068 +1 2 54 0.051 0.920 2.049 +1 2 55 0.057 0.894 2.078 +1 2 56 0.050 0.883 2.086 +1 2 57 0.034 0.894 2.087 +1 2 58 0.016 0.879 2.072 +1 2 59 0.058 0.874 2.071 +1 2 60 0.057 0.898 2.016 +1 2 61 0.047 0.902 2.085 +1 2 62 0.038 0.899 2.071 +1 2 63 0.077 0.876 2.064 +1 2 64 0.039 0.878 2.073 +1 2 65 0.090 0.912 2.082 +1 2 66 0.060 0.876 2.065 +1 2 67 0.050 0.925 2.029 +1 2 68 0.070 0.860 2.080 +1 2 69 0.029 0.888 2.040 +1 2 70 0.063 0.913 2.089 +1 2 71 0.015 0.876 2.087 +1 2 72 0.038 0.934 2.040 +1 2 73 0.053 0.910 2.071 +1 2 74 0.047 0.928 2.050 +1 2 75 0.061 0.921 2.080 +1 2 76 0.029 0.899 2.070 +1 2 77 0.033 0.899 2.046 +1 2 78 0.066 0.876 2.052 +1 2 79 0.092 0.902 2.051 +1 2 80 0.050 0.931 2.068 +1 2 81 0.024 0.893 2.064 +1 2 82 0.051 0.914 2.063 +1 2 83 0.031 0.888 2.048 +1 2 84 0.075 0.920 2.063 +1 2 85 0.045 0.886 2.060 +1 2 86 0.041 0.946 2.068 +1 2 87 0.053 0.893 2.062 +1 2 88 0.053 0.904 2.061 +1 2 89 0.029 0.901 2.046 +1 2 90 0.020 0.861 2.077 +1 2 91 -0.025 0.877 2.063 +1 2 92 0.027 0.897 2.065 +1 2 93 0.010 0.741 2.051 +1 2 94 0.017 0.888 2.088 +1 2 95 0.010 0.912 2.068 +1 2 96 0.043 0.871 2.024 +1 2 97 0.038 0.886 2.047 +1 2 98 0.018 0.866 2.052 +1 2 99 0.030 0.876 2.070 +1 2 100 0.046 0.840 2.057 +1 2 101 0.151 0.895 2.066 +1 2 102 0.032 0.878 2.061 +1 2 103 0.047 0.893 2.034 +1 2 104 0.048 0.893 2.064 +1 2 105 0.037 0.858 2.047 +1 2 106 0.046 0.891 2.064 +1 2 107 -0.006 0.907 2.083 +1 2 108 0.032 0.891 2.050 +1 2 109 0.049 0.906 2.047 +1 2 110 0.025 0.887 2.041 +1 2 111 0.017 0.919 2.053 +1 2 112 0.035 0.892 2.058 +1 2 113 0.033 0.880 2.056 +1 2 114 0.044 0.897 2.036 +1 2 115 0.032 0.869 2.041 +1 2 116 0.021 0.885 2.044 +1 2 117 -0.006 0.877 2.059 +1 2 118 0.037 0.919 2.063 +1 2 119 0.045 0.880 2.056 +1 2 120 0.041 0.887 2.032 +1 2 121 0.004 0.895 2.041 +1 2 122 0.034 0.873 2.042 +1 2 123 -0.006 0.891 2.024 +1 2 124 0.024 0.892 2.037 +1 2 125 0.036 0.899 2.043 +1 2 126 0.027 0.914 2.044 +1 2 127 0.040 0.881 2.042 +1 3 0 -0.022 0.852 2.002 +1 3 1 -0.036 0.851 2.008 +1 3 2 0.003 0.400 1.984 +1 3 3 -0.029 0.841 2.001 +1 3 4 -0.008 0.848 2.016 +1 3 5 0.018 0.863 2.016 +1 3 6 -0.022 0.870 1.990 +1 3 7 -0.009 0.822 1.999 +1 3 8 -0.012 0.848 2.008 +1 3 9 -0.045 0.851 2.017 +1 3 10 0.009 0.844 2.015 +1 3 11 -0.040 0.796 1.970 +1 3 12 -0.013 0.841 1.993 +1 3 13 -0.005 0.861 1.982 +1 3 14 -0.024 0.786 2.008 +1 3 15 -0.029 0.866 1.991 +1 3 16 -0.049 0.790 2.005 +1 3 17 -0.030 0.861 1.926 +1 3 18 -0.045 0.859 2.011 +1 3 19 -0.032 0.814 2.006 +1 3 20 -0.027 0.834 1.982 +1 3 21 0.012 0.836 1.991 +1 3 22 -0.074 0.794 2.005 +1 3 23 -0.031 0.848 1.988 +1 3 24 0.028 0.841 1.953 +1 3 25 -0.086 0.858 2.015 +1 3 26 -0.040 0.842 1.965 +1 3 27 -0.020 0.856 2.014 +1 3 28 -0.007 0.825 1.963 +1 3 29 -0.057 0.823 1.995 +1 3 30 -0.024 0.869 1.991 +1 3 31 -0.069 0.830 1.968 +1 3 32 0.009 0.773 1.990 +1 3 33 -0.057 0.832 1.979 +1 3 34 -0.073 0.854 2.007 +1 3 35 0.033 0.858 2.018 +1 3 36 0.004 0.815 2.009 +1 3 37 -0.004 0.861 2.020 +1 3 38 -0.126 0.835 1.985 +1 3 39 -0.020 0.820 2.024 +1 3 40 -0.039 0.852 1.943 +1 3 41 -0.055 0.846 1.990 +1 3 42 0.002 0.786 1.965 +1 3 43 -0.043 0.811 1.992 +1 3 44 -0.041 0.847 2.003 +1 3 45 0.000 0.829 1.985 +1 3 46 -0.024 0.844 1.972 +1 3 47 -0.056 0.819 1.963 +1 3 48 -0.021 0.812 1.895 +1 3 49 -0.045 0.819 1.980 +1 3 50 -0.047 0.831 1.947 +1 3 51 -0.035 0.798 1.974 +1 3 52 -0.008 0.866 1.964 +1 3 53 0.002 0.818 1.970 +1 3 54 -0.031 0.850 1.945 +1 3 55 -0.063 0.827 2.004 +1 3 56 -0.062 0.852 1.961 +1 3 57 -0.028 0.809 1.971 +1 3 58 -0.054 0.855 1.983 +1 3 59 -0.016 0.829 1.974 +1 3 60 -0.050 0.775 1.985 +1 3 61 -0.008 0.833 2.010 +1 3 62 -0.094 0.786 1.919 +1 3 63 -0.049 0.730 2.007 +1 3 64 -0.180 0.810 1.981 +1 3 65 -0.060 0.823 1.971 +1 3 66 -0.079 0.814 1.942 +1 3 67 -0.063 1.003 1.961 +1 3 68 -0.016 0.852 1.956 +1 3 69 0.827 1.986 +1 3 70 -0.008 0.791 1.956 +1 3 71 -0.059 0.838 1.976 +1 3 72 -0.075 0.838 2.008 +1 3 73 0.014 0.811 1.987 +1 3 74 -0.036 0.803 1.986 +1 3 75 -0.035 0.725 1.986 +1 3 76 -0.033 0.813 1.961 +1 3 77 -0.045 0.821 1.992 +1 3 78 -0.048 0.843 +1 3 79 0.005 0.848 1.941 +1 3 80 -0.085 0.829 1.972 +1 3 81 0.771 1.972 +1 3 82 -0.035 0.787 1.972 +1 3 83 -0.075 0.775 1.975 +1 3 84 -0.054 0.772 1.954 +1 3 85 -0.171 0.796 1.998 +1 3 86 -0.039 0.759 1.937 +1 3 87 -0.047 0.777 1.945 +1 3 88 -0.017 0.823 1.969 +1 3 89 0.787 1.960 +1 3 90 0.804 1.966 +1 3 91 0.766 1.947 +1 3 92 -0.025 0.818 1.962 +1 3 93 0.853 1.959 +1 3 94 0.793 1.982 +1 3 95 -0.094 0.800 1.986 +1 3 96 0.778 1.935 +1 3 97 0.831 1.993 +1 3 98 -0.013 0.867 1.944 +1 3 99 0.807 1.878 +1 3 100 0.116 1.895 +1 3 101 -0.026 0.802 1.897 +1 3 102 -0.070 0.735 1.996 +1 3 103 0.791 1.918 +1 3 104 0.767 1.941 +1 3 105 0.834 1.959 +1 3 106 0.812 2.017 +1 3 107 0.124 0.751 1.950 +1 3 108 0.786 1.927 +1 3 109 0.745 1.150 1.950 +1 3 110 -0.001 0.758 1.932 +1 3 111 -0.175 1.916 +1 3 112 0.812 1.965 +1 3 113 0.856 1.976 +1 3 114 0.854 1.905 +1 3 115 0.713 1.944 +1 3 116 0.825 1.975 +1 3 117 0.861 1.939 +1 3 118 0.775 1.924 +1 3 119 0.775 1.925 +1 3 120 0.825 1.925 +1 3 121 0.725 1.911 +1 3 122 0.923 1.865 +1 3 123 0.775 1.880 +1 3 124 0.883 1.884 +1 3 125 0.883 1.889 +1 3 126 0.775 1.950 +1 3 127 0.725 2.028 +2 1 30 0.003 +2 1 31 0.024 +2 1 32 0.026 +2 1 33 0.038 +2 1 34 0.034 +2 1 35 0.050 +2 1 36 0.060 +2 1 37 0.073 +2 1 38 0.019 +2 1 39 0.034 +2 1 40 0.051 +2 1 41 0.031 +2 1 42 0.030 +2 1 43 0.037 +2 1 44 0.067 +2 1 45 0.056 +2 1 46 -0.001 +2 1 47 0.047 +2 1 48 0.026 +2 1 49 0.052 +2 1 50 0.012 +2 1 51 0.004 +2 1 52 0.032 +2 1 53 0.021 +2 1 54 0.014 +2 1 55 0.041 +2 1 56 0.037 +2 1 57 0.012 +2 1 58 0.014 +2 1 59 0.023 +2 1 60 -0.038 +2 1 61 0.011 +2 1 62 0.032 +2 1 63 0.054 +2 1 64 0.029 +2 1 65 0.060 +2 1 66 0.019 +2 1 67 0.026 +2 1 68 0.030 +2 1 69 0.005 +2 1 70 0.015 0.925 +2 1 71 0.014 0.754 +2 1 72 0.013 0.879 +2 1 73 0.034 0.833 +2 1 74 0.013 0.972 +2 1 75 0.024 0.878 +2 1 76 0.027 0.890 +2 1 77 0.015 0.858 +2 1 78 0.032 0.920 +2 1 79 0.056 0.902 +2 1 80 -0.011 0.851 +2 1 81 0.017 0.863 +2 1 82 0.025 0.892 +2 1 83 0.059 0.907 +2 1 84 0.038 0.895 +2 1 85 0.005 0.853 +2 1 86 0.048 0.921 +2 1 87 0.032 0.893 +2 1 88 0.016 0.871 +2 1 89 -0.002 0.881 +2 1 90 0.017 0.922 +2 1 91 0.035 0.890 +2 1 92 -0.029 0.896 +2 1 93 -0.010 0.884 +2 1 94 0.005 0.885 +2 1 95 0.025 0.860 +2 1 96 0.064 0.893 +2 1 97 0.025 0.875 +2 1 98 0.012 0.888 +2 1 99 0.027 0.900 +2 1 100 0.009 0.886 +2 1 101 0.041 0.902 +2 1 102 0.020 0.946 +2 1 103 0.046 0.865 +2 1 104 0.040 0.914 +2 1 105 0.042 0.882 +2 1 106 0.029 0.874 +2 1 107 0.030 0.841 +2 1 108 0.041 0.857 +2 1 109 0.024 0.942 +2 1 110 0.014 0.843 +2 1 111 0.030 0.874 +2 1 112 0.022 0.903 +2 1 113 0.023 0.860 +2 1 114 0.025 0.910 +2 1 115 0.057 0.929 +2 1 116 0.062 0.900 +2 1 117 0.022 0.862 +2 1 118 0.024 0.844 +2 1 119 0.045 0.870 +2 2 0 0.004 0.858 2.008 +2 2 1 0.025 0.842 1.941 +2 2 2 0.006 0.843 1.998 +2 2 3 -0.020 0.851 1.958 +2 2 4 -0.007 0.834 1.918 +2 2 5 -0.020 0.860 2.035 +2 2 6 0.004 0.848 2.004 +2 2 7 -0.003 0.835 2.029 +2 2 8 -0.011 0.813 1.981 +2 2 9 -0.014 0.867 2.007 +2 2 10 -0.010 0.857 2.022 +2 2 11 -0.009 0.861 2.020 +2 2 12 -0.018 0.848 2.011 +2 2 13 0.019 0.872 2.021 +2 2 14 0.048 0.855 2.002 +2 2 15 0.003 0.838 2.000 +2 2 16 -0.001 0.866 2.028 +2 2 17 -0.003 0.834 2.000 +2 2 18 0.002 0.849 1.996 +2 2 19 -0.005 0.841 2.011 +2 2 20 -0.005 0.857 2.030 +2 2 21 -0.007 0.843 2.018 +2 2 22 -0.027 0.821 2.017 +2 2 23 -0.018 0.846 2.031 +2 2 24 0.007 0.792 2.008 +2 2 25 -0.010 0.846 2.018 +2 2 26 -0.026 0.843 1.989 +2 2 27 -0.007 0.850 2.006 +2 2 28 -0.007 0.828 2.005 +2 2 29 -0.005 0.848 1.996 +2 2 30 -0.001 0.833 2.022 +2 2 31 0.007 0.842 2.004 +2 2 32 -0.002 0.816 2.026 +2 2 33 -0.013 0.839 1.997 +2 2 34 -0.020 0.913 2.025 +2 2 35 -0.010 0.870 2.011 +2 2 36 -0.001 0.854 2.006 +2 2 37 -0.010 0.828 1.957 +2 2 38 -0.001 0.827 1.998 +2 2 39 0.007 0.822 1.999 +2 2 40 -0.005 0.818 2.020 +2 2 41 -0.011 0.846 2.016 +2 2 42 0.022 0.875 2.016 +2 2 43 -0.021 0.846 2.009 +2 2 44 -0.014 0.831 1.988 +2 2 45 -0.004 0.748 2.008 +2 2 46 -0.012 0.834 2.003 +2 2 47 -0.014 0.851 2.017 +2 2 48 -0.001 0.838 2.016 +2 2 49 0.000 0.818 1.993 +2 2 50 0.021 0.800 2.008 +2 2 51 -0.005 0.800 2.022 +2 2 52 -0.012 0.828 2.024 +2 2 53 -0.025 0.825 2.025 +2 2 54 -0.007 0.842 2.008 +2 2 55 0.010 0.836 2.012 +2 2 56 -0.003 0.849 1.989 +2 2 57 -0.007 0.843 2.027 +2 2 58 -0.025 0.825 1.993 +2 2 59 0.002 0.845 2.002 +2 2 60 -0.017 0.835 1.992 +2 2 61 -0.001 0.851 1.998 +2 2 62 -0.020 0.814 2.021 +2 2 63 -0.025 0.852 2.007 +2 2 64 -0.034 0.856 2.016 +2 2 65 -0.005 0.841 1.989 +2 2 66 -0.015 0.815 2.018 +2 2 67 -0.015 0.793 1.984 +2 2 68 -0.011 0.840 2.018 +2 2 69 -0.013 0.846 1.995 +2 2 70 -0.009 0.838 2.025 +2 2 71 -0.021 0.812 2.008 +2 2 72 -0.016 0.844 2.005 +2 2 73 0.002 0.837 1.978 +2 2 74 0.036 0.860 2.024 +2 2 75 -0.013 0.859 2.002 +2 2 76 -0.008 0.870 1.990 +2 2 77 0.009 0.825 1.994 +2 2 78 0.029 0.864 1.995 +2 2 79 -0.011 0.842 1.993 +2 2 80 -0.001 0.826 1.973 +2 2 81 -0.013 0.853 1.996 +2 2 82 -0.003 0.817 1.980 +2 2 83 0.008 0.854 1.975 +2 2 84 -0.032 0.862 2.018 +2 2 85 0.005 0.852 2.026 +2 2 86 0.001 0.804 1.962 +2 2 87 -0.016 0.844 2.018 +2 2 88 0.003 0.843 1.976 +2 2 89 0.006 0.813 1.991 +2 2 90 0.027 0.865 2.037 +2 2 91 -0.000 0.861 2.011 +2 2 92 -0.026 0.835 2.000 +2 2 93 0.013 0.867 2.010 +2 2 94 0.022 0.861 2.019 +2 2 95 -0.020 0.850 2.017 +2 2 96 -0.005 0.858 1.997 +2 2 97 -0.000 0.852 2.018 +2 2 98 -0.025 0.853 1.982 +2 2 99 -0.013 0.844 2.017 +2 2 100 -0.008 0.702 1.997 +2 2 101 -0.004 0.855 2.004 +2 2 102 -0.023 0.851 1.970 +2 2 103 -0.018 0.860 1.982 +2 2 104 -0.036 0.803 1.963 +2 2 105 -0.017 0.754 2.006 +2 2 106 -0.003 0.832 1.985 +2 2 107 -0.020 0.836 1.985 +2 2 108 -0.001 0.834 2.000 +2 2 109 -0.015 0.854 1.977 +2 2 110 -0.006 0.834 1.984 +2 2 111 0.021 0.817 2.022 +2 2 112 -0.004 0.801 1.986 +2 2 113 0.004 0.832 1.997 +2 2 114 -0.041 0.854 1.980 +2 2 115 -0.029 0.830 1.981 +2 2 116 0.019 0.865 1.978 +2 2 117 0.011 0.829 1.946 +2 2 118 -0.017 0.855 1.972 +2 2 119 -0.039 0.808 1.999 +2 2 120 -0.022 0.859 1.993 +2 2 121 0.033 0.857 2.003 +2 2 122 0.004 0.834 1.993 +2 2 123 -0.007 0.837 2.003 +2 2 124 -0.037 0.821 2.000 +2 2 125 0.006 0.812 1.998 +2 2 126 -0.028 0.836 2.020 +2 2 127 -0.020 0.825 1.943 +2 3 0 -0.029 0.853 1.992 +2 3 1 -0.035 0.838 1.984 +2 3 2 0.014 0.864 2.019 +2 3 3 -0.006 0.868 1.964 +2 3 4 -0.026 0.844 1.989 +2 3 5 -0.010 0.825 1.991 +2 3 6 0.057 0.704 1.997 +2 3 7 0.015 0.825 2.002 +2 3 8 -0.016 0.672 1.984 +2 3 9 -0.038 0.840 1.966 +2 3 10 -0.023 0.851 2.014 +2 3 11 -0.023 0.822 1.977 +2 3 12 -0.059 0.872 1.972 +2 3 13 -0.038 0.834 1.985 +2 3 14 -0.057 0.849 1.984 +2 3 15 -0.025 0.852 1.986 +2 3 16 -0.025 0.831 1.997 +2 3 17 0.014 0.830 1.981 +2 3 18 -0.008 0.812 2.013 +2 3 19 -0.010 0.833 1.981 +2 3 20 -0.004 0.855 2.008 +2 3 21 0.818 1.989 +2 3 22 -0.029 0.792 1.985 +2 3 23 -0.027 0.824 1.985 +2 3 24 -0.011 0.875 2.016 +2 3 25 -0.035 0.839 1.985 +2 3 26 -0.036 0.825 1.955 +2 3 27 -0.035 0.804 1.997 +2 3 28 -0.011 0.867 1.976 +2 3 29 -0.021 0.842 2.004 +2 3 30 0.819 1.987 +2 3 31 -0.030 0.832 1.974 +2 3 32 -0.032 0.700 2.000 +2 3 33 -0.031 0.809 1.951 +2 3 34 -0.058 0.841 1.975 +2 3 35 0.002 0.796 1.951 +2 3 36 -0.047 0.810 1.991 +2 3 37 -0.035 0.876 2.008 +2 3 38 -0.029 0.811 2.008 +2 3 39 0.002 0.830 1.939 +2 3 40 0.048 0.844 1.972 +2 3 41 -0.041 0.837 2.006 +2 3 42 -0.045 0.808 1.981 +2 3 43 0.022 0.846 1.977 +2 3 44 -0.034 0.834 2.000 +2 3 45 -0.021 0.842 1.980 +2 3 46 0.825 1.980 +2 3 47 -0.069 0.815 1.997 +2 3 48 -0.031 0.841 1.965 +2 3 49 -0.049 0.822 1.976 +2 3 50 0.009 2.013 +2 3 51 -0.037 0.814 1.957 +2 3 52 -0.033 0.827 1.945 +2 3 53 -0.034 0.798 1.974 +2 3 54 0.005 0.832 1.934 +2 3 55 0.795 2.012 +2 3 56 0.020 0.826 1.964 +2 3 57 -0.045 0.830 1.956 +2 3 58 -0.086 0.806 1.959 +2 3 59 -0.043 0.844 1.977 +2 3 60 -0.013 0.818 1.960 +2 3 61 -0.023 0.797 1.997 +2 3 62 -0.108 0.793 2.034 +2 3 63 -0.039 0.786 1.970 +2 3 64 -0.053 0.816 1.978 +2 3 65 -0.069 0.845 1.981 +2 3 66 -0.049 0.851 1.951 +2 3 67 -0.021 0.852 2.047 +2 3 68 0.002 0.864 1.962 +2 3 69 -0.067 0.777 1.930 +2 3 70 -0.011 0.790 1.922 +2 3 71 -0.049 0.800 1.974 +2 3 72 -0.032 0.827 1.971 +2 3 73 -0.030 0.786 1.967 +2 3 74 -0.037 0.790 1.917 +2 3 75 0.831 1.964 +2 3 76 -0.056 0.799 1.914 +2 3 77 -0.039 0.770 1.986 +2 3 78 -0.051 0.773 1.958 +2 3 79 -0.086 0.822 1.933 +2 3 80 0.794 1.943 +2 3 81 -0.083 0.801 2.023 +2 3 82 -0.030 0.783 1.948 +2 3 83 -0.059 0.686 1.983 +2 3 84 -0.079 0.773 1.950 +2 3 85 -0.018 0.782 1.940 +2 3 86 -0.083 0.814 1.969 +2 3 87 -0.042 0.804 1.945 +2 3 88 0.030 0.828 1.928 +2 3 89 -0.102 0.811 1.938 +2 3 90 0.745 1.954 +2 3 91 -0.043 0.817 1.979 +2 3 92 -0.122 0.765 1.950 +2 3 93 -0.075 0.800 1.937 +2 3 94 -0.014 0.788 1.930 +2 3 95 -0.048 0.744 1.896 +2 3 96 -0.045 0.826 1.963 +2 3 97 -0.105 0.851 1.953 +2 3 98 -0.008 0.751 1.961 +2 3 99 -0.101 0.807 1.926 +2 3 100 -0.037 0.847 1.947 +2 3 101 -0.068 0.766 1.920 +2 3 102 0.796 1.940 +2 3 103 0.023 0.795 1.957 +2 3 104 -0.087 0.788 1.921 +2 3 105 -0.099 0.811 1.968 +2 3 106 -0.059 0.771 1.971 +2 3 107 -0.125 0.808 1.928 +2 3 108 0.798 1.938 +2 3 109 -0.075 0.880 2.000 +2 3 110 -0.066 0.785 1.908 +2 3 111 -0.075 0.823 1.999 +2 3 112 -0.050 0.744 1.958 +2 3 113 -0.073 0.712 1.970 +2 3 114 -0.020 0.795 1.977 +2 3 115 -0.091 0.789 1.969 +2 3 116 -0.140 0.904 1.925 +2 3 117 -0.055 0.796 1.925 +2 3 118 1.923 +2 3 119 -0.148 0.743 1.925 +2 3 120 -0.075 0.755 1.974 +2 3 121 0.000 0.775 1.925 +2 3 122 -0.075 0.725 2.029 +2 3 123 -0.029 0.642 1.901 +2 3 124 -0.117 0.846 1.863 +2 3 125 -0.025 0.812 1.975 +2 3 126 -0.025 0.745 1.897 +2 3 127 -0.127 1.856 \ No newline at end of file diff --git a/getPeaks.cpp b/getPeaks.cpp new file mode 100644 index 0000000..07353ba --- /dev/null +++ b/getPeaks.cpp @@ -0,0 +1,130 @@ +#include "TFile.h" +#include "TTree.h" +#include "TH1F.h" +#include "TStyle.h" +#include "TH2F.h" +#include "TCanvas.h" +#include "TTreeReader.h" + +#include "AutoFit.C" + +// mod = 0, 1, 2 +// row = 1, 2, 3 +// pstrip = 0 - 127 + +#define MAXMOD 3 +#define MAXROW 3 +#define MAXPStrip 128 + +const double ExRange[2] = {-0.5, 2.5}; +const double ExRel = 0.05; //MeV + + +TH1F * haha[3][3][128]; + +TH2F * dudu[3][3]; + +TH1F * hEx; + +void getPeaks(){ + + + TFile * file = new TFile("R9-172_summed_tree.root"); + + TTree * tree = (TTree *) file->Get("rxtree"); + + const int nBin = (ExRange[1] - ExRange[0])/ExRel; + + for( int i = 0; i < MAXMOD; i++ ){ + for( int j = 0; j < MAXROW; j++ ){ + + dudu[i][j] = new TH2F(Form("dudu_%d%d", i, j+1), Form("%d-%d; pStrip; Ex", i, j+1), 128, 0, 128, nBin, ExRange[0], ExRange[1]); + + for( int p = 0; p < MAXPStrip; p ++ ){ + haha[i][j][p] = new TH1F(Form("h_%d%d_%d", i, j+1, p), Form("%d-%d-%d", i, j+1, p), nBin, ExRange[0], ExRange[1]); + } + } + } + + hEx = new TH1F("hEx", "Ex", 10*nBin, ExRange[0], ExRange[1]); + + TTreeReader reader(tree); + + TTreeReaderValue Ex = {reader, "Ex"}; + TTreeReaderValue mod = {reader, "mod"}; + TTreeReaderValue row = {reader, "row"}; + TTreeReaderValue pStrip = {reader, "pstrip"}; + TTreeReaderValue ebis = {reader, "ebis_on"}; + + while( reader.Next() ){ + + if( *ebis == 0 ) continue; + if( *row == 0 ) continue; + + int m = *mod; + int r = *row -1; + int p = *pStrip; + + // printf("%d %d %d | %f\n", m, r, p, *Ex/1000.); + + haha[m][r][p]->Fill(*Ex/1000.); + + dudu[m][r]->Fill(p, *Ex/1000); + + hEx->Fill( *Ex/1000); + + } + + gStyle->SetOptStat(""); + + // TCanvas * canvas = new TCanvas("canvas", "canvas", 2100, 1500); + // haha[0][2][10]->Draw(); + + + // canvas->Divide(3,3); + + // canvas->cd(1); dudu[0][0]->Draw("colz"); + // canvas->cd(2); dudu[0][1]->Draw("colz"); + // canvas->cd(3); dudu[0][2]->Draw("colz"); + + // canvas->cd(4); dudu[1][0]->Draw("colz"); + // canvas->cd(5); dudu[1][1]->Draw("colz"); + // canvas->cd(6); dudu[1][2]->Draw("colz"); + + // canvas->cd(7); dudu[2][0]->Draw("colz"); + // canvas->cd(8); dudu[2][1]->Draw("colz"); + // canvas->cd(9); dudu[2][2]->Draw("colz"); + + // hEx->Draw(); + + // FILE * outList = fopen("fit_result.txt", "w"); + + // for( int i = 0; i < MAXMOD; i ++ ){ + // printf("############################## mod:%d\n", i); + // for( int j = 0; j < MAXROW; j++ ){ + // printf("============================== row:%d\n",j); + // for( int k = 0; k < MAXPStrip; k ++ ){ + // if( j == 0 && k < 30 ) continue; + + // std::vector exList = fitAuto(haha[i][j][k], -1, 0.1, 0.13, 3); + + // fprintf(outList, "%d %d %3d", i, j+1, k); + // for( size_t p = 0; p < exList.size(); p++ ) { + // if( 1.8 > exList[p] && exList[p] > 1.2 ) continue; + // fprintf(outList, " %6.3f", exList[p]); + // } + // fprintf(outList, "\n"); + + // // printf("%3d |", k); + // // for( size_t p = 0; p < exList.size(); p++ ) printf(" %.3f", exList[p]); + // // printf("\n"); + + // } + // } + // } + + // fclose(outList); + + + +} \ No newline at end of file diff --git a/raw_ex_50keV.txt b/raw_ex_50keV.txt new file mode 100644 index 0000000..435373c --- /dev/null +++ b/raw_ex_50keV.txt @@ -0,0 +1,800 @@ + -0.49750000 52.000000 + -0.49250000 61.000000 + -0.48750000 51.000000 + -0.48250000 56.000000 + -0.47750000 42.000000 + -0.47250000 49.000000 + -0.46750000 61.000000 + -0.46250000 63.000000 + -0.45750000 43.000000 + -0.45250000 54.000000 + -0.44750000 52.000000 + -0.44250000 56.000000 + -0.43750000 37.000000 + -0.43250000 44.000000 + -0.42750000 53.000000 + -0.42250000 62.000000 + -0.41750000 51.000000 + -0.41250000 44.000000 + -0.40750000 53.000000 + -0.40250000 61.000000 + -0.39750000 42.000000 + -0.39250000 52.000000 + -0.38750000 65.000000 + -0.38250000 56.000000 + -0.37750000 54.000000 + -0.37250000 72.000000 + -0.36750000 52.000000 + -0.36250000 56.000000 + -0.35750000 44.000000 + -0.35250000 49.000000 + -0.34750000 67.000000 + -0.34250000 44.000000 + -0.33750000 51.000000 + -0.33250000 57.000000 + -0.32750000 70.000000 + -0.32250000 51.000000 + -0.31750000 69.000000 + -0.31250000 63.000000 + -0.30750000 54.000000 + -0.30250000 61.000000 + -0.29750000 68.000000 + -0.29250000 66.000000 + -0.28750000 57.000000 + -0.28250000 67.000000 + -0.27750000 46.000000 + -0.27250000 56.000000 + -0.26750000 58.000000 + -0.26250000 54.000000 + -0.25750000 60.000000 + -0.25250000 56.000000 + -0.24750000 63.000000 + -0.24250000 68.000000 + -0.23750000 87.000000 + -0.23250000 55.000000 + -0.22750000 59.000000 + -0.22250000 63.000000 + -0.21750000 71.000000 + -0.21250000 55.000000 + -0.20750000 69.000000 + -0.20250000 84.000000 + -0.19750000 80.000000 + -0.19250000 79.000000 + -0.18750000 100.000000 + -0.18250000 91.000000 + -0.17750000 103.000000 + -0.17250000 121.000000 + -0.16750000 106.000000 + -0.16250000 120.000000 + -0.15750000 129.000000 + -0.15250000 137.000000 + -0.14750000 128.000000 + -0.14250000 169.000000 + -0.13750000 174.000000 + -0.13250000 193.000000 + -0.12750000 247.000000 + -0.12250000 231.000000 + -0.11750000 251.000000 + -0.11250000 300.000000 + -0.10750000 330.000000 + -0.10250000 355.000000 + -0.09750000 421.000000 + -0.09250000 471.000000 + -0.08750000 460.000000 + -0.08250000 531.000000 + -0.07750000 579.000000 + -0.07250000 558.000000 + -0.06750000 645.000000 + -0.06250000 674.000000 + -0.05750000 753.000000 + -0.05250000 749.000000 + -0.04750000 770.000000 + -0.04250000 889.000000 + -0.03750000 878.000000 + -0.03250000 888.000000 + -0.02750000 916.000000 + -0.02250000 917.000000 + -0.01750000 935.000000 + -0.01250000 1018.000000 + -0.00750000 1009.000000 + -0.00250000 1043.000000 + 0.00250000 1024.000000 + 0.00750000 977.000000 + 0.01250000 963.000000 + 0.01750000 1005.000000 + 0.02250000 963.000000 + 0.02750000 942.000000 + 0.03250000 917.000000 + 0.03750000 891.000000 + 0.04250000 844.000000 + 0.04750000 809.000000 + 0.05250000 788.000000 + 0.05750000 818.000000 + 0.06250000 772.000000 + 0.06750000 748.000000 + 0.07250000 696.000000 + 0.07750000 658.000000 + 0.08250000 639.000000 + 0.08750000 604.000000 + 0.09250000 542.000000 + 0.09750000 553.000000 + 0.10250000 533.000000 + 0.10750000 501.000000 + 0.11250000 436.000000 + 0.11750000 418.000000 + 0.12250000 396.000000 + 0.12750000 419.000000 + 0.13250000 343.000000 + 0.13750000 342.000000 + 0.14250000 301.000000 + 0.14750000 294.000000 + 0.15250000 268.000000 + 0.15750000 240.000000 + 0.16250000 244.000000 + 0.16750000 212.000000 + 0.17250000 208.000000 + 0.17750000 194.000000 + 0.18250000 185.000000 + 0.18750000 175.000000 + 0.19250000 151.000000 + 0.19750000 147.000000 + 0.20250000 132.000000 + 0.20750000 123.000000 + 0.21250000 117.000000 + 0.21750000 143.000000 + 0.22250000 120.000000 + 0.22750000 130.000000 + 0.23250000 96.000000 + 0.23750000 110.000000 + 0.24250000 98.000000 + 0.24750000 111.000000 + 0.25250000 89.000000 + 0.25750000 83.000000 + 0.26250000 109.000000 + 0.26750000 105.000000 + 0.27250000 87.000000 + 0.27750000 123.000000 + 0.28250000 112.000000 + 0.28750000 95.000000 + 0.29250000 81.000000 + 0.29750000 105.000000 + 0.30250000 77.000000 + 0.30750000 91.000000 + 0.31250000 91.000000 + 0.31750000 75.000000 + 0.32250000 99.000000 + 0.32750000 80.000000 + 0.33250000 90.000000 + 0.33750000 89.000000 + 0.34250000 91.000000 + 0.34750000 97.000000 + 0.35250000 84.000000 + 0.35750000 85.000000 + 0.36250000 92.000000 + 0.36750000 98.000000 + 0.37250000 93.000000 + 0.37750000 89.000000 + 0.38250000 88.000000 + 0.38750000 70.000000 + 0.39250000 79.000000 + 0.39750000 85.000000 + 0.40250000 92.000000 + 0.40750000 78.000000 + 0.41250000 78.000000 + 0.41750000 94.000000 + 0.42250000 58.000000 + 0.42750000 88.000000 + 0.43250000 92.000000 + 0.43750000 84.000000 + 0.44250000 92.000000 + 0.44750000 84.000000 + 0.45250000 97.000000 + 0.45750000 72.000000 + 0.46250000 105.000000 + 0.46750000 84.000000 + 0.47250000 80.000000 + 0.47750000 94.000000 + 0.48250000 73.000000 + 0.48750000 78.000000 + 0.49250000 97.000000 + 0.49750000 70.000000 + 0.50250000 94.000000 + 0.50750000 93.000000 + 0.51250000 91.000000 + 0.51750000 79.000000 + 0.52250000 81.000000 + 0.52750000 87.000000 + 0.53250000 98.000000 + 0.53750000 94.000000 + 0.54250000 87.000000 + 0.54750000 102.000000 + 0.55250000 99.000000 + 0.55750000 81.000000 + 0.56250000 94.000000 + 0.56750000 112.000000 + 0.57250000 75.000000 + 0.57750000 93.000000 + 0.58250000 77.000000 + 0.58750000 77.000000 + 0.59250000 89.000000 + 0.59750000 92.000000 + 0.60250000 96.000000 + 0.60750000 99.000000 + 0.61250000 115.000000 + 0.61750000 92.000000 + 0.62250000 98.000000 + 0.62750000 93.000000 + 0.63250000 108.000000 + 0.63750000 102.000000 + 0.64250000 95.000000 + 0.64750000 101.000000 + 0.65250000 116.000000 + 0.65750000 115.000000 + 0.66250000 120.000000 + 0.66750000 134.000000 + 0.67250000 136.000000 + 0.67750000 134.000000 + 0.68250000 133.000000 + 0.68750000 162.000000 + 0.69250000 191.000000 + 0.69750000 175.000000 + 0.70250000 203.000000 + 0.70750000 227.000000 + 0.71250000 261.000000 + 0.71750000 297.000000 + 0.72250000 290.000000 + 0.72750000 347.000000 + 0.73250000 391.000000 + 0.73750000 377.000000 + 0.74250000 430.000000 + 0.74750000 462.000000 + 0.75250000 521.000000 + 0.75750000 513.000000 + 0.76250000 577.000000 + 0.76750000 599.000000 + 0.77250000 626.000000 + 0.77750000 658.000000 + 0.78250000 700.000000 + 0.78750000 747.000000 + 0.79250000 779.000000 + 0.79750000 764.000000 + 0.80250000 819.000000 + 0.80750000 822.000000 + 0.81250000 783.000000 + 0.81750000 843.000000 + 0.82250000 826.000000 + 0.82750000 850.000000 + 0.83250000 835.000000 + 0.83750000 912.000000 + 0.84250000 802.000000 + 0.84750000 885.000000 + 0.85250000 875.000000 + 0.85750000 837.000000 + 0.86250000 885.000000 + 0.86750000 857.000000 + 0.87250000 816.000000 + 0.87750000 802.000000 + 0.88250000 786.000000 + 0.88750000 794.000000 + 0.89250000 764.000000 + 0.89750000 712.000000 + 0.90250000 684.000000 + 0.90750000 718.000000 + 0.91250000 623.000000 + 0.91750000 669.000000 + 0.92250000 635.000000 + 0.92750000 583.000000 + 0.93250000 605.000000 + 0.93750000 518.000000 + 0.94250000 530.000000 + 0.94750000 493.000000 + 0.95250000 454.000000 + 0.95750000 458.000000 + 0.96250000 446.000000 + 0.96750000 379.000000 + 0.97250000 340.000000 + 0.97750000 352.000000 + 0.98250000 348.000000 + 0.98750000 313.000000 + 0.99250000 333.000000 + 0.99750000 292.000000 + 1.00250000 283.000000 + 1.00750000 260.000000 + 1.01250000 227.000000 + 1.01750000 233.000000 + 1.02250000 222.000000 + 1.02750000 217.000000 + 1.03250000 187.000000 + 1.03750000 189.000000 + 1.04250000 166.000000 + 1.04750000 172.000000 + 1.05250000 156.000000 + 1.05750000 164.000000 + 1.06250000 148.000000 + 1.06750000 149.000000 + 1.07250000 129.000000 + 1.07750000 131.000000 + 1.08250000 138.000000 + 1.08750000 135.000000 + 1.09250000 130.000000 + 1.09750000 130.000000 + 1.10250000 117.000000 + 1.10750000 104.000000 + 1.11250000 98.000000 + 1.11750000 122.000000 + 1.12250000 113.000000 + 1.12750000 116.000000 + 1.13250000 131.000000 + 1.13750000 101.000000 + 1.14250000 103.000000 + 1.14750000 109.000000 + 1.15250000 107.000000 + 1.15750000 106.000000 + 1.16250000 129.000000 + 1.16750000 108.000000 + 1.17250000 113.000000 + 1.17750000 125.000000 + 1.18250000 115.000000 + 1.18750000 136.000000 + 1.19250000 131.000000 + 1.19750000 137.000000 + 1.20250000 152.000000 + 1.20750000 140.000000 + 1.21250000 153.000000 + 1.21750000 187.000000 + 1.22250000 198.000000 + 1.22750000 190.000000 + 1.23250000 202.000000 + 1.23750000 222.000000 + 1.24250000 205.000000 + 1.24750000 254.000000 + 1.25250000 286.000000 + 1.25750000 268.000000 + 1.26250000 268.000000 + 1.26750000 306.000000 + 1.27250000 359.000000 + 1.27750000 353.000000 + 1.28250000 363.000000 + 1.28750000 391.000000 + 1.29250000 350.000000 + 1.29750000 414.000000 + 1.30250000 412.000000 + 1.30750000 417.000000 + 1.31250000 431.000000 + 1.31750000 416.000000 + 1.32250000 407.000000 + 1.32750000 421.000000 + 1.33250000 462.000000 + 1.33750000 438.000000 + 1.34250000 403.000000 + 1.34750000 417.000000 + 1.35250000 427.000000 + 1.35750000 458.000000 + 1.36250000 430.000000 + 1.36750000 456.000000 + 1.37250000 436.000000 + 1.37750000 437.000000 + 1.38250000 421.000000 + 1.38750000 410.000000 + 1.39250000 391.000000 + 1.39750000 403.000000 + 1.40250000 384.000000 + 1.40750000 396.000000 + 1.41250000 383.000000 + 1.41750000 392.000000 + 1.42250000 368.000000 + 1.42750000 359.000000 + 1.43250000 351.000000 + 1.43750000 358.000000 + 1.44250000 359.000000 + 1.44750000 357.000000 + 1.45250000 344.000000 + 1.45750000 342.000000 + 1.46250000 355.000000 + 1.46750000 315.000000 + 1.47250000 344.000000 + 1.47750000 312.000000 + 1.48250000 361.000000 + 1.48750000 272.000000 + 1.49250000 309.000000 + 1.49750000 323.000000 + 1.50250000 285.000000 + 1.50750000 303.000000 + 1.51250000 282.000000 + 1.51750000 272.000000 + 1.52250000 241.000000 + 1.52750000 259.000000 + 1.53250000 252.000000 + 1.53750000 247.000000 + 1.54250000 255.000000 + 1.54750000 234.000000 + 1.55250000 248.000000 + 1.55750000 227.000000 + 1.56250000 215.000000 + 1.56750000 231.000000 + 1.57250000 212.000000 + 1.57750000 225.000000 + 1.58250000 201.000000 + 1.58750000 193.000000 + 1.59250000 196.000000 + 1.59750000 198.000000 + 1.60250000 177.000000 + 1.60750000 180.000000 + 1.61250000 167.000000 + 1.61750000 161.000000 + 1.62250000 156.000000 + 1.62750000 190.000000 + 1.63250000 168.000000 + 1.63750000 169.000000 + 1.64250000 137.000000 + 1.64750000 161.000000 + 1.65250000 162.000000 + 1.65750000 172.000000 + 1.66250000 157.000000 + 1.66750000 143.000000 + 1.67250000 137.000000 + 1.67750000 119.000000 + 1.68250000 126.000000 + 1.68750000 129.000000 + 1.69250000 114.000000 + 1.69750000 132.000000 + 1.70250000 122.000000 + 1.70750000 132.000000 + 1.71250000 113.000000 + 1.71750000 111.000000 + 1.72250000 110.000000 + 1.72750000 121.000000 + 1.73250000 120.000000 + 1.73750000 110.000000 + 1.74250000 126.000000 + 1.74750000 120.000000 + 1.75250000 111.000000 + 1.75750000 116.000000 + 1.76250000 111.000000 + 1.76750000 107.000000 + 1.77250000 111.000000 + 1.77750000 139.000000 + 1.78250000 118.000000 + 1.78750000 101.000000 + 1.79250000 121.000000 + 1.79750000 143.000000 + 1.80250000 152.000000 + 1.80750000 136.000000 + 1.81250000 141.000000 + 1.81750000 162.000000 + 1.82250000 169.000000 + 1.82750000 181.000000 + 1.83250000 219.000000 + 1.83750000 240.000000 + 1.84250000 246.000000 + 1.84750000 287.000000 + 1.85250000 281.000000 + 1.85750000 293.000000 + 1.86250000 355.000000 + 1.86750000 401.000000 + 1.87250000 469.000000 + 1.87750000 480.000000 + 1.88250000 515.000000 + 1.88750000 590.000000 + 1.89250000 628.000000 + 1.89750000 661.000000 + 1.90250000 723.000000 + 1.90750000 801.000000 + 1.91250000 812.000000 + 1.91750000 877.000000 + 1.92250000 848.000000 + 1.92750000 887.000000 + 1.93250000 925.000000 + 1.93750000 920.000000 + 1.94250000 882.000000 + 1.94750000 906.000000 + 1.95250000 914.000000 + 1.95750000 989.000000 + 1.96250000 1009.000000 + 1.96750000 977.000000 + 1.97250000 997.000000 + 1.97750000 1006.000000 + 1.98250000 990.000000 + 1.98750000 1048.000000 + 1.99250000 983.000000 + 1.99750000 1058.000000 + 2.00250000 1007.000000 + 2.00750000 984.000000 + 2.01250000 987.000000 + 2.01750000 924.000000 + 2.02250000 1035.000000 + 2.02750000 899.000000 + 2.03250000 973.000000 + 2.03750000 861.000000 + 2.04250000 892.000000 + 2.04750000 844.000000 + 2.05250000 837.000000 + 2.05750000 856.000000 + 2.06250000 816.000000 + 2.06750000 708.000000 + 2.07250000 738.000000 + 2.07750000 706.000000 + 2.08250000 677.000000 + 2.08750000 673.000000 + 2.09250000 656.000000 + 2.09750000 604.000000 + 2.10250000 561.000000 + 2.10750000 531.000000 + 2.11250000 539.000000 + 2.11750000 538.000000 + 2.12250000 498.000000 + 2.12750000 489.000000 + 2.13250000 437.000000 + 2.13750000 403.000000 + 2.14250000 403.000000 + 2.14750000 384.000000 + 2.15250000 353.000000 + 2.15750000 347.000000 + 2.16250000 332.000000 + 2.16750000 322.000000 + 2.17250000 293.000000 + 2.17750000 270.000000 + 2.18250000 266.000000 + 2.18750000 260.000000 + 2.19250000 271.000000 + 2.19750000 244.000000 + 2.20250000 239.000000 + 2.20750000 205.000000 + 2.21250000 174.000000 + 2.21750000 186.000000 + 2.22250000 170.000000 + 2.22750000 171.000000 + 2.23250000 173.000000 + 2.23750000 150.000000 + 2.24250000 156.000000 + 2.24750000 162.000000 + 2.25250000 156.000000 + 2.25750000 142.000000 + 2.26250000 161.000000 + 2.26750000 126.000000 + 2.27250000 140.000000 + 2.27750000 129.000000 + 2.28250000 108.000000 + 2.28750000 145.000000 + 2.29250000 115.000000 + 2.29750000 113.000000 + 2.30250000 119.000000 + 2.30750000 104.000000 + 2.31250000 95.000000 + 2.31750000 116.000000 + 2.32250000 123.000000 + 2.32750000 104.000000 + 2.33250000 100.000000 + 2.33750000 117.000000 + 2.34250000 109.000000 + 2.34750000 103.000000 + 2.35250000 97.000000 + 2.35750000 102.000000 + 2.36250000 81.000000 + 2.36750000 106.000000 + 2.37250000 113.000000 + 2.37750000 69.000000 + 2.38250000 104.000000 + 2.38750000 96.000000 + 2.39250000 112.000000 + 2.39750000 97.000000 + 2.40250000 94.000000 + 2.40750000 101.000000 + 2.41250000 99.000000 + 2.41750000 88.000000 + 2.42250000 97.000000 + 2.42750000 109.000000 + 2.43250000 92.000000 + 2.43750000 102.000000 + 2.44250000 89.000000 + 2.44750000 89.000000 + 2.45250000 81.000000 + 2.45750000 88.000000 + 2.46250000 89.000000 + 2.46750000 101.000000 + 2.47250000 85.000000 + 2.47750000 93.000000 + 2.48250000 87.000000 + 2.48750000 77.000000 + 2.49250000 75.000000 + 2.49750000 90.000000 + 2.50250000 85.000000 + 2.50750000 72.000000 + 2.51250000 86.000000 + 2.51750000 85.000000 + 2.52250000 88.000000 + 2.52750000 87.000000 + 2.53250000 90.000000 + 2.53750000 94.000000 + 2.54250000 98.000000 + 2.54750000 95.000000 + 2.55250000 84.000000 + 2.55750000 70.000000 + 2.56250000 94.000000 + 2.56750000 80.000000 + 2.57250000 101.000000 + 2.57750000 96.000000 + 2.58250000 85.000000 + 2.58750000 86.000000 + 2.59250000 107.000000 + 2.59750000 92.000000 + 2.60250000 79.000000 + 2.60750000 93.000000 + 2.61250000 87.000000 + 2.61750000 77.000000 + 2.62250000 86.000000 + 2.62750000 85.000000 + 2.63250000 99.000000 + 2.63750000 80.000000 + 2.64250000 94.000000 + 2.64750000 98.000000 + 2.65250000 85.000000 + 2.65750000 112.000000 + 2.66250000 109.000000 + 2.66750000 87.000000 + 2.67250000 103.000000 + 2.67750000 119.000000 + 2.68250000 129.000000 + 2.68750000 109.000000 + 2.69250000 105.000000 + 2.69750000 109.000000 + 2.70250000 111.000000 + 2.70750000 108.000000 + 2.71250000 125.000000 + 2.71750000 127.000000 + 2.72250000 125.000000 + 2.72750000 119.000000 + 2.73250000 134.000000 + 2.73750000 156.000000 + 2.74250000 112.000000 + 2.74750000 133.000000 + 2.75250000 145.000000 + 2.75750000 137.000000 + 2.76250000 137.000000 + 2.76750000 153.000000 + 2.77250000 149.000000 + 2.77750000 136.000000 + 2.78250000 135.000000 + 2.78750000 135.000000 + 2.79250000 125.000000 + 2.79750000 162.000000 + 2.80250000 136.000000 + 2.80750000 115.000000 + 2.81250000 153.000000 + 2.81750000 131.000000 + 2.82250000 142.000000 + 2.82750000 156.000000 + 2.83250000 121.000000 + 2.83750000 115.000000 + 2.84250000 134.000000 + 2.84750000 139.000000 + 2.85250000 129.000000 + 2.85750000 137.000000 + 2.86250000 142.000000 + 2.86750000 133.000000 + 2.87250000 126.000000 + 2.87750000 125.000000 + 2.88250000 136.000000 + 2.88750000 149.000000 + 2.89250000 127.000000 + 2.89750000 151.000000 + 2.90250000 111.000000 + 2.90750000 121.000000 + 2.91250000 117.000000 + 2.91750000 122.000000 + 2.92250000 109.000000 + 2.92750000 118.000000 + 2.93250000 115.000000 + 2.93750000 106.000000 + 2.94250000 110.000000 + 2.94750000 124.000000 + 2.95250000 104.000000 + 2.95750000 114.000000 + 2.96250000 92.000000 + 2.96750000 95.000000 + 2.97250000 98.000000 + 2.97750000 95.000000 + 2.98250000 101.000000 + 2.98750000 91.000000 + 2.99250000 87.000000 + 2.99750000 94.000000 + 3.00250000 102.000000 + 3.00750000 110.000000 + 3.01250000 81.000000 + 3.01750000 83.000000 + 3.02250000 77.000000 + 3.02750000 87.000000 + 3.03250000 97.000000 + 3.03750000 82.000000 + 3.04250000 67.000000 + 3.04750000 82.000000 + 3.05250000 81.000000 + 3.05750000 82.000000 + 3.06250000 76.000000 + 3.06750000 78.000000 + 3.07250000 93.000000 + 3.07750000 80.000000 + 3.08250000 77.000000 + 3.08750000 80.000000 + 3.09250000 70.000000 + 3.09750000 85.000000 + 3.10250000 87.000000 + 3.10750000 77.000000 + 3.11250000 86.000000 + 3.11750000 75.000000 + 3.12250000 78.000000 + 3.12750000 91.000000 + 3.13250000 82.000000 + 3.13750000 87.000000 + 3.14250000 85.000000 + 3.14750000 74.000000 + 3.15250000 82.000000 + 3.15750000 77.000000 + 3.16250000 82.000000 + 3.16750000 69.000000 + 3.17250000 87.000000 + 3.17750000 98.000000 + 3.18250000 82.000000 + 3.18750000 73.000000 + 3.19250000 64.000000 + 3.19750000 91.000000 + 3.20250000 67.000000 + 3.20750000 77.000000 + 3.21250000 87.000000 + 3.21750000 70.000000 + 3.22250000 75.000000 + 3.22750000 74.000000 + 3.23250000 69.000000 + 3.23750000 65.000000 + 3.24250000 60.000000 + 3.24750000 77.000000 + 3.25250000 82.000000 + 3.25750000 74.000000 + 3.26250000 68.000000 + 3.26750000 90.000000 + 3.27250000 77.000000 + 3.27750000 80.000000 + 3.28250000 69.000000 + 3.28750000 76.000000 + 3.29250000 70.000000 + 3.29750000 85.000000 + 3.30250000 65.000000 + 3.30750000 79.000000 + 3.31250000 74.000000 + 3.31750000 77.000000 + 3.32250000 74.000000 + 3.32750000 75.000000 + 3.33250000 77.000000 + 3.33750000 89.000000 + 3.34250000 82.000000 + 3.34750000 91.000000 + 3.35250000 80.000000 + 3.35750000 80.000000 + 3.36250000 76.000000 + 3.36750000 83.000000 + 3.37250000 69.000000 + 3.37750000 67.000000 + 3.38250000 77.000000 + 3.38750000 80.000000 + 3.39250000 69.000000 + 3.39750000 87.000000 + 3.40250000 72.000000 + 3.40750000 76.000000 + 3.41250000 72.000000 + 3.41750000 67.000000 + 3.42250000 55.000000 + 3.42750000 78.000000 + 3.43250000 86.000000 + 3.43750000 73.000000 + 3.44250000 67.000000 + 3.44750000 72.000000 + 3.45250000 81.000000 + 3.45750000 68.000000 + 3.46250000 74.000000 + 3.46750000 81.000000 + 3.47250000 72.000000 + 3.47750000 82.000000 + 3.48250000 81.000000 + 3.48750000 82.000000 + 3.49250000 69.000000 + 3.49750000 72.000000