Commit 65ba8bb5 authored by Jérémie Dudouet's avatar Jérémie Dudouet
Browse files

Add utility for energy calibration

parent 456ab31a
Pipeline #43843 passed with stage
in 9 minutes and 53 seconds
......@@ -41,7 +41,7 @@ LINK_DIRECTORIES( ${ROOT_LIBRARY_DIR} )
SET(EXTRA_EXTERNAL_LIBRARIES ${EXTRA_EXTERNAL_LIBRARIES} Matrix Rint)
# gw
set(EXTRA_INTERNAL_LIBRARIES GWPHYSICS )
set(EXTRA_INTERNAL_LIBRARIES GWPHYSICS GWTOOLS)
####################
### define files ###
......
53.156
79.623
80.999
160.609
223.116
276.404
302.858
356.014
383.859
475.36
563.27
569.30
604.68
795.78
801.86
1038.53
1167.89
1365.17
121.7793
244.6927
344.2724
411.111
443.979
778.890
964.014
1085.793
#1089.700
1112.070
1407.993
186.211
241.981
295.213
351.921
609.312
768.356
934.061
1120.287
1238.110
1377.669
1509.228
1729.595
1764.494
1847.420
2118.551
2204.215
2447.810
30.6382
83.899
212.017
271.175
831.45
846.764
983.02
1622.87
1778.969
1810.726
2282.773
2590.244
3033.89
3465.07
4133.408
4259.539
4733.847
6702.034
7213.034
7724.034
846.772
1037.840
1175.102
1238.282
1360.250
1771.351
2015.181
2034.755
2598.458
3201.962
3253.416
3272.990
3451.152
3547.925
#include "CXHist1DCalib.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include "TGNumberEntry.h"
#include "TGButton.h"
#include "TGLabel.h"
#include "TROOT.h"
#include "TObjArray.h"
#include "TF1.h"
#include "TMath.h"
#include "TGComboBox.h"
#include "CXMainWindow.h"
using namespace std;
CXHist1DCalib::CXHist1DCalib(const TGCompositeFrame *MotherFrame, UInt_t w, UInt_t h) : TGVerticalFrame(MotherFrame, w, h, kFixedWidth)
{
TGGroupFrame *fGroupFrame = new TGGroupFrame(MotherFrame, "Energy calibration", kVerticalFrame);
fGroupFrame->SetTextColor(CXblue);
fGroupFrame->SetTitlePos(TGGroupFrame::kLeft); // right aligned
AddFrame(fGroupFrame, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 3, 3, 0, 0));
fGroupFrame->AddFrame(new TGLabel(fGroupFrame,"Sources"), new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 0, 0, 5, 0));
fSources = new TGTextEntry(fGroupFrame, "");
fSources->SetToolTipText("List of sources (files in DataBase/Sources)");
fSources->Connect("TextChanged(const char *)", "CXHist1DCalib", this, "UpdateText()");
fSources->Connect("ReturnPressed()", "CXHist1DCalib", this, "UpdateSources()");
fGroupFrame->AddFrame(fSources,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,0,0));
TGCompositeFrame *fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Verbose level"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,20,0,0));
fVerboseLevel = new TGNumberEntry(fHorizontalFrame, 1, 3, 0, TGNumberFormat::kNESInteger, TGNumberFormat::kNEANonNegative ,TGNumberFormat::kNELLimitMinMax,0,3);
fHorizontalFrame->AddFrame(fVerboseLevel,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,5,5));
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);
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "FWHM"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,40,0,0));
fFWHMSPEntry = new TGNumberEntry(fHorizontalFrame, 15, 3, 0, TGNumberFormat::kNESReal, TGNumberFormat::kNEANonNegative ,TGNumberFormat::kNELNoLimits);
fHorizontalFrame->AddFrame(fFWHMSPEntry,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,0,0));
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Amplitude"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,15,10,0,0));
fThresholdSPEntry = new TGNumberEntry(fHorizontalFrame, 5, 4, 0, TGNumberFormat::kNESReal, TGNumberFormat::kNEANonNegative ,TGNumberFormat::kNELNoLimits);
fHorizontalFrame->AddFrame(fThresholdSPEntry,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,0,0));
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Range: From "), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,5,0,0));
fRangeMin = new TGNumberEntry(fHorizontalFrame, 0, 3, 0, TGNumberFormat::kNESReal, TGNumberFormat::kNEANonNegative ,TGNumberFormat::kNELNoLimits);
fHorizontalFrame->AddFrame(fRangeMin,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,-1,0,0));
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "To "), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,15,5,0,0));
fRangeMax = new TGNumberEntry(fHorizontalFrame, 32000, 4, 0, TGNumberFormat::kNESReal, TGNumberFormat::kNEANonNegative ,TGNumberFormat::kNELNoLimits);
fHorizontalFrame->AddFrame(fRangeMax,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,0,0));
TGTextButton *DoCurrentRange = new TGTextButton(fHorizontalFrame, "Current");
DoCurrentRange->Connect("Clicked()", "CXHist1DCalib", this, "GetCurrentRange()");
fHorizontalFrame->AddFrame(DoCurrentRange,new TGLayoutHints(kLHintsCenterY | kLHintsExpandX,5,3,0,0));
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Left tail"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,20,0,0));
fLeftTail = new TGCheckButton(fHorizontalFrame);
fLeftTail->SetState(kButtonDown);
fHorizontalFrame->AddFrame(fLeftTail,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,-1,0,0));
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Right tail"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,15,20,0,0));
fRightTail = new TGCheckButton(fHorizontalFrame);
fRightTail->SetState(kButtonDown);
fHorizontalFrame->AddFrame(fRightTail,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,0,0));
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
fHorizontalFrame->AddFrame(new TGLabel(fHorizontalFrame, "Use 2nd-derivative search"), new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,20,0,0));
f2DSearch = new TGCheckButton(fHorizontalFrame);
fHorizontalFrame->AddFrame(f2DSearch,new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX ,1,3,0,0));
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,-10,-10,5,5));
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
TGTextButton *DoCalib = new TGTextButton(fHorizontalFrame, "Calibrate");
DoCalib->SetTextColor(CXred);
DoCalib->Connect("Clicked()", "CXHist1DCalib", this, "Calibrate()");
fHorizontalFrame->AddFrame(DoCalib,new TGLayoutHints(kLHintsCenterY | kLHintsExpandX,5,10,0,0));
fGroupFrame->AddFrame(fHorizontalFrame,new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX,50,50,10,10));
fListOfObjects = new TList;
fListOfObjects->SetOwner();
fSourcesFolder = Form("%s/DataBase/Sources",getenv("GWSYS"));
fRecalEnergy = new GwRecalEnergy;
}
CXHist1DCalib::~CXHist1DCalib()
{
}
void CXHist1DCalib::SetMainWindow(CXMainWindow *w)
{
fMainWindow = w;
}
void CXHist1DCalib::HandleMouse(Int_t EventType,Int_t EventX,Int_t EventY, TObject* selected)
{
TPad *pad = dynamic_cast<TPad*>(gROOT->GetSelectedPad());
if(pad==nullptr)
return;
}
void CXHist1DCalib::HandleMyButton()
{
}
void CXHist1DCalib::UpdateText()
{
fSources->SetTextColor(CXred);
}
void CXHist1DCalib::UpdateSources()
{
fEnergies.clear();
fERef = 0.;
TString tmp = fSources->GetText();
tmp.ReplaceAll(",", " ");
TObjArray *arr = tmp.Tokenize(" ");
if(arr->GetEntries()==0) {
delete arr;
return;
}
fSources->SetTextColor(CXblack);
INFO_MESS << "Energies used for calibration: " << ENDL;
for(int i=0 ; i<arr->GetEntries() ; i++) {
TString FileName = Form("%s/%s.source",fSourcesFolder.Data(),arr->At(i)->GetName());
if(gSystem->IsFileInIncludePath(FileName)) {
ifstream file(FileName);
cout << "Source: " << arr->At(i)->GetName() << endl;
string line;
TString Buffer;
while(file) {
getline(file,line);
Buffer=line;
if(Buffer.BeginsWith("#"))
continue;
Buffer.ReplaceAll("\t"," ");
TObjArray *arr2 = Buffer.Tokenize(" ");
if(arr2->GetEntries()>0) {
bool ref = false;
TString textline = (TString)arr2->First()->GetName();
if(textline.Contains("*")) {
textline.ReplaceAll("*","");
ref = true;
}
Double_t E = ((TString)arr2->First()->GetName()).Atof();
fEnergies.push_back(E);
if(ref) fERef = E;
cout << " --> E: " << E << " keV" << endl;
}
delete arr2;
}
file.close();
}
else if(((TString)arr->At(i)->GetName()) == "-ener" && i<arr->GetEntries()-1) {
TString manualvalue = ((TString)arr->At(i+1)->GetName());
if(manualvalue.IsFloat()) {
fEnergies.push_back(manualvalue.Atof());
cout << "Value: " << arr->At(i+1)->GetName() << " (keV) manually added." << endl;
i++;
continue;
}
else {
WARN_MESS << " Error in manually added value" << ENDL;
fSources->SetTextColor(CXred);
continue;
}
}
else {
WARN_MESS << FileName << " not found " << ENDL;
fSources->SetTextColor(CXred);
}
}
delete arr;
if(fEnergies.size() && fERef>0.) {
INFO_MESS << "Reference energy for printouts: " << fERef << "keV" << ENDL;
}
cout<<endl;
}
void CXHist1DCalib::CleanCalib()
{
fListOfObjects->Clear();
}
void CXHist1DCalib::GetCurrentRange()
{
TH1 *hist = fMainWindow->GetCanvas()->FindHisto();
if(hist == nullptr || hist->GetDimension()>1) {
WARN_MESS << "No 1D histogram in the current pad, ignored " << ENDL;
return;
}
if(fRangeMin->GetNumber()<hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetFirst()))
fRangeMin->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetFirst()));
if(fRangeMax->GetNumber()>hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetLast()))
fRangeMax->SetNumber(hist->GetXaxis()->GetBinLowEdge(hist->GetXaxis()->GetLast()));
}
void CXHist1DCalib::Calibrate()
{
CleanCalib();
TH1 *hist = fMainWindow->GetCanvas()->FindHisto();
if(hist == nullptr || hist->GetDimension()>1) {
WARN_MESS << "No 1D histogram in the current pad, ignored " << ENDL;
return;
}
if(fEnergies.size()==0) {
WARN_MESS << "No source defined, ignored " << ENDL;
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()));
fRecalEnergy->Reset();
fRecalEnergy->SetDataFromHistTH1(hist,0);
for (auto ie : fEnergies)
fRecalEnergy->AddPeak(ie);
if(fERef>0.) fRecalEnergy->SetRefPeak(fERef);
fRecalEnergy->SetGain(1.); // scaling factor for the slope [1]
// 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->UseLeftTail(fLeftTail->GetState());
fRecalEnergy->UseRightTail(fRightTail->GetState());
if(f2DSearch->GetState() == kButtonUp)
fRecalEnergy->UseFirstDerivativeSearch();
else
fRecalEnergy->UseSecondDerivativeSearch();
//fRecalEnergy->UseSecondDerivativeSearch(); // use the 2nd-derivative search
fRecalEnergy->SetGlobalChannelLimits(fRangeMin->GetNumber(),fRangeMax->GetNumber()); // limit the search to this range in channels
fRecalEnergy->SetGlobalPeaksLimits(fFWHMSPEntry->GetNumber(),fFWHMSPEntry->GetNumber()); // default fwhm and minmum amplitude for the peaksearch [15 5]
fRecalEnergy->StartCalib();
vector < Fitted > FitResults = fRecalEnergy->GetFitResults();
double xmin=-1;
double xmax=-1;
int NGoodPeak=0;
for(size_t i=0 ; i<FitResults.size() ; i++)
{
Fitted FitRes = FitResults[i];
if(!FitRes.good) continue;
NGoodPeak++;
// Calc n good sub peaks
int NSubPeaks = FitRes.NSubPeaks;
int NGoodSubPeaks = 0;
for(int j=0 ; j<NSubPeaks ; j++) {
FitRes = FitResults[i+j];
if(FitRes.good) NGoodSubPeaks++;
}
FitRes = FitResults[i];
TF1 *f = GetDinoFct(Form("Peak%d_%.1f",NGoodPeak,FitRes.eref),FitRes.BgFrom,FitRes.BgTo,5+6*NGoodSubPeaks);
f->SetNpx(10000);
f->SetParameter(0,NGoodSubPeaks);
f->SetParameter(1,FitRes.BgFrom);
f->SetParameter(2,FitRes.BgTo);
f->SetParameter(3,FitRes.BgdOff);
f->SetParameter(4,FitRes.BgdSlope);
if(xmin == -1) xmin = FitRes.BgFrom;
xmax = FitRes.BgTo;
int peakid=0;
for(int j=0 ; j<NSubPeaks ; j++) {
FitRes = FitResults[i+j];
if(!FitRes.good) continue;
f->SetParameter(5+peakid*6+0,FitRes.ampli);
f->SetParameter(5+peakid*6+1,FitRes.posi);
f->SetParameter(5+peakid*6+2,FitRes.fwhm);
f->SetParameter(5+peakid*6+3,FitRes.Lambda);
f->SetParameter(5+peakid*6+4,FitRes.Rho);
f->SetParameter(5+peakid*6+5,FitRes.S);
if(fLeftTail->GetState()==kButtonUp)
f->SetParameter(5+peakid*6+3,-50);
if(fRightTail->GetState()==kButtonUp)
f->SetParameter(5+peakid*6+4,50);
peakid++;
}
f->Draw("same");
fListOfObjects->Add(f);
i += NSubPeaks;
}
hist->GetXaxis()->SetRangeUser(xmin-(xmax-xmin)*0.1,xmax+(xmax-xmin)*0.1);
fMainWindow->RefreshPads();
}
double CXHist1DCalib::DinoFct(double*xx,double*pp)
{
double x = xx[0];
int NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
double BgFrom = pp[1]; //First Channel for the Bg estimation
double BgTo = pp[2]; //Last Channel for the Bg estimation
double BgdOff = pp[3]; //Bg offset
double BgdSlope = pp[4]; //Bg slope
double f_tot = 0.;
if(x<BgFrom || x>BgTo) return 0.;
else
{
double BGd = BgdSlope*(x-BgFrom) + BgdOff;
f_tot += BGd;
}
// cout<<NSubPeaks<<" "<<BgFrom<<" "<<BgTo<<" "<<BgdOff<<" "<<BgdSlope<<endl;
for(int i=0 ; i<NSubPeaks ; i++)
{
double Ampli = pp[5+i*6+0];
double Mean = pp[5+i*6+1];
double Sigma = pp[5+i*6+2]*1./sqrt(8.*log(2.));;
double Lambda = pp[5+i*6+3];
double Rho = pp[5+i*6+4];
double S = pp[5+i*6+5];
// cout<<Ampli<<" "<<Mean<<" "<<Sigma<<" "<<Lambda<<" "<<Rho<<" "<<S<<endl;
double U = (x-Mean)/Sigma;
double f_g = Ampli*TMath::Exp(-U*U*0.5);
double f_lambda = Ampli*TMath::Exp(-0.5*Lambda*(2.*U-Lambda));
double f_rho = Ampli*TMath::Exp(-0.5*Rho*(2.*U-Rho));
double f_S = Ampli*S*1./((1+TMath::Exp(U))*(1+TMath::Exp(U)));
if(U<Lambda) f_tot += f_lambda;
else if(U>Rho) f_tot += f_rho;
else f_tot += f_g;
f_tot += f_S;
}
return f_tot;
}
///******************************************************************************************///
TF1 *CXHist1DCalib::GetDinoFct(TString Name,double min, double max, int Npar)
{
TF1 *f = new TF1(Name,this,&CXHist1DCalib::DinoFct,min,max,Npar,"CXHist1DCalib","DinoFct");
return f;
}
ClassImp(CXHist1DCalib)
#ifndef CXHist1DCalib_H
#define CXHist1DCalib_H
#include "RQ_OBJECT.h"
#include "TGFrame.h"
#include "TArrow.h"
#include "TObject.h"
#include "GwRecalEnergy.h"
using namespace std;
class TGTextEntry;
class CXMainWindow;
class TGNumberEntry;
class TF1;
class TGComboBox;
class TGCheckButton;
class CXHist1DCalib : public TGVerticalFrame
{
RQ_OBJECT("CXHist1DCalib");
public:
private:
TGTextEntry *fSources = nullptr;
CXMainWindow *fMainWindow = nullptr;
TGNumberEntry *fRangeMin = nullptr;
TGNumberEntry *fRangeMax = nullptr;
TGNumberEntry *fFWHMSPEntry = nullptr;
TGNumberEntry *fThresholdSPEntry = nullptr;
TGNumberEntry *fVerboseLevel = nullptr;
TGComboBox *fCalibType = nullptr;
TGCheckButton *fLeftTail = nullptr;
TGCheckButton *fRightTail = nullptr;
TGCheckButton *f2DSearch = nullptr;
TList *fListOfObjects = nullptr;