2788 lines
85 KiB
C
2788 lines
85 KiB
C
/************************************************************
|
|
#
|
|
# Created by Tsz Leung (Ryan) Tang
|
|
# goluckyryan@gmail.com
|
|
# rtang@fsu.edu
|
|
#
|
|
*************************************************************/
|
|
|
|
#ifndef AutoFit_C
|
|
#define AutoFit_C
|
|
|
|
#include <TF1.h>
|
|
#include <TGraph.h>
|
|
#include <TColor.h>
|
|
#include <TSpectrum.h>
|
|
#include <TMath.h>
|
|
#include <vector>
|
|
|
|
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(" fitDecay() - fit A * Exp(-t/T) + B\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(" saveFitPara() - Save the initial guess parameters.\n");
|
|
printf("---------------------------------------------------------\n");
|
|
}
|
|
|
|
std::vector<std::string> SplitStrAF(std::string tempLine, std::string splitter, int shift = 0){
|
|
|
|
std::vector<std::string> 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 = max(0., (1+2*cos(ang))/3.);
|
|
double g = max(0., (1 - cos(ang) + sqrt(3)* sin(ang))/3.);
|
|
double b = 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;p<nPeaks;p++) {
|
|
Double_t norm = par[3*p+0];
|
|
Double_t mean = par[3*p+1];
|
|
Double_t sigma = par[3*p+2];
|
|
result += norm * TMath::Gaus(x[0],mean,sigma, 1); // normalized Gaussian
|
|
}
|
|
return result;
|
|
}
|
|
|
|
int nPols = 1;
|
|
|
|
Double_t nPolFunc(Double_t *x, Double_t *par) {
|
|
Double_t result = 0;
|
|
|
|
for( int p = 0; p < nPols+1; p++){
|
|
result += par[p]*TMath::Power(x[0], p);
|
|
}
|
|
return result;
|
|
}
|
|
|
|
Double_t nGaussPol(Double_t *x, Double_t *par) {
|
|
|
|
Double_t result = 0;
|
|
for (Int_t p=0;p<nPeaks;p++) {
|
|
Double_t norm = par[3*p+0];
|
|
Double_t mean = par[3*p+1];
|
|
Double_t sigma = par[3*p+2];
|
|
result += norm * TMath::Gaus(x[0],mean,sigma, 1); // normalized Gaussian
|
|
}
|
|
|
|
for( int p = 0; p < nPols+1; p++){
|
|
result += par[3*nPeaks + p]*TMath::Power(x[0], p);
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
Double_t nGF3(Double_t *x, Double_t *par){
|
|
|
|
/** this is the fitting function for gamma peak from gf3, RadWare **/
|
|
Double_t result = 0;
|
|
for( Int_t p = 0; p < nPeaks; p++){
|
|
|
|
Double_t norm = par[6*p+0];
|
|
Double_t mean = par[6*p+1];
|
|
Double_t sigma = par[6*p+2];
|
|
Double_t ratio = par[6*p+3];
|
|
Double_t beta = par[6*p+4]; // skewness
|
|
Double_t step = par[6*p+5];
|
|
|
|
result += norm * (1-ratio)* TMath::Gaus(x[0], mean, sigma, 1) ;
|
|
result += norm * ratio * exp( sigma * sigma/2/beta/beta) / (2* beta )* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ;
|
|
result += norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) );
|
|
}
|
|
return result;
|
|
}
|
|
|
|
|
|
Double_t nGF3Pol(Double_t *x, Double_t *par){
|
|
|
|
/** this is the fitting function for gamma peak from gf3, RadWare **/
|
|
Double_t result = 0;
|
|
for( Int_t p = 0; p < nPeaks; p++){
|
|
Double_t norm = par[6*p+0];
|
|
Double_t mean = par[6*p+1];
|
|
Double_t sigma = par[6*p+2];
|
|
Double_t ratio = par[6*p+3];
|
|
Double_t beta = par[6*p+4]; // skewness
|
|
Double_t step = par[6*p+5];
|
|
|
|
result += norm * (1.0-ratio)* TMath::Gaus(x[0], mean, sigma, 1) ;
|
|
result += norm * ratio * exp( sigma * sigma/2/beta/beta) / (2* beta )* exp((x[0]-mean)/beta) * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) + sigma/beta/sqrt(2)) ;
|
|
result += norm * step * TMath::Erfc( (x[0]-mean)/(sigma * sqrt(2)) );
|
|
}
|
|
|
|
for( int p = 0; p < nPols + 1; p++){
|
|
result += par[6*nPeaks + p]*TMath::Power(x[0], p);
|
|
}
|
|
|
|
return result;
|
|
}
|
|
|
|
void PrintPar(TF1 * fit, int numParPerPeak){
|
|
int totPar = fit->GetNpar();
|
|
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");
|
|
}
|
|
|
|
vector<double> energy, height, sigma, lowE, highE ;
|
|
vector<int> 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());
|
|
ifstream file;
|
|
file.open(fitParaFile.Data());
|
|
|
|
if( !file){
|
|
printf("\ncannot read file : %s \n", fitParaFile.Data());
|
|
return 0;
|
|
}
|
|
|
|
while( file.good()) {
|
|
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());
|
|
|
|
vector<string> 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);
|
|
double xMin = hRes->GetXaxis()->GetXmin();
|
|
double xMax = hRes->GetXaxis()->GetXmax();
|
|
fit->SetRange(xMin, xMax);
|
|
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");
|
|
|
|
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);
|
|
|
|
}
|
|
|
|
//########################################
|
|
//#### fit 2 gauss + pol-1 // not updated
|
|
//########################################
|
|
vector<double> 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");
|
|
|
|
|
|
vector<double> 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);
|
|
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");
|
|
|
|
|
|
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
|
|
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
|
|
//##############################################
|
|
vector<double> 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");
|
|
|
|
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 );
|
|
vector<double> 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);
|
|
|
|
vector<double> 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"));
|
|
|
|
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]));
|
|
}
|
|
|
|
return exPos;
|
|
|
|
}
|
|
|
|
|
|
//########################################
|
|
//###### NOT tested
|
|
//########################################
|
|
vector<double> 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();
|
|
}
|
|
|
|
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 );
|
|
vector<double> 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);
|
|
|
|
vector<double> 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;
|
|
|
|
}
|
|
|
|
//########################################
|
|
//###### fir N Gauss with estimated BG
|
|
//########################################
|
|
void fitNGauss(TH1F * hist, int bgEst = 10, 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");
|
|
|
|
|
|
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"));
|
|
|
|
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]));
|
|
}
|
|
|
|
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");
|
|
|
|
|
|
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"));
|
|
|
|
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]));
|
|
}
|
|
|
|
}
|
|
|
|
|
|
//########################################
|
|
//#### 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");
|
|
|
|
|
|
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"));
|
|
|
|
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]));
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
|
|
//########################################
|
|
//#### 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");
|
|
|
|
|
|
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), paraA[3*nPeaks+i], paraE[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"));
|
|
|
|
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]));
|
|
}
|
|
|
|
}
|
|
|
|
//########################################
|
|
//###### fit A * Exp(-t/T) + B
|
|
//########################################
|
|
void fitDecay(TH1F * hist){
|
|
|
|
gStyle->SetOptStat("");
|
|
TCanvas * cFitDecay = NULL;
|
|
if( gROOT->FindObjectAny("cFitDecay") == NULL ){
|
|
cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000);
|
|
}else{
|
|
delete gROOT->FindObjectAny("cFitDecay") ;
|
|
cFitDecay = new TCanvas("cFitDecay", Form("fit A * Exp(-t/T) + B | fitDecay"), 800, 1000);
|
|
}
|
|
|
|
cFitDecay->Divide(1, 2);
|
|
|
|
double yMax = hist->GetBinContent(hist->GetMaximumBin());
|
|
hist->SetMaximum(yMax*1.2);
|
|
cFitDecay->cd(1); hist->Draw();
|
|
double estHalf = hist->GetBinCenter(hist->FindLastBinAbove(yMax/2.)) *1.3;
|
|
|
|
//printf("estimator : A = %.2f, T = %.2f \n", yMax, estHalf);
|
|
|
|
TF1 * fit = new TF1("fit", "[0]*TMath::Exp(-x/[1]/TMath::Log(2))+[2]");
|
|
fit->SetLineWidth(3);
|
|
fit->SetLineColor(1);
|
|
fit->SetNpx(1000);
|
|
|
|
fit->SetParameter(0, yMax );
|
|
fit->SetParameter(1, estHalf );
|
|
fit->SetParameter(2, 1);
|
|
|
|
fit->SetParLimits(0, 0, 1e9);
|
|
fit->SetParLimits(1, 0, 1e9);
|
|
fit->SetParLimits(2, 0, 1e9);
|
|
|
|
gStyle->SetOptFit(000000);
|
|
hist->Fit("fit", "q");
|
|
|
|
const Double_t* paraE = fit->GetParErrors();
|
|
const Double_t* paraA = fit->GetParameters();
|
|
|
|
TLatex text;
|
|
text.SetNDC();
|
|
text.SetTextFont(82);
|
|
text.SetTextSize(0.04);
|
|
text.SetTextColor(1);
|
|
|
|
text.DrawLatex(0.6, 0.8, "A e^{#left(-#frac{t}{T Ln(2)}#right)}+B");
|
|
text.DrawLatex(0.6, 0.75, Form("A : %.f(%.f)", paraA[0], paraE[0]));
|
|
text.DrawLatex(0.6, 0.7, Form("T : %.f(%.f)", paraA[1], paraE[1]));
|
|
text.DrawLatex(0.6, 0.65, Form("B : %.f(%.f)", paraA[2], paraE[2]));
|
|
|
|
double chi2 = fit->GetChisquare();
|
|
int ndf = fit->GetNDF();
|
|
text.DrawLatex(0.6, 0.6, Form("#bar{#chi^{2}} : %5.3f", chi2/ndf));
|
|
|
|
|
|
TF1 * fite = new TF1("fite", "[0]*TMath::Exp(-x/[1])" );
|
|
fite->SetParameter(0, paraA[0]);
|
|
fite->SetParameter(1, paraA[1]);
|
|
fite->Draw("same");
|
|
|
|
TF1 * fitbg = new TF1("fitbg", "pol0" );
|
|
fitbg->SetParameter(0, paraA[2]);
|
|
fitbg->Draw("same");
|
|
|
|
cFitDecay->cd(2);
|
|
PlotResidual(hist, fit);
|
|
|
|
|
|
}
|
|
|
|
//#################################################
|
|
//#### fit N-Gauss with pol-n BG using mouse click
|
|
//#################################################
|
|
|
|
int nClick = 0;
|
|
bool peakFlag = 1;
|
|
vector<double> xPeakList;
|
|
vector<double> yPeakList;
|
|
vector<double> xBGList;
|
|
vector<double> 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(TString fileName = "AutoFit_para.txt"){
|
|
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");
|
|
|
|
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){
|
|
|
|
printf("=========================================================\n");
|
|
printf("======== fit n-Gauss + Pol-%d BG using mouse click =====\n", degPol );
|
|
printf("==========================================================\n");
|
|
|
|
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);
|
|
|
|
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 );
|
|
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"));
|
|
|
|
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]));
|
|
}
|
|
}
|
|
|
|
|
|
|
|
void clickFitNGaussPolSub(TH1F * hist, int degPol, double sigmaMax = 0){
|
|
|
|
printf("=========================================================\n");
|
|
printf("= fit n-Gauss + Pol-%d BG using mouse click (method-2) =\n", degPol );
|
|
printf("==========================================================\n");
|
|
|
|
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();
|
|
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("--------------- 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);
|
|
|
|
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 );
|
|
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"));
|
|
|
|
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]));
|
|
}
|
|
}
|
|
|
|
|
|
//########################################
|
|
//#### 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
|