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

Correct errors on calibration procedure

Correct a bug in fit result
Add possibility to select the minimizer
parent bd57b806
......@@ -43,10 +43,15 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib_)
set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib_)
set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin_)
#if(APPLE)
# set(CMAKE_MACOSX_RPATH 1)
#endif()
#set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
#cmake_policy(SET CMP0077 NEW)
#cmake_policy(SET CMP0068 NEW)
# FIX Warning with rpath on macosx
if(APPLE)
set(CMAKE_MACOSX_RPATH 1)
set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
endif()
#
# ROOT MANDATORY
#
......
......@@ -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})
#####################
......
How to create/update the database
1) dowload, in the DataBase/LevelScheme/Downloads folder (create it if not existing), the zip files (3 files) from: https://www.nndc.bnl.gov/ensarchivals/
1) download, in the DataBase/LevelScheme/Downloads folder (create it if not existing), the zip files (3 files) from: https://www.nndc.bnl.gov/ensarchivals/
ex: for the data from 01/11/2018
ex: for the data from 04/10/2019
wget https://www.nndc.bnl.gov/ensarchivals/distributions/dist18/ensdf_181101_099.zip -P Downloads/
wget https://www.nndc.bnl.gov/ensarchivals/distributions/dist18/ensdf_181101_199.zip -P Downloads/
wget https://www.nndc.bnl.gov/ensarchivals/distributions/dist18/ensdf_181101_299.zip -P Downloads/
mkdir -p $GWSYS/DataBase/LevelScheme/Downloads
cd $GWSYS/DataBase/LevelScheme/
wget https://www.nndc.bnl.gov/ensarchivals/distributions/dist19/ensdf_191004_099.zip -P Downloads/
wget https://www.nndc.bnl.gov/ensarchivals/distributions/dist19/ensdf_191004_199.zip -P Downloads/
wget https://www.nndc.bnl.gov/ensarchivals/distributions/dist19/ensdf_191004_300.zip -P Downloads/
2) execute the ExtractFiles script, giving the files to process in argument:
chmod u+w src/compile.sh
src/compile.sh
chmod u+x src/compile.sh
./src/compile.sh
./ExtractFiles Downloads/ensdf_181101_*.zip
......@@ -2,4 +2,4 @@
NAME=ExtractFiles
NAME2=src/ExtractFiles.cpp
g++ $NAME2 -o $NAME -O6 -g `root-config --cflags` `root-config --cflags --libs` -lProof
g++ $NAME2 -o $NAME -O3 -g `root-config --cflags` `root-config --cflags --libs` -lProof
......@@ -4,22 +4,30 @@
#include "TBox.h"
#include "TH1.h"
#include "TVirtualPad.h"
#include "TFrame.h"
#include "CXFit.h"
#include "CXBgdFit.h"
CXArrow::CXArrow(CXFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize,Option_t *option) : TArrow(E, y1, E, y2, arrowsize, option)
CXArrow::CXArrow(CXFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize, Float_t textsize, Option_t *option) : TArrow(E, y1, E, y2, arrowsize, option)
{
fFit = fit;
fTextSize = textsize;
}
CXArrow::~CXArrow()
CXArrow::CXArrow(CXBgdFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize, Float_t textsize, Option_t *option) : TArrow(E, y1, E, y2, arrowsize, option)
{
fBgdFit = fit;
fTextSize = textsize;
}
void CXArrow::SetEnergy(Float_t E)
{
SetX1(E);
SetX2(E);
if(fFit) fFit->Update();
if(fBgdFit) fBgdFit->Update();
}
void CXArrow::Set(Double_t X, Double_t Y1, Double_t Y2)
......@@ -32,33 +40,34 @@ void CXArrow::Set(Double_t X, Double_t Y1, Double_t Y2)
void CXArrow::RemoveArrow()
{
if(fFit)
fFit->RemoveArrow(this);
if(fFit) fFit->RemoveArrow(this);
if(fBgdFit) fBgdFit->RemoveArrow(this);
}
void CXArrow::RemoveFit()
{
if(fFit)
delete fFit;
delete fFit;
delete fBgdFit;
}
Int_t CXArrow::Compare(const TObject *obj) const
{
if(obj->InheritsFrom(CXArrow::Class())) {
CXArrow *arr = (CXArrow*)obj;
auto *arr = dynamic_cast<const CXArrow*>(obj);
return (GetEnergy() > arr->GetEnergy()) ? 1 : -1;
}
else return 0;
return 0;
}
void CXArrow::SetText(TH1 *hist, TString text, TString tooltip)
void CXArrow::SetText(TH1 *hist, const TString &text, const TString &tooltip)
{
Int_t XMinbin = hist->GetXaxis()->GetFirst();
Int_t XMaxbin = hist->GetXaxis()->GetLast();
Int_t XMin = hist->GetBinLowEdge(XMinbin);
Int_t XMax = hist->GetBinLowEdge(XMaxbin);
if(fBox) delete fBox;
delete fBox;
fBox = new CXArrowBox(this);
fBox->SetX1(fX1-(XMax-XMin)/300.);
......@@ -72,12 +81,14 @@ void CXArrow::SetText(TH1 *hist, TString text, TString tooltip)
fBox->SetLineColor(0);
fBox->SetLineStyle(3);
if(fLatex) delete fLatex;
fLatex = new TLatex(fX1,fY2,text);
Double_t height = gPad->GetFrame()->GetY2() - gPad->GetFrame()->GetY1();
delete fLatex;
fLatex = new TLatex(fX1,fY2+fArrowSize*height,text);
fLatex->SetTextAngle(90);
fLatex->SetTextFont(132);
fLatex->SetTextSize(fTextSize);
fLatex->SetTextColor(hist->GetLineColor());
fLatex->SetTextAlign(12);
fLatex->SetBit(TObject::kCannotPick);
fLatex->Draw();
......
......@@ -5,6 +5,8 @@
#include "TBox.h"
class CXFit;
class CXBgdFit;
class TLatex;
class TBox;
class TH1;
......@@ -19,6 +21,7 @@ private:
TList *fList = nullptr;
CXFit *fFit = nullptr;
CXBgdFit *fBgdFit = nullptr;
TLatex *fLatex = nullptr;
CXArrowBox *fBox = nullptr;
......@@ -27,10 +30,13 @@ private:
public:
CXArrow(CXFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize=0.05,Option_t *option=">");
~CXArrow();
CXArrow(CXFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize=0.05, Float_t textsize = 0.03, Option_t *option=">");
CXArrow(CXBgdFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize=0.05, Float_t textsize = 0.03, Option_t *option=">");
~CXArrow() = default;
CXFit *GetFit(){return fFit;}
CXBgdFit *GetBgdFit(){return fBgdFit;}
//! Energy
void SetEnergy(Float_t E); // *MENU* *ARGS={E=>fX1}*
......@@ -41,7 +47,7 @@ public:
void RemoveArrow(); // *MENU*
void RemoveFit(); // *MENU*
void SetText(TH1 *hist, TString text, TString tooltip);
void SetText(TH1 *hist, const TString &text, const TString &tooltip);
void ClearPad(TVirtualPad *pad = nullptr, bool refresh = true);
......
#include "CXBgdFit.h"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include "TMath.h"
#include "TF1.h"
#include "TH1.h"
#include "TFitResultPtr.h"
#include "TFitResult.h"
#include "TPad.h"
#include "TGNumberEntry.h"
#include "TGButton.h"
#include "Fit/Fitter.h"
#include "CXHist1DPlayer.h"
#include "CXArrow.h"
#include "CXMainWindow.h"
using namespace std;
CXBgdFit::CXBgdFit(TH1 *hist, TVirtualPad *pad, CXHist1DPlayer *player) : TObject()
{
fHistogram = hist;
fPad = pad;
fPlayer = player;
fListOfArrows = new TList;
fListOfArrows->SetOwner();
}
CXBgdFit::~CXBgdFit()
{
Clear(fPad);
delete fListOfArrows;
delete fBackFunction;
fPlayer->GetMainWindow()->RefreshPads();
}
void CXBgdFit::AddArrow(Double_t Energy)
{
Int_t Bin = fHistogram->FindBin(Energy);
Double_t Value = fHistogram->GetBinContent(Bin);
for(int i=fHistogram->FindBin(Energy)-2 ; i<=fHistogram->FindBin(Energy)+2 ; i++){
if(i>0 && fHistogram->GetBinContent(i)>Value){
Value = fHistogram->GetBinContent(i);
Bin = i;
}
}
Energy = fHistogram->GetBinCenter(Bin);
Double_t MaxGlob = fHistogram->GetMaximum();
auto *arrow = new CXArrow(this,Energy,(Value + MaxGlob/100.) ,(Value +MaxGlob/15.),0.01,0.03,"<|");
arrow->SetAngle(30);
arrow->SetLineWidth(2);
arrow->Draw();
fListOfArrows->Add(arrow);
Update();
}
void CXBgdFit::RemoveArrow(CXArrow *arrow)
{
if(arrow == nullptr) {
fListOfArrows->RemoveLast();
fPad->GetListOfPrimitives()->RemoveLast();
}
else {
fListOfArrows->Remove(arrow);
fPad->GetListOfPrimitives()->Remove(arrow);
}
Update();
}
void CXBgdFit::Update()
{
fBackgd.clear();
fListOfArrows->Sort();
TList back,Ener;
if(!fPlayer->DoNewBgdFit && (fListOfArrows->GetEntries()%2)==1) {
Clear(fPad);
WARN_MESS<<"Setts of two ranges are needed for a bgd fit --> command ignored"<<ENDL;
return;
}
for(int i=0 ; i<fListOfArrows->GetEntries() ; i++) {
auto *arr = dynamic_cast<CXArrow*>(fListOfArrows->At(i));
back.Add(arr);
}
for(int i=0 ; i<back.GetEntries() ; i++) {
auto *arr = dynamic_cast<CXArrow*>(back.At(i));
Double_t E = arr->GetEnergy();
Double_t MaxGlob = fHistogram->GetMaximum();
Double_t Value = fHistogram->GetBinContent(fHistogram->FindBin(E));
arr->Set(E,(Value + MaxGlob/100.), (Value+MaxGlob/15.));
arr->SetLineColor(kBlue);
arr->SetFillColor(kBlue);
fBackgd.push_back(E);
cout<<E<<endl;
}
fPlayer->GetMainWindow()->RefreshPads();
}
void CXBgdFit::Clear(TVirtualPad *pad)
{
fPlayer->EndFit();
if(pad == nullptr) pad = fPad;
TList *list = pad->GetListOfPrimitives();
list->Remove(fBackFunction);
for(int i=0 ; i<fListOfArrows->GetEntries() ; i++)
list->Remove(fListOfArrows->At(i));
fPlayer->RemoveBgdFit(this);
fPlayer->GetMainWindow()->RefreshPads();
}
void CXBgdFit::Fit()
{
if(fListOfArrows->GetEntries()<4 || (fListOfArrows->GetEntries()%2) !=0) {
return;
}
if(fPad==nullptr) {
WARN_MESS<<"No selected pad, ignored"<<ENDL;
return;
}
fPad->cd();
if(fPlayer == nullptr) {
WARN_MESS<<"1DPlayer not defined, ignored"<<ENDL;
return;
}
if(fHistogram==nullptr || fHistogram->InheritsFrom("TH2")) {
WARN_MESS<<"No 1D histogram found, ignored"<<ENDL;
return;
}
ROOT::Math::MinimizerOptions::SetDefaultMinimizer(fPlayer->GetMinimizer(),fPlayer->GetAlgorithm());
ROOT::Math::MinimizerOptions::SetDefaultTolerance(fPlayer->GetTolerance());
ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(fPlayer->GetPrintLevel());
Int_t NPars = 3;
delete fBackFunction;
fBackFunction = new TF1("MyFit", this, &CXBgdFit::FuncBackground, fBackgd.front(), fBackgd.back(), NPars, "CXBgdFit", "FuncBackground");
fBackFunction->SetParName(0, "BkgConst");
fBackFunction->SetParName(1, "BkgSlope");
fBackFunction->SetParName(2, "BkgExp");
fBackFunction->SetNpx(1000);
fBackFunction->SetLineColor(kRed);
// Copy only the ranges contains within the arrows
TH1 *HistoToFit = dynamic_cast<TH1*>(fHistogram->Clone());
HistoToFit->Reset();
for(size_t i=0; i<fBackgd.size() ; i+=2) {
Int_t binmin = fHistogram->GetXaxis()->FindBin(fBackgd.at(i));
Int_t binmax = fHistogram->GetXaxis()->FindBin(fBackgd.at(i+1));
for(int ibin=binmin ; ibin<=binmax ; ibin++) HistoToFit->SetBinContent(ibin,fHistogram->GetBinContent(ibin));
}
HistoToFit->GetXaxis()->SetRangeUser(fBackgd.front(), fBackgd.back());
//Calc Bckd
fBackFunction->SetParameter(0, HistoToFit->GetBinContent(HistoToFit->FindBin(fBackgd.front())));
fBackFunction->SetParLimits(0, 0.,HistoToFit->GetMaximum());
fBackFunction->SetParameter(1, 0);
fBackFunction->SetParLimits(1, -50., 0.);
fBackFunction->SetParameter(2, 0.);
// fBackFunction->SetParLimits(2, -1., 0.);
if(!fPlayer->fUseBgdPol1)
fBackFunction->FixParameter(1,0);
if(!fPlayer->fUseBgdExp)
fBackFunction->FixParameter(2,0);
TString FitOpt = "R0S";
if(fPlayer->GetPrintLevel()>0) FitOpt +="V";
TFitResultPtr r = HistoToFit->Fit(fBackFunction,FitOpt.Data());
ostringstream text;
text << "Fit results :";
cout<<text.str()<<endl;fPlayer->PrintInListBox(text.str(),kPrint);text.str("");
text << "Status: ";
if(r->Status()==0)
text << " Successeful" << endl;
else
text << " Failed" << endl;
cout<<text.str();
if(r->Status()==0)
fPlayer->PrintInListBox(text.str(),kPrint);
else
fPlayer->PrintInListBox(text.str(),kError);
text.str("");
Float_t Area = fHistogram->Integral(fHistogram->GetXaxis()->FindBin(fBackgd.front()),fHistogram->GetXaxis()->FindBin(fBackgd.back()));
Float_t BgdArea = fBackFunction->Integral(fBackgd.front(),fBackgd.back(),1e-6);
Float_t CorrArea = Area-BgdArea;
text<<left<<setw(11)<<"Integral"<<": "<<setprecision(7)<<setw(10)<<Area<<" ("<<setprecision(7)<<setw(10)<<sqrt(Area)<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"Bgd Area"<<": "<<setprecision(7)<<setw(10)<<BgdArea<<" ("<<setprecision(7)<<setw(10)<<sqrt(BgdArea)<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"Peak Area"<<": "<<setprecision(7)<<setw(10)<<CorrArea<<" ("<<setprecision(7)<<setw(10)<<sqrt(CorrArea)<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
fBackFunction->Draw("same");
fPlayer->GetMainWindow()->RefreshPads();
delete HistoToFit;
}
Double_t CXBgdFit::FuncBackground(Double_t*xx,Double_t*pp)
{
Double_t x = xx[0];
Double_t Back_const = pp[0];
Double_t Back_slope = pp[1];
Double_t Back_Exp = pp[2];
Double_t f_tot = (Back_const + (x-fBackgd.front())*Back_slope)*exp((x-fBackgd.front())*Back_Exp);
return f_tot;
}
#ifndef CXBgdFit_H
#define CXBgdFit_H
#include "TObject.h"
#include <vector>
class CXHist1DPlayer;
class TVirtualPad;
class TH1;
class TF1;
class CXArrow;
class CXBgdFit : public TObject
{
private:
CXHist1DPlayer *fPlayer = nullptr;
TList *fListOfArrows = nullptr;
TVirtualPad *fPad = nullptr;
TH1 *fHistogram = nullptr;
TF1 *fBackFunction = nullptr;
std::vector<Double_t> fBackgd;
public:
CXBgdFit(TH1 *hist, TVirtualPad *pad, CXHist1DPlayer *player);
~CXBgdFit();
void AddArrow(Double_t Energy);
void RemoveArrow(CXArrow *arrow = nullptr);
void Update();
void Fit();
Double_t FuncBackground(Double_t*xx,Double_t*pp);
void Clear(TVirtualPad *pad = nullptr);
ClassDef(CXBgdFit,0);
};
#endif
</
......@@ -21,15 +21,14 @@ CXBgdUtility::CXBgdUtility(const TGCompositeFrame *MotherFrame, UInt_t w, UInt_t
fListOfButtons1D = new TList();
TGLayoutHints *GroupHints = new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX, -10, -10, 0, 0);
TGLayoutHints *HFHints = new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,0,8,0);
auto *GroupHints = new TGLayoutHints(kLHintsCenterY | kLHintsLeft | kLHintsExpandX, -10, -10, 0, 0);
auto *HFHints = new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,0,8,0);
TGGroupFrame *f1DGroupFrame = new TGGroupFrame(MotherFrame, "1D background", kVerticalFrame);
auto *f1DGroupFrame = new TGGroupFrame(MotherFrame, "1D background", kVerticalFrame);
f1DGroupFrame->SetTextColor(red);
f1DGroupFrame->SetTitlePos(TGGroupFrame::kLeft); // right aligned
AddFrame(f1DGroupFrame, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 0, 0, 0, 0));
//Activation du mode background supress
fGroupFrame = new TGGroupFrame(f1DGroupFrame, "Activate utility", kVerticalFrame);
......@@ -142,16 +141,15 @@ CXBgdUtility::CXBgdUtility(const TGCompositeFrame *MotherFrame, UInt_t w, UInt_t
fGroupFrame->AddFrame(fHorizontalFrame,HFHints);
fHorizontalFrame = new TGCompositeFrame(fGroupFrame, 60, 20, kHorizontalFrame);
for(int i=0 ; i<7 ; i++)
{
fSmoothingButton1D[i] = new TGRadioButton(fHorizontalFrame, Form("%d",3+2*i), i);
if(i==0) fSmoothingButton1D[i]->SetState(kButtonDown);
fSmoothingButton1D[i]->Connect("Clicked()", "CXBgdUtility", this, "HandleSmoothButtons()");
fSmoothingButton1D[i]->Connect("Clicked()","CXBgdUtility", this, "GetParams()");
if(i==0) fHorizontalFrame->AddFrame(fSmoothingButton1D[i],new TGLayoutHints(kLHintsCenterY | kLHintsLeft,20,0,0,0));
else fHorizontalFrame->AddFrame(fSmoothingButton1D[i],new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,0,0,0));
fListOfButtons1D->Add(fSmoothingButton1D[i]);
for(int i=0 ; i<7 ; i++) {
fSmoothingButton1D.at(i) = new TGRadioButton(fHorizontalFrame, Form("%d",3+2*i), i);
if(i==0) fSmoothingButton1D.at(i)->SetState(kButtonDown);
fSmoothingButton1D.at(i)->Connect("Clicked()", "CXBgdUtility", this, "HandleSmoothButtons()");
fSmoothingButton1D.at(i)->Connect("Clicked()","CXBgdUtility", this, "GetParams()");
if(i==0) fHorizontalFrame->AddFrame(fSmoothingButton1D.at(i),new TGLayoutHints(kLHintsCenterY | kLHintsLeft,20,0,0,0));
else fHorizontalFrame->AddFrame(fSmoothingButton1D.at(i),new TGLayoutHints(kLHintsCenterY | kLHintsLeft,5,0,0,0));
fListOfButtons1D->Add(fSmoothingButton1D.at(i));
}
fGroupFrame->AddFrame(fHorizontalFrame,HFHints);
......@@ -185,11 +183,11 @@ CXBgdUtility::CXBgdUtility(const TGCompositeFrame *MotherFrame, UInt_t w, UInt_t
fListOfButtons1D->Add(fSubtractButton);
SetButtonsStatus(0);
SetButtonsStatus(false);
/// Background 2D
TGGroupFrame *f2DGroupFrame = new TGGroupFrame(MotherFrame, "2D background", kVerticalFrame);
auto *f2DGroupFrame = new TGGroupFrame(MotherFrame, "2D background", kVerticalFrame);
f2DGroupFrame->SetTextColor(red);
f2DGroupFrame->SetTitlePos(TGGroupFrame::kLeft); // right aligned
AddFrame(f2DGroupFrame, new TGLayoutHints(kLHintsTop | kLHintsLeft | kLHintsExpandX, 0, 0, 0, 0));
......@@ -258,10 +256,7 @@ CXBgdUtility::CXBgdUtility(const TGCompositeFrame *MotherFrame, UInt_t w, UInt_t
fGroupFrame->AddFrame(fHorizontalFrame, GroupHints);
}
CXBgdUtility::~CXBgdUtility()
{
}
CXBgdUtility::~CXBgdUtility() = default;
void CXBgdUtility::DoSubtract()
{
......@@ -275,7 +270,7 @@ void CXBgdUtility::DoSubtract()
return;
}
TH1 *bkgrd = (TH1*) gPad->GetListOfPrimitives()->FindObject(Form("%s_BG",hist->GetName()));
TH1 *bkgrd = dynamic_cast<TH1*>(gPad->GetListOfPrimitives()->FindObject(Form("%s_BG",hist->GetName())));
if(bkgrd == nullptr)
{
......@@ -305,12 +300,12 @@ void CXBgdUtility::ToggleBckSupp()
if(fActivate->GetState() == kButtonDown)
{
SetButtonsStatus(1);
SetButtonsStatus(true);
GetParams();
}
else
{
SetButtonsStatus(0);
SetButtonsStatus(false);
for(int i=gPad->GetListOfPrimitives()->GetEntries()-1 ; i>=0 ; i--)
{
......@@ -333,7 +328,7 @@ void CXBgdUtility::SetButtonsStatus(bool on)
for(int i=1 ; i<fListOfButtons1D->GetEntries() ; i++)