From cd5df8686bd067f545e8789b7b0ad53e0ba330d5 Mon Sep 17 00:00:00 2001
From: flavigny <flavigny@lpccaen.in2p3.fr>
Date: Tue, 22 Feb 2022 14:35:08 +0100
Subject: [PATCH] * Update 3-alpha source fit functions and procedure
---
NPLib/Calibration/NPCalibrationSource.cxx | 243 +++++++++++++++++++---
NPLib/Calibration/NPCalibrationSource.h | 9 +
NPLib/Calibration/NPSiliconCalibrator.cxx | 190 ++++++++++++++++-
NPLib/Calibration/NPSiliconCalibrator.h | 2 +-
4 files changed, 409 insertions(+), 35 deletions(-)
diff --git a/NPLib/Calibration/NPCalibrationSource.cxx b/NPLib/Calibration/NPCalibrationSource.cxx
index e0a342291..482747079 100644
--- a/NPLib/Calibration/NPCalibrationSource.cxx
+++ b/NPLib/Calibration/NPCalibrationSource.cxx
@@ -7,6 +7,7 @@
////////////////////////////////////////////////////////////////////////////////
NPL::CalibrationSource::CalibrationSource(){
m_SourceSignal = 0 ;
+m_FitFunctionType = 1 ;
}
////////////////////////////////////////////////////////////////////////////////
@@ -15,36 +16,224 @@ NPL::CalibrationSource::~CalibrationSource(){
////////////////////////////////////////////////////////////////////////////////
void NPL::CalibrationSource::ComputeSourceSignal(){
- TString contribution;
- unsigned int arg_number = 0;
-
- for(unsigned int i = 0 ; i < m_Energies.size() ; i++){
- for(unsigned int j = 0 ; j < m_Energies[i].size() ; j++){
- contribution +=
- Form("(%f)*[%i]*exp(-(x-[%i]*%f)*(x-[%i]*%f)/(2*[%i]*[%i]))",
- m_BranchingRatio[i][j]/m_BranchingRatio[i][0], arg_number,
- arg_number+1, m_Energies[i][j]/m_Energies[i][0],
- arg_number+1, m_Energies[i][j]/m_Energies[i][0],
- arg_number+2,arg_number+2);
-
- if(j!=m_Energies[i].size()-1) contribution+="+";
- }
-
- arg_number+=3;
- if(i!=m_Energies.size()-1) contribution+="+";
- }
-
- m_SourceSignal = new TF1("np_source_signal",contribution,0,1e9);
- m_SourceSignal->SetNpx(20000);
- arg_number = 0;
- for(unsigned int i = 0 ; i < m_Energies.size() ; i++){
- m_SourceSignal->SetParameter(arg_number++,1000);
- m_SourceSignal->SetParameter(arg_number++,5000+i*1000);
- m_SourceSignal->SetParameter(arg_number++,20);
- }
+
+
+ if (m_FitFunctionType == 0){ //3-alpha source golobal (8 gaussians)
+ m_SourceSignal = new TF1("AlphaSourceGauss",AlphaSourceGaus,0,1e9,23);
+ m_SourceSignal->SetNpx(20000);
+
+ unsigned int arg_number = 0;
+ m_SourceSignal->SetParameter(arg_number++,10);
+
+ for(unsigned int i = 0 ; i < m_Energies.size() ; i++){
+ m_SourceSignal->SetParameter(arg_number++, 1000); //main peak amp
+ m_SourceSignal->SetParameter(arg_number++, 5000+i*1000); //main peak mean
+ for(unsigned int j = 0 ; j < m_Energies[i].size() ; j++){
+ m_SourceSignal->SetParameter(arg_number++,m_BranchingRatio[i][j]/m_BranchingRatio[i][0]);
+ m_SourceSignal->SetParameter(arg_number++,m_Energies[i][j]/m_Energies[i][0]);
+ }
+ }
+
+ } else if(m_FitFunctionType==1){// triplet of alpha peaks (3 gaussian with tail)
+
+ //m_SourceSignal = new TF1("AlphaTriplet",AlphaTriplet,0,1e9,9);
+ m_SourceSignal = new TF1("AlphaTripletWithTail",AlphaTripletWithTwoTail,0,1e9,12);
+ m_SourceSignal->SetNpx(20000);
+
+ unsigned int arg_number = 0;
+
+ m_SourceSignal->SetParameter(arg_number++,1000);
+ m_SourceSignal->SetParameter(arg_number++,5000+1*1000);
+ m_SourceSignal->SetParameter(arg_number++,20);
+ m_SourceSignal->SetParameter(arg_number++,20);
+ m_SourceSignal->FixParameter(arg_number++,m_BranchingRatio[1][0]/m_BranchingRatio[1][0]);
+ m_SourceSignal->FixParameter(arg_number++,m_BranchingRatio[1][1]/m_BranchingRatio[1][0]);
+ m_SourceSignal->SetParameter(arg_number++,m_BranchingRatio[1][2]/m_BranchingRatio[1][0]);
+ m_SourceSignal->SetParameter(arg_number++,m_Energies[1][0]/m_Energies[1][0]);
+ m_SourceSignal->SetParameter(arg_number++,m_Energies[1][1]/m_Energies[1][0]);
+ m_SourceSignal->SetParameter(arg_number++,m_Energies[1][2]/m_Energies[1][0]);
+ m_SourceSignal->SetParLimits(9,0.95*m_Energies[1][2]/m_Energies[1][0],1.05*m_Energies[1][2]/m_Energies[1][0]);
+ m_SourceSignal->SetParameter(arg_number++,100);
+ //m_SourceSignal->SetParLimits(11,10,1000);
+ m_SourceSignal->FixParameter(arg_number++,0.4);
+ //m_SourceSignal->SetParLimits(11,0.35,0.45);
+ }
+ else if (m_FitFunctionType == 2){//old way, inline i definition, 8 gaussians.
+ TString contribution;
+ unsigned int arg_number = 0;
+
+ for(unsigned int i = 0 ; i < m_Energies.size() ; i++){
+ for(unsigned int j = 0 ; j < m_Energies[i].size() ; j++){
+ contribution +=
+ Form("(%f)*[%i]*exp(-(x-[%i]*%f)*(x-[%i]*%f)/(2*[%i]*[%i]))",
+ m_BranchingRatio[i][j]/m_BranchingRatio[i][0], arg_number,
+ arg_number+1, m_Energies[i][j]/m_Energies[i][0],
+ arg_number+1, m_Energies[i][j]/m_Energies[i][0],
+ arg_number+2,arg_number+2);
+ if(j!=m_Energies[i].size()-1) contribution+="+";
+ }
+
+ arg_number+=3;
+ if(i!=m_Energies.size()-1) contribution+="+";
+ }
+
+ m_SourceSignal = new TF1("np_source_signal",contribution,0,1e9);
+ m_SourceSignal->SetNpx(20000);
+ arg_number = 0;
+ for(unsigned int i = 0 ; i < m_Energies.size() ; i++){
+ m_SourceSignal->SetParameter(arg_number++,1000);
+ m_SourceSignal->SetParameter(arg_number++,5000+i*1000);
+ m_SourceSignal->SetParameter(arg_number++,20);
+ }
+ }
}
+////////////////////////////////////////////////////////////////////////////////
+double AlphaSourceGaus(double *x,double *p)
+{
+ // eight gaussians: (3 for 239Pu, 3 for 241Am, 2 for Cm244)
+ // same sigma for all gaussians
+
+ //p0 = sigma of all peaks
+
+ //239 Pu
+ //p1 = amplitude of most intense peak
+ //p2 = mean of most intense peak
+ //p3 = relative scaling in amplitude
+ //p4 = relative scaling in mean for main peak
+ //p5 = relative scaling in amplitude for second peak
+ //p6 = relative scaling in mean for second peak
+ //p7 = relative scaling in amplitude for third peak
+ //p8 = relative scaling in mean for third peak
+ //241 Am
+ //p9 = amplitude of most intense peak
+ //p10 = mean of most intense peak
+ //p2 = sigma of most intense peak
+ //p11 = relative scaling in amplitude
+ //p12 = relative scaling in mean for main peak
+ //p13 = relative scaling in amplitude for second peak
+ //p14 = relative scaling in mean for second peak
+ //p15 = relative scaling in amplitude for third peak
+ //p16 = relative scaling in mean for third peak
+ //244 Cm
+ //p17 = amplitude of most intense peak
+ //p18 = mean of most intense peak
+ //p2 = sigma of most intense peak
+ //p19 = relative scaling in amplitude
+ //p20 = relative scaling in mean for main peak
+ //p21 = relative scaling in amplitude for second peak
+ //p22 = relative scaling in mean for second peak
+ double arg = 0;
+ double fitval = 0;
+ int ngaus = 8;
+ int k = 0;
+ int rest = 0;
+
+ if (p[2]==0) {cout<<"sigma=0 in AlphaTripletGaus"<<endl; return fitval;}
+
+ for(int i=0; i<ngaus; i++){
+ rest = i%3;
+ if (i!=0 && rest==0) k++;
+ arg = (x[0] - p[2+k*ngaus]*p[4+k*ngaus+2*rest])/p[0];
+ fitval += p[1+k*ngaus]*p[3+k*ngaus+2*rest]*TMath::Exp(-0.5*arg*arg);
+ }
+
+ return fitval;
+}
+////////////////////////////////////////////////////////////////////////////////
+double AlphaTripletGaus(double *x,double *p)
+{
+ // eight gaussians, same sigma
+ //p0 = amplitude of most intense peak
+ //p1 = mean of most intense peak
+ //p2 = sigma most intense peak
+ //p3 = relative scaling in amplitude
+ //p4 = relative scaling in mean for main peak
+ //p5 = relative scaling in amplitude for second peak
+ //p6 = relative scaling in mean for second peak
+ //p7 = relative scaling in amplitude for third peak
+ //p8 = relative scaling in mean for third peak
+ double arg = 0;
+ if (p[2]!=0) arg = (x[0] - p[1]*p[4])/p[2];
+ double fitval = p[0]*p[3]*TMath::Exp(-0.5*arg*arg);
+ if (p[2]!=0) arg = (x[0] - p[1]*p[6])/p[2];
+ fitval += p[0]*p[5]*TMath::Exp(-0.5*arg*arg);
+ if (p[2]!=0) arg = (x[0] - p[1]*p[8])/p[2];
+ fitval += p[0]*p[7]*TMath::Exp(-0.5*arg*arg);
+ return fitval;
+}
+////////////////////////////////////////////////////////////////////////////////
+double AlphaTripletWithTail(double *x,double *p)
+{
+ //p0 = area
+ //p1 = mean
+ //p2 = sigma
+ //p3 = tail
+ //p4 = amplitude scaling first peak
+ //p5 = amplitude scaling second peak
+ //p6 = amplitude scaling third peak
+ //p7 = energy scaling first peak
+ //p8 = energy scaling second peak
+ //p9 = energy scaling third peak
+ double arg1 = 0, arg2 = 0, arg3=0, arg4=0, arg5=0, arg6=0;
+ if (p[2]!=0 && p[3]!=0) {
+ arg1 = (x[0] - p[1]*p[7])/p[3] + p[2]*p[2]/(2*p[3]*p[3]);
+ arg2 = ( (x[0] - p[1]*p[7])/p[2] + p[2]/p[3] )/TMath::Sqrt(2);
+
+ arg3 = (x[0] - p[1]*p[8])/p[3] + p[2]*p[2]/(2*p[3]*p[3]);
+ arg4 = ( (x[0] - p[1]*p[8])/p[2] + p[2]/p[3] )/TMath::Sqrt(2);
+
+ arg5 = (x[0] - p[1]*p[9])/p[3] + p[2]*p[2]/(2*p[3]*p[3]);
+ arg6 = ( (x[0] - p[1]*p[9])/p[2] + p[2]/p[3] )/TMath::Sqrt(2);
+ }
+ double fitval = p[0]/(2*p[3])*p[4]*TMath::Exp(arg1)*TMath::Erfc(arg2)
+ + p[0]/(2*p[3])*p[5]*TMath::Exp(arg3)*TMath::Erfc(arg4)
+ + p[0]/(2*p[3])*p[6]*TMath::Exp(arg5)*TMath::Erfc(arg6);
+ return fitval;
+}
+////////////////////////////////////////////////////////////////////////////////
+double AlphaTripletWithTwoTail(double *x,double *p)
+{
+ //p0 = area
+ //p1 = mean
+ //p2 = sigma
+ //p3 = tail
+ //p4 = amplitude scaling first peak
+ //p5 = amplitude scaling second peak
+ //p6 = amplitude scaling third peak
+ //p7 = energy scaling first peak
+ //p8 = energy scaling second peak
+ //p9 = energy scaling third peak
+ //p10 = tail2
+ //p11 = eta
+ double arg1 = 0, arg2 = 0, arg3=0, arg4=0, arg5=0, arg6=0;
+ double arg1b = 0, arg2b = 0, arg3b=0, arg4b=0, arg5b=0, arg6b=0;
+ if (p[2]!=0 && p[3]!=0) {
+ arg1 = (x[0] - p[1]*p[7])/p[3] + p[2]*p[2]/(2*p[3]*p[3]);
+ arg2 = ( (x[0] - p[1]*p[7])/p[2] + p[2]/p[3] )/TMath::Sqrt(2);
+ arg1b = (x[0] - p[1]*p[7])/p[10] + p[2]*p[2]/(2*p[10]*p[10]);
+ arg2b = ( (x[0] - p[1]*p[7])/p[2] + p[2]/p[10] )/TMath::Sqrt(2);
+
+ arg3 = (x[0] - p[1]*p[8])/p[3] + p[2]*p[2]/(2*p[3]*p[3]);
+ arg4 = ( (x[0] - p[1]*p[8])/p[2] + p[2]/p[3] )/TMath::Sqrt(2);
+ arg3b = (x[0] - p[1]*p[8])/p[10] + p[2]*p[2]/(2*p[10]*p[10]);
+ arg4b = ( (x[0] - p[1]*p[8])/p[2] + p[2]/p[10] )/TMath::Sqrt(2);
+
+ arg5 = (x[0] - p[1]*p[9])/p[3] + p[2]*p[2]/(2*p[3]*p[3]);
+ arg6 = ( (x[0] - p[1]*p[9])/p[2] + p[2]/p[3] )/TMath::Sqrt(2);
+ arg5b = (x[0] - p[1]*p[9])/p[10] + p[2]*p[2]/(2*p[10]*p[10]);
+ arg6b = ( (x[0] - p[1]*p[9])/p[2] + p[2]/p[10] )/TMath::Sqrt(2);
+ }
+ double fitval =
+ p[0]/2*p[4]*((1-p[11])/p[3]*TMath::Exp(arg1)*TMath::Erfc(arg2)
+ + p[11]/p[10]*TMath::Exp(arg1b)*TMath::Erfc(arg2b))
+ + p[0]/2*p[5]*((1-p[11])/p[3]*TMath::Exp(arg3)*TMath::Erfc(arg4)
+ + p[11]/p[10]*TMath::Exp(arg3b)*TMath::Erfc(arg4b))
+ + p[0]/2*p[6]*((1-p[11])/p[3]*TMath::Exp(arg5)*TMath::Erfc(arg6)
+ + p[11]/p[10]*TMath::Exp(arg5b)*TMath::Erfc(arg6b));
+ return fitval;
+}
////////////////////////////////////////////////////////////////////////////////
TF1* NPL::CalibrationSource::GetSourceSignal(){
if(!m_SourceSignal)
diff --git a/NPLib/Calibration/NPCalibrationSource.h b/NPLib/Calibration/NPCalibrationSource.h
index b7919854e..9bca6fa24 100644
--- a/NPLib/Calibration/NPCalibrationSource.h
+++ b/NPLib/Calibration/NPCalibrationSource.h
@@ -14,6 +14,11 @@ using namespace std;
// NPL
#include "NPEnergyLoss.h"
+double AlphaSourceGaus(double *,double *);
+double AlphaTripletGaus(double *,double *);
+double AlphaTripletWithTail(double *,double *);
+double AlphaTripletWithTwoTail(double *,double *);
+
namespace NPL{
class CalibrationSource{
@@ -26,6 +31,7 @@ namespace NPL{
vector< vector<double> > m_Energies;
vector< vector<double> > m_BranchingRatio;
vector< vector<double> > m_ErrEnergies;
+ int m_FitFunctionType;
private: // Computed source signal (a TF1 ready for fit use, set up with physical value)
TF1* m_SourceSignal;
@@ -33,6 +39,8 @@ namespace NPL{
public: // Getter and Setter
TF1* GetSourceSignal();
+ void SetFitFunctionType(int type){m_FitFunctionType=type;};
+ int GetFitFunctionType(){return m_FitFunctionType;};
public: // Get and Set the energies
// for particle like source (one gaussian per energies)
@@ -50,5 +58,6 @@ namespace NPL{
void Set_244Cm();
};
+
}
#endif
diff --git a/NPLib/Calibration/NPSiliconCalibrator.cxx b/NPLib/Calibration/NPSiliconCalibrator.cxx
index 2cfb5961a..3658708db 100644
--- a/NPLib/Calibration/NPSiliconCalibrator.cxx
+++ b/NPLib/Calibration/NPSiliconCalibrator.cxx
@@ -24,7 +24,7 @@ NPL::SiliconCalibrator::~SiliconCalibrator(){
}
//////////////////////////////////////////////
-double NPL::SiliconCalibrator::ZeroExtrapolation(TH1* histo, NPL::CalibrationSource* CS, NPL::EnergyLoss* EL, vector<double>& coeff, unsigned int pedestal, unsigned int max_iteration, double rmin,double rmax){
+double NPL::SiliconCalibrator::ZeroExtrapolation(TH1* histo, NPL::CalibrationSource* CS, NPL::EnergyLoss* EL, vector<double>& coeff, double pedestal, unsigned int max_iteration, double rmin,double rmax){
if(histo->GetEntries()==0){
coeff.clear();
coeff.push_back(0);
@@ -63,7 +63,8 @@ double NPL::SiliconCalibrator::ZeroExtrapolation(TH1* histo, NPL::CalibrationSou
bool LastStepDecrease = false;
// 0.1% precision on the total scale (typically 6 MeV )
- double my_precision = (rmax-pedestal)*0.001;
+ //double my_precision = (rmax-pedestal)*0.001;
+ double my_precision = 1;
// Energy of the Source and sigma on the value (MeV)
double Assume_Thickness = 0 ; // micrometer
vector<double> Assume_E; // Energie calculated assuming Assume_Thickness deadlayer of Al
@@ -96,15 +97,17 @@ double NPL::SiliconCalibrator::ZeroExtrapolation(TH1* histo, NPL::CalibrationSou
g->SetPoint(i,x[i] ,Assume_E[i]);
}
- DistanceToPedestal = FitPoints(g,&Assume_E[0], Source_Sig, coeff, 0 );
-
+ DistanceToPedestal = FitPoints(g,&Assume_E[0], Source_Sig, coeff, pedestal );
+/*
+ cout << "test " << histo->GetName() << ", Assume_Thick =" << Assume_Thickness / micrometer << "um with step "<< Step/micrometer << ", Dst to Pedestal = " << DistanceToPedestal <<", pedestal= "<<pedestal << endl;
+*/
if(abs(DistanceToPedestal) < my_precision || Step < 0.01*micrometer){
break;
}
else if(DistanceToPedestal < 0 ){
if(LastStepIncrease)
- Step *= 5; //multiply step by 5
+ //Step *= 5; //multiply step by 5
Assume_Thickness += Step;
LastStepIncrease = true;
@@ -112,7 +115,7 @@ double NPL::SiliconCalibrator::ZeroExtrapolation(TH1* histo, NPL::CalibrationSou
else if(DistanceToPedestal > 0 ){
if(LastStepDecrease)
- Step *= 0.5; //decrease step by half
+ //Step *= 0.5; //decrease step by half
Assume_Thickness -= Step;
LastStepDecrease = true;
@@ -249,7 +252,179 @@ TGraphErrors* NPL::SiliconCalibrator::FitSpectrum(TH1* histo, double rmin, doubl
// sigma: must be >=1, higher sgma higher resolution
// option: nobackground => Don't subtract back ground
// threshold: look for peaks within range [HighespeakAmp*threshold ; HighespeakAmp]
- Int_t nfound = sp.Search(histo,2,"",0.15);
+ Int_t nfound = sp.Search(histo,2,"",0.6);
+ //std::cout << "Number of peaks found = " << nfound << std::endl; // used for debugging
+ TSpectrumPosition_t xpeaks = sp.GetPositionX();
+
+ // order list of peaks
+ sort(xpeaks, xpeaks+nfound);
+
+ m_FitFunction = m_CalibrationSource->GetSourceSignal();
+ int FitFunctionType = m_CalibrationSource->GetFitFunctionType();
+ if (nfound == 3) {
+
+ if(FitFunctionType==0){
+ // 3 * (triple gaussian (no tail))
+ vector< vector<double> > Energies = m_CalibrationSource->GetEnergies();
+ vector< vector<double> > ErrEnergies = m_CalibrationSource->GetEnergiesErrors();
+ vector< vector<double> > BranchingRatio = m_CalibrationSource->GetBranchingRatio();
+
+ // ballpark the sigma
+ double sigma = (xpeaks[1]-xpeaks[0])/10;
+ m_FitFunction->FixParameter(0,sigma);
+ m_FitFunction->SetParLimits(0,0.1,sigma*3);
+
+ int arg_number=1;
+ double H=0, Bratio=0, Eratio=0;
+ for(unsigned int i = 0 ; i < Energies.size() ; i++){
+ H = histo->GetBinContent(histo->FindBin(xpeaks[i]));
+ m_FitFunction->SetParameter(arg_number,H ); //main peak amp
+ m_FitFunction->SetParLimits(arg_number, H*0.8, H*1.2 );
+ arg_number++;
+ m_FitFunction->FixParameter(arg_number, xpeaks[i]); //main peak mean
+ //m_FitFunction->SetParLimits(arg_number, xpeaks[i]*0.99, xpeaks[i]*1.01);
+ arg_number++;
+ for(unsigned int j = 0 ; j < Energies[i].size() ; j++){
+ Bratio = BranchingRatio[i][j]/BranchingRatio[i][0];
+ Eratio = Energies[i][j]/Energies[i][0];
+ m_FitFunction->FixParameter(arg_number,Bratio);
+ //m_FitFunction->SetParLimits(arg_number,Bratio*0.95,Bratio*1.05);
+ arg_number++;
+ m_FitFunction->SetParameter(arg_number,Eratio);
+ m_FitFunction->SetParLimits(arg_number,Eratio*0.998,Eratio*1.01);
+ arg_number++;
+ }
+ }
+
+ // Set range and fit spectrum
+ histo->GetXaxis()->SetRangeUser(xpeaks[0]-xpeaks[0]/20.,xpeaks[2]+xpeaks[2]/20.);
+ histo->Fit(m_FitFunction,"MQR");
+ TF1* fit = histo->GetFunction(m_FitFunction->GetName());
+ std::cout << "sigma: " << fit->GetParameter(0)<< std::endl;
+
+ // Graph energies(MeV) VS fitted peak channel and associated uncertainties
+ TGraphErrors* gre = new TGraphErrors();
+ int point = 0 ;
+ for (unsigned int i = 0; i < Energies.size(); i++) {
+ gre->SetPoint(point, fit->GetParameter(2+8*i), Energies[i][0]);
+ gre->SetPointError(point++, fit->GetParError(2+8*i), ErrEnergies[i][0]);
+ }
+ //std::cout << "gre->GetN() = " << gre->GetN() << std::endl; // used for debugging
+ return gre;
+
+
+ } else if(FitFunctionType==1){// 3 * (triple gaussian with tail)
+
+ //WORK IN PROGRESS (only one triplet fit)
+ double H2 = histo->GetBinContent(histo->FindBin(xpeaks[1]));
+
+ m_FitFunction->SetParameter(0,H2*30);
+ m_FitFunction->SetParameter(1,xpeaks[1]);
+ double sigma = (xpeaks[1]-xpeaks[0])/50;
+ m_FitFunction->SetParameter(2,sigma);
+ m_FitFunction->SetParameter(3,sigma); //tail=sigma as a start
+ m_FitFunction->SetParameter(10,5*sigma); //tail2=5*sigma as a start
+ m_FitFunction->SetParLimits(10,2*sigma,10*sigma); //tail2=5*sigma as a start
+ // Set range and fit spectrum
+ //histo->GetXaxis()->SetRangeUser(xpeaks[0]-xpeaks[0]/20.,xpeaks[2]+xpeaks[2]/20.);
+ histo->GetXaxis()->SetRangeUser(xpeaks[1]-2000.,xpeaks[1]+1000);
+ histo->Fit(m_FitFunction,"MQR");
+ TF1* fit = histo->GetFunction(m_FitFunction->GetName());
+
+ vector< vector<double> > Energies = m_CalibrationSource->GetEnergies();
+ vector< vector<double> > ErrEnergies = m_CalibrationSource->GetEnergiesErrors();
+
+ TGraphErrors* gre = new TGraphErrors();
+ int point = 0 ;
+ for (unsigned int i = 0; i < Energies.size(); i++) {
+ gre->SetPoint(point, fit->GetParameter(3*i+1), Energies[i][0]);
+ gre->SetPointError(point++, fit->GetParError(3*i+1), ErrEnergies[i][0]);
+ }
+ return gre;
+
+ }else if(FitFunctionType==2){//old way
+
+ // Set initial mean
+ m_FitFunction->SetParameter(1,xpeaks[0]);
+ m_FitFunction->SetParameter(4,xpeaks[1]);
+ m_FitFunction->SetParameter(7,xpeaks[2]);
+
+ // Set Limit
+ m_FitFunction->SetParLimits(1,xpeaks[0]*0.8,xpeaks[0]*1.2);
+ m_FitFunction->SetParLimits(4,xpeaks[1]*0.8,xpeaks[1]*1.2);
+ m_FitFunction->SetParLimits(7,xpeaks[2]*0.8,xpeaks[2]*1.2);
+
+
+ // Set initial height
+ double H1 = histo->GetBinContent(histo->FindBin(xpeaks[0]));
+ double H2 = histo->GetBinContent(histo->FindBin(xpeaks[1]));
+ double H3 = histo->GetBinContent(histo->FindBin(xpeaks[2]));
+
+ m_FitFunction->SetParameter(0,H1);
+ m_FitFunction->SetParameter(3,H2);
+ m_FitFunction->SetParameter(6,H3);
+
+ // Set Limit
+ m_FitFunction->SetParLimits(0,H1*0.4,H1*2);
+ m_FitFunction->SetParLimits(3,H2*0.4,H2*2);
+ m_FitFunction->SetParLimits(6,H3*0.4,H3*2);
+
+ // ballpark the sigma
+ double sigma = (xpeaks[1]-xpeaks[0])/10;
+ m_FitFunction->SetParameter(2,sigma);
+ m_FitFunction->SetParameter(5,sigma);
+ m_FitFunction->SetParameter(8,sigma);
+
+ // Set Limit
+ m_FitFunction->SetParLimits(2,0,sigma*3);
+ m_FitFunction->SetParLimits(5,0,sigma*3);
+ m_FitFunction->SetParLimits(8,0,sigma*3);
+
+ // Set range and fit spectrum
+ histo->GetXaxis()->SetRangeUser(xpeaks[0]-xpeaks[0]/20.,xpeaks[2]+xpeaks[2]/20.);
+ histo->Fit(m_FitFunction,"MQR");
+ TF1* fit = histo->GetFunction(m_FitFunction->GetName());
+ // energies and associated uncertainties
+ vector< vector<double> > Energies = m_CalibrationSource->GetEnergies();
+ vector< vector<double> > ErrEnergies = m_CalibrationSource->GetEnergiesErrors();
+
+ // unsigned int mysize = Energies.size();
+ TGraphErrors* gre = new TGraphErrors();
+ int point = 0 ;
+ for (unsigned int i = 0; i < Energies.size(); i++) {
+ gre->SetPoint(point, fit->GetParameter(3*i+1), Energies[i][0]);
+ gre->SetPointError(point++, fit->GetParError(3*i+1), ErrEnergies[i][0]);
+ }
+ //std::cout << "gre->GetN() = " << gre->GetN() << std::endl; // used for debugging
+ return gre;
+
+ } else{ //Fit functiuon type not 0,1 or 2 (undefined)
+ return (new TGraphErrors());
+ }
+
+ }
+
+ else{ //nfound!=3
+ return (new TGraphErrors());
+ }
+}
+/*
+TGraphErrors* NPL::SiliconCalibrator::FitSpectrum(TH1* histo, double rmin, double rmax)
+{
+// for(unsigned int i = 0 ; i < histo->GetNbinsX() ; i++)
+// if(histo->GetBinCenter(i)<rmin || histo->GetBinCenter(i)>rmax )
+// histo->SetBinContent(i,0);
+
+ // apply range for peak search
+ histo->GetXaxis()->SetRangeUser(rmin, rmax);
+ // Perform a peak search to get a hint of where are the peaks
+ TSpectrum sp(4,1); // number of peaks is 4 to avoid warning of peak buffer
+ //arguments:
+ // hist : input histogram
+ // sigma: must be >=1, higher sgma higher resolution
+ // option: nobackground => Don't subtract back ground
+ // threshold: look for peaks within range [HighespeakAmp*threshold ; HighespeakAmp]
+ Int_t nfound = sp.Search(histo,2,"",0.6);
//std::cout << "Number of peaks found = " << nfound << std::endl; // used for debugging
TSpectrumPosition_t xpeaks = sp.GetPositionX();
@@ -317,5 +492,6 @@ TGraphErrors* NPL::SiliconCalibrator::FitSpectrum(TH1* histo, double rmin, doubl
return (new TGraphErrors());
}
}
+*/
diff --git a/NPLib/Calibration/NPSiliconCalibrator.h b/NPLib/Calibration/NPSiliconCalibrator.h
index eab6a935d..1b359a573 100644
--- a/NPLib/Calibration/NPSiliconCalibrator.h
+++ b/NPLib/Calibration/NPSiliconCalibrator.h
@@ -24,7 +24,7 @@ namespace NPL{
public:
// Use the Zero Extrapolation method to perform fit and return the dead layer thickness
- double ZeroExtrapolation(TH1* histo, NPL::CalibrationSource* CS, NPL::EnergyLoss* EL, vector<double>& coeff, unsigned int pedestal, unsigned int max_iteration = 10000 , double rmin=-1,double rmax=-1);
+ double ZeroExtrapolation(TH1* histo, NPL::CalibrationSource* CS, NPL::EnergyLoss* EL, vector<double>& coeff, double pedestal, unsigned int max_iteration = 10000 , double rmin=-1,double rmax=-1);
double SimpleCalibration(TH1* histo, NPL::CalibrationSource* CS, NPL::EnergyLoss* EL, vector<double>& coeff, double Assume_Thickness=0,double rmin=-1,double rmax=-1);
--
GitLab