Commit 4177d91b authored by Jérémie Dudouet's avatar Jérémie Dudouet
Browse files

Add Generic polynomial calibration

parent ea53e003
Pipeline #54525 passed with stage
in 8 minutes and 54 seconds
......@@ -72,7 +72,7 @@ ROOT_GENERATE_DICTIONARY(G__${Lib_NAME} ${headers} MODULE ${Lib_NAME} LINKDEF ${
#####################
#---Create a shared library with geneated dictionary
add_library(${Lib_NAME} SHARED ${sources} G__${Lib_NAME}.cxx)
add_library(${Lib_NAME} SHARED ${sources} ${headers} G__${Lib_NAME}.cxx)
target_link_libraries(${Lib_NAME} ${EXTRA_EXTERNAL_LIBRARIES} ${EXTRA_INTERNAL_LIBRARIES})
#####################
......
......@@ -39,14 +39,13 @@ CXHist1DCalib::CXHist1DCalib(const TGCompositeFrame *MotherFrame, UInt_t w, UInt
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Calibration type"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,10,0,0));
fCalibType = new TGComboBox(fHorizontalFrame);
fHorizontalFrame->AddFrame(fCalibType,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX | kLHintsExpandY ,1,3,0,0));
fCalibType->AddEntry("Slope only (a1*x)",0);
fCalibType->AddEntry("Affine (a0 + a1*x)",1);
fCalibType->AddEntry("Parabolic (a0 + a1*x + a2*x*x)",2);
fCalibType->Select(0);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Calibration order"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,10,0,0));
fCalibOrder = new TGNumberEntry(fHorizontalFrame, 1, 5, 0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative);
fHorizontalFrame->AddFrame(fCalibOrder,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,0,0));
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "No offset"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,10,0,0));
fNoOffset = new TGCheckButton(fHorizontalFrame);
fNoOffset->SetState(kButtonUp);
fHorizontalFrame->AddFrame(fNoOffset,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,-1,0,0));
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
......@@ -221,8 +220,8 @@ void CXHist1DCalib::GetCurrentRange()
return;
}
fRangeMin->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetFirst()));
fRangeMax->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetLast()));
fRangeMin->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetFirst()));
fRangeMax->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetLast()));
}
void CXHist1DCalib::Calibrate()
......@@ -240,10 +239,10 @@ void CXHist1DCalib::Calibrate()
return;
}
if(fRangeMin->GetNumber()<hist->GetXaxis()->GetBinLowEdge(1))
fRangeMin->SetNumber(hist->GetXaxis()->GetBinLowEdge(1));
if(fRangeMax->GetNumber()>hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins()))
fRangeMax->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins()));
if(fRangeMin->GetNumber()<hist->GetXaxis()->GetBinLowEdge(1))
fRangeMin->SetNumber(hist->GetXaxis()->GetBinLowEdge(1));
if(fRangeMax->GetNumber()>hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins()))
fRangeMax->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetNbins()));
fRecalEnergy->Reset();
fRecalEnergy->SetDataFromHistTH1(hist,0);
......@@ -256,12 +255,8 @@ void CXHist1DCalib::Calibrate()
// fRecalEnergy->SetChannelOffset(0); // channel offset to subtract to the position of the peaks [0]
fRecalEnergy->SetVerbosityLevel(fVerboseLevel->GetNumber()-1); // verbosity -1=noprint 0=fit_details, 1=calib_details, 2=more_calib_details [-1]
if(fCalibType->GetSelected()==0)
fRecalEnergy->UseSlopeCalib();
else if(fCalibType->GetSelected()==1)
fRecalEnergy->UseAffineCalib();
else if(fCalibType->GetSelected()==2)
fRecalEnergy->UseParabolicCalib();
fRecalEnergy->SetFitPlynomialOrder(fCalibOrder->GetNumber());
fRecalEnergy->SetNoOffset(fNoOffset->GetState());
fRecalEnergy->UseLeftTail(fLeftTail->GetState());
fRecalEnergy->UseRightTail(fRightTail->GetState());
......@@ -280,12 +275,14 @@ void CXHist1DCalib::Calibrate()
vector < Fitted > FitResults = fRecalEnergy->GetFitResults();
if(fVerboseLevel->GetNumber()==0) {
cout<< setprecision(6) << hist->GetName()<<": ";
if(fCalibType->GetSelected() == 0) cout<< left << setw(10) << fRecalEnergy->slope*fRecalEnergy->hGain << endl;
else if(fCalibType->GetSelected() == 1) cout<< left << setw(10) << fRecalEnergy->offset1 << " " << setw(10) << fRecalEnergy->slope1*fRecalEnergy->hGain << endl;
else if(fCalibType->GetSelected() == 2) cout<< left << setw(10) << fRecalEnergy->offset2 << " " << setw(10) << fRecalEnergy->slope2*fRecalEnergy->hGain << " " << setw(10) << fRecalEnergy->curv2*fRecalEnergy->hGain*fRecalEnergy->hGain << endl;
cout<< left << scientific << setprecision(6);
cout<< hist->GetName()<<": ";
cout << setw(14) << fRecalEnergy->fCalibFunction->GetParameter(0);
for(int i=1 ; i<=fRecalEnergy->fCalibOrder ; i++) cout << setw(14) << fRecalEnergy->fCalibFunction->GetParameter(i)*TMath::Power(fRecalEnergy->hGain,i);
cout<<endl;
}
double xmin=-1;
double xmax=-1;
......@@ -351,6 +348,18 @@ void CXHist1DCalib::Calibrate()
fMainWindow->RefreshPads();
if(fCalibCanvas != nullptr) delete fCalibCanvas;
fCalibCanvas = new TCanvas;
fCalibCanvas->SetName("CalibrationResults");
fCalibCanvas->SetTitle("Calibration Results");
fCalibCanvas->Divide(1,2,0.0001,0.0001);
fCalibCanvas->cd(1);
fRecalEnergy->fCalibGraph->Draw("ap");
fRecalEnergy->fCalibFunction->Draw("same");
fCalibCanvas->cd(2);
fRecalEnergy->fResidueGraph->Draw("ape");
fMainWindow->GetCanvas()->cd();
}
double CXHist1DCalib::DinoFct(double*xx,double*pp)
......
......@@ -16,6 +16,7 @@ class TGNumberEntry;
class TF1;
class TGComboBox;
class TGCheckButton;
class CXCanvas;
class CXHist1DCalib : public TGVerticalFrame
{
......@@ -34,7 +35,8 @@ private:
TGNumberEntry *fThresholdSPEntry = nullptr;
TGNumberEntry *fVerboseLevel = nullptr;
TGComboBox *fCalibType = nullptr;
TGNumberEntry *fCalibOrder = nullptr;
TGCheckButton *fNoOffset = nullptr;
TGCheckButton *fLeftTail = nullptr;
TGCheckButton *fRightTail = nullptr;
......@@ -48,6 +50,7 @@ private:
Double_t fERef=0.;
GwRecalEnergy *fRecalEnergy = nullptr;
TCanvas *fCalibCanvas = nullptr;
public:
......
This diff is collapsed.
......@@ -18,6 +18,8 @@
#include "TString.h"
#include "GwFitSpek.h"
#include "TF1.h"
#include "TGraphErrors.h"
using namespace std;
......@@ -97,6 +99,13 @@ public:
Int_t fId;
// Calibration
Int_t fCalibOrder=1;
Bool_t fNoOffset = false;
TF1 *fCalibFunction = nullptr;
TGraphErrors *fCalibGraph = nullptr;
TGraphErrors *fResidueGraph = nullptr;
public:
void SetDataFromHistTH1(TH1 *hist, Int_t Id=0);
......@@ -123,10 +132,10 @@ public:
void UseLeftTail(bool on){useTL = on;} // disable using the Left Tail in peak fit
void UseRightTail(bool on){useTR = on;} // disable using the Right Tail in peak fit
void UseLinearCalib(){numZero++;} // add a fake peak at (0,0) to push calibration through origin
void UseSlopeCalib(){doSlope = true; doPoly1 = false; doPoly2 = false;} // affine calibration y = offset + slope*x
void UseAffineCalib(){doPoly1 = true; doPoly2 = false;} // affine calibration y = offset + slope*x
void UseParabolicCalib(){doPoly1 = true; doPoly2 = true;} //parabolic calibration y = offset + slope*x + curv*x*x
void SetFitPlynomialOrder(Int_t order){fCalibOrder = order;}
void SetNoOffset(Bool_t on) { fNoOffset = on;}
void SetGain(float Gain){hGain = Gain;} //scaling factor for the slope [1]
void SetVerbosityLevel(int Verb){Verbosity = Verb;} // verbosity bit0=fit_details, bit1=calib_details, bit2=more_calib_details [0]
......@@ -136,6 +145,8 @@ public:
Float_t GetOffset(){return offset1;}
Float_t GetSlope(){return slope1*hGain;}
Double_t PolynomialFunc(Double_t*xx,Double_t*pp);
public:
vector<Fitted> Peaks;
......@@ -207,7 +218,7 @@ public:
int xP_Next1(float *yVal, int yLen);
int xP_Next2(float *yVal, int yLen);
int FitPeaks(int verbose);
double ECalibration(double& offset1, double& slope1_, double& offset2_, double& slope2_, double& curv2_, int verbose_); // energy calibration from position-energy pair(s)
Int_t EROOTCalibration();
double TCalibration(); // shift of the largest peak from its reference position (to be done)
bool InvertMatrix3(const double m[9], double invOut[9]);
double Calibrated(double x);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment