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

Update Cubix

parent 2139fb07
ENCAL OUTPUT FILE
TMEUHS.STO EU152.SOU TM 159 EU152 HIGH SPIN CALIBRATION
2 0.000000000000D+00 1.000000000000D+00 0.000000000000D+00
0.000000000000D+00 0.000000000000D+00 0.000000000000D+00
### Cube file
CubeFileName MyCube.cub
### 2D projection file
2DProj MyCube.2dp
### LUT file
LUT MyCube.tab
### Contraction factor (default 1)
CompressFact 1
### Total projection file
TotalProj MyCube.spe
### Backgroud file
BackgroudFile MyCubeBG.spe
### Calibration file (not mandatory)
#CalFile MyCube.aca
### Efficiency file (not mandatory)
#EffFile MyCube.aef
#######################
### Cube Properties ###
#######################
# Name of the 3D cube file
CubeFileName MyCube.cub
# name of the look-up table file to be used to map ADC chs to cube chs
TabFile test.tab
# Size to be used for the "scratch" disk file, in MB
ScratchSize 512
# Compression factor (to obtain binning <1 keV, one needs to scale the energy)
# For example, to obtain a minimal 0.5 KeV/Channel, we fill the Cube on 8k channel for 4keV range ==> Compression factor = 2
# Must be similar to the value in the GenBinning conf file
CompressionFactor 1
################################
## Input Root Tree Properties ##
################################
#Name of the input TTree (regular expressions can be used as for building TChain)
TreeFile MyTree.root
#Name of the TTree
TreeName TreeName
#Events to process (0 => All)
NEvents 1000
#Name of the branch corresponding to the gamma-ray energy
EGamma EGamma
#Name of the branch corresponding to the gamma-ray muliplicity
GammaMult GammaMult
#Cut to be apllied (see TTreeFormulaExpression)
Cut None
EFFIT PARAMETER FILE 19-FEB-92 10:14:04
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+00
0.00000000E+00 0.00000000E+00 0.00000000E+00 0.00000000E+02 0.00000000E+03
#Name of the tab file
TabFile test.tab
#To fit the FWHM vs Energy function:
#Fit function: sqrt([0]*[0] + [1]*[1]*x/1000. + [2]*[2]*x/1000.*x/1000.)
#For each new point: FWHM Energy (keV) Value (keV)
FWHM 122. 1.35
FWHM 1333. 2.35
FWHM 2500. 3.00
# Compression factor (to obtain binning <1 keV, one needs to scale the energy)
# For example, to obtain a minimal 0.5 KeV/Channel, we fill the Cube on 8k channel for 4keV range ==> Compression factor = 2
CompressionFactor 1
#Number of bin per FWHM
Precision 1
#RangeMax (keV)
RangeMax 8192
#include "CXArrow.h"
#include "CXFit.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)
{
fFit = fit;
}
CXArrow::~CXArrow()
{
}
void CXArrow::SetEnergy(Float_t E)
{
SetX1(E);
SetX2(E);
}
void CXArrow::Set(Double_t X, Double_t Y1, Double_t Y2)
{
fX1=X;
fX2=X ;
fY1=Y1;
fY2=Y2;
}
void CXArrow::RemoveArrow()
{
if(fFit)
fFit->RemoveArrow(this);
}
void CXArrow::RemoveFit()
{
if(fFit)
delete fFit;
}
Int_t CXArrow::Compare(const TObject *obj) const
{
if(obj->InheritsFrom(CXArrow::Class())) {
CXArrow *arr = (CXArrow*)obj;
return (GetEnergy() > arr->GetEnergy()) ? 1 : -1;
}
else return 0;
}
ClassImp(CXArrow);
#ifndef CXArrow_H
#define CXArrow_H
#include "TArrow.h"
class CXFit;
class CXArrow : public TArrow
{
private:
TList *fList = nullptr;
CXFit *fFit = nullptr;
public:
CXArrow(CXFit *fit, Double_t E,Double_t y1 ,Double_t y2,Float_t arrowsize=0.05,Option_t *option=">");
~CXArrow();
CXFit *GetFit(){return fFit;}
//! Energy
void SetEnergy(Float_t E); // *MENU* *ARGS={E=>fX1}*
Double_t GetEnergy() const {return fX1;}
void Set(Double_t X, Double_t Y1, Double_t Y2);
void RemoveArrow(); // *MENU*
void RemoveFit(); // *MENU*
//! Sort
virtual Bool_t IsSortable() const {return kTRUE;}
virtual Int_t Compare(const TObject *obj) const;
ClassDef(CXArrow,0);
};
#endif
#include "CXFit.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 "CXHist1DPlayer.h"
#include "CXArrow.h"
#include "CXMainWindow.h"
using namespace std;
CXFit::CXFit(TH1 *hist, TVirtualPad *pad, CXHist1DPlayer *player) : TObject()
{
fHistogram = hist;
fPad = pad;
fPlayer = player;
fListOfArrows = new TList;
fListOfArrows->SetOwner();
fListOfPeaks = new TList;
fListOfPeaks->SetOwner();
}
CXFit::~CXFit()
{
Clear(fPad);
if(fListOfArrows) delete fListOfArrows;
if(fFitFunction) delete fFitFunction;
if(fBackFunction) delete fBackFunction;
if(fResidue) delete fResidue;
if(fListOfPeaks) delete fListOfPeaks;
fPlayer->GetMainWindow()->RefreshPads();
}
void CXFit::AddArrow(Double_t Energy)
{
Int_t Bin = fHistogram->FindBin(Energy);
Float_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);
Float_t MaxGlob = fHistogram->GetMaximum();
CXArrow *arrow = new CXArrow(this,Energy,(Value + MaxGlob/100.) ,(Value +MaxGlob/15.),0.01,"<|");
arrow->SetAngle(30);
arrow->SetLineWidth(2);
arrow->Draw();
fListOfArrows->Add(arrow);
Update();
}
void CXFit::RemoveArrow(CXArrow *arrow)
{
if(arrow == nullptr) {
fListOfArrows->RemoveLast();
fPad->GetListOfPrimitives()->RemoveLast();
}
else {
fListOfArrows->Remove(arrow);
fPad->GetListOfPrimitives()->Remove(arrow);
}
Update();
}
void CXFit::Update()
{
fEnergies.clear();
fBackgd.clear();
fListOfArrows->Sort();
TList back,Ener;
if(fPlayer->DoNewFit == false && fListOfArrows->GetEntries()<3) {
Clear(fPad);
WARN_MESS<<"At least one peak and two background are needed for a fit --> command ignored"<<ENDL;
return;
}
for(int i=0 ; i<fListOfArrows->GetEntries() ; i++) {
CXArrow *arr = dynamic_cast<CXArrow*>(fListOfArrows->At(i));
if(i==0) back.Add(arr);
else if(i==fListOfArrows->GetEntries()-1) back.Add(arr);
else Ener.Add(arr);
}
for(uint i=0 ; i<back.GetEntries() ; i++) {
CXArrow *arr = (CXArrow*)back.At(i);
Double_t E = arr->GetEnergy();
Float_t MaxGlob = fHistogram->GetMaximum();
Float_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);
}
for(uint i=0 ; i<Ener.GetEntries() ; i++) {
CXArrow *arr = (CXArrow*)Ener.At(i);
Double_t E = arr->GetEnergy();
Float_t MaxGlob = fHistogram->GetMaximum();
Float_t Value = fHistogram->GetBinContent(fHistogram->FindBin(E));
arr->Set(E,(Value + MaxGlob/100.), (Value+MaxGlob/15.));
arr->SetLineColor(kRed);
arr->SetFillColor(kRed);
fEnergies.push_back(E);
}
fPlayer->GetMainWindow()->RefreshPads();
}
void CXFit::Clear(TVirtualPad *pad)
{
fPlayer->EndFit();
if(pad == nullptr) pad = fPad;
TList *list = pad->GetListOfPrimitives();
list->Remove(fFitFunction);
list->Remove(fBackFunction);
list->Remove(fResidue);
for(int i=0 ; i<fListOfPeaks->GetEntries() ; i++)
list->Remove(fListOfPeaks->At(i));
for(int i=0 ; i<fListOfArrows->GetEntries() ; i++)
list->Remove(fListOfArrows->At(i));
fPlayer->RemoveFit(this);
fPlayer->GetMainWindow()->RefreshPads();
}
void CXFit::Fit()
{
if(fListOfArrows->GetEntries()<3) {
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;
}
// fPlayer->GetFitResultsBox()->RemoveAll();
Float_t DefFWHM = fPlayer->fNE_FWHM[0]->GetNumber();
Float_t DefFWHM_min = fPlayer->fNE_FWHM[1]->GetNumber();
Float_t DefFWHM_max = fPlayer->fNE_FWHM[2]->GetNumber();
Float_t LeftTailVal = fPlayer->fNE_LT[0]->GetNumber();;
Float_t LeftTailValMin = fPlayer->fNE_LT[1]->GetNumber();;
Float_t LeftTailValMax = fPlayer->fNE_LT[2]->GetNumber();;
Float_t RightTailVal = fPlayer->fNE_RT[0]->GetNumber();;
Float_t RightTailValMin = fPlayer->fNE_RT[1]->GetNumber();;
Float_t RightTailValMax = fPlayer->fNE_RT[2]->GetNumber();;
Float_t StepVal = 0.01;
Float_t StepValMin = -1.;
Float_t StepValMax = 1.;
Int_t NPars = 3+6*fEnergies.size();
if(fFitFunction) delete fFitFunction;
fFitFunction = new TF1("MyFit", this, &CXFit::DoubleTailedStepedGaussian, fBackgd[0], fBackgd[1], NPars, "CXFit", "DoubleTailedStepedGaussian");
fFitFunction->SetParName(0, "NumberOfPeaks");
fFitFunction->SetParName(1, "BkgConst");
fFitFunction->SetParName(2, "BkgSlope");
for(uint i=0 ; i<fEnergies.size() ; i++) {
fFitFunction->SetParName(3+i*6+0, Form("Height_%d",i));
fFitFunction->SetParName(3+i*6+1, Form("Position_%d",i));
fFitFunction->SetParName(3+i*6+2, Form("FWHM_%d",i));
fFitFunction->SetParName(3+i*6+3, Form("LeftTail_%d",i));
fFitFunction->SetParName(3+i*6+4, Form("RightTail_%d",i));
fFitFunction->SetParName(3+i*6+5, Form("AmplitudeStep_%d",i));
}
fFitFunction->SetNpx(1000);
fFitFunction->SetLineColor(kRed);
fFitFunction->FixParameter(0, fEnergies.size()); // 1 peak
Float_t x,y;
x = fPad->GetUxmin();
y = fPad->GetUxmax();
fHistogram->GetXaxis()->SetRangeUser(fBackgd[0],fBackgd[1]);
//Calc Bckd
fFitFunction->SetParameter(1, fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[0])));
fFitFunction->SetParLimits(1, fHistogram->GetMinimum(),fHistogram->GetMaximum());
fFitFunction->SetParameter(2, 0);
fFitFunction->SetParLimits(2, -10, 10);
if(fPlayer->fUsePol1==false)
fFitFunction->FixParameter(2,0);
for(uint i=0 ; i<fEnergies.size() ; i++) {
//Height
fFitFunction->SetParameter(3+i*6+0, fHistogram->GetBinContent(fHistogram->FindBin(fEnergies[i])) - (fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[0]))+fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[1])))*0.5 );
fFitFunction->SetParLimits(3+i*6+0, fHistogram->GetBinContent((fHistogram->FindBin(fEnergies[i]))-(fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[0]))+fHistogram->GetBinContent(fHistogram->FindBin(fBackgd[1])))*0.5)*0.5, fHistogram->GetBinContent(fHistogram->GetMaximumBin()));
//Position
fFitFunction->SetParameter(3+i*6+1, fEnergies[i]);
fFitFunction->SetParLimits(3+i*6+1, fEnergies[i]-DefFWHM, fEnergies[i]+DefFWHM);
if(fPlayer->fFixMean->GetState() == kButtonDown)
fFitFunction->FixParameter(3+i*6+1,fEnergies[i]);
//FWHM
fFitFunction->SetParameter(3+i*6+2, DefFWHM);
fFitFunction->SetParLimits(3+i*6+2, DefFWHM_min, DefFWHM_max);
if(fPlayer->fFixFWHM->GetState() == kButtonDown)
fFitFunction->FixParameter(3+i*6+2,DefFWHM);
//LeftTail
fFitFunction->SetParameter(3+i*6+3, LeftTailVal);
fFitFunction->SetParLimits(3+i*6+3, LeftTailValMin, LeftTailValMax);
if(fPlayer->fUseLT->GetState() == kButtonUp)
fFitFunction->FixParameter(3+i*6+3,-5);
else if(fPlayer->fFixLT->GetState() == kButtonDown)
fFitFunction->FixParameter(3+i*6+3,LeftTailVal);
//RightTail
fFitFunction->SetParameter(3+i*6+4, RightTailVal);
fFitFunction->SetParLimits(3+i*6+4, RightTailValMin, RightTailValMax);
if(fPlayer->fUseRT->GetState() == kButtonUp)
fFitFunction->FixParameter(3+i*6+4,5);
else if(fPlayer->fFixRT->GetState() == kButtonDown)
fFitFunction->FixParameter(3+i*6+4,RightTailVal);
//AmplitudeStep
fFitFunction->SetParameter(3+i*6+5, StepVal);
fFitFunction->SetParLimits(3+i*6+5, StepValMin, StepValMax);
if(fPlayer->fUseStep==false)
fFitFunction->FixParameter(3+i*6+5,0);
}
fHistogram->GetXaxis()->SetRangeUser(x,y);
TFitResultPtr r = fHistogram->Fit(fFitFunction,"R0S");
if(fPlayer->fFixAmpli->GetState() == kButtonDown) {
//Extract Background
if(fBackFunction) delete fBackFunction;
fBackFunction = new TF1("Background", this, &CXFit::StepedBackground, fBackgd[0], fBackgd[1], NPars, "CXFit", "StepedBackground");
fBackFunction->SetParameters(fFitFunction->GetParameters());
for(uint i=0 ; i<fEnergies.size() ; i++) {
Float_t Ampli = fHistogram->GetBinContent(fHistogram->FindBin(fEnergies[i]))- fBackFunction->Eval(fEnergies[i]);
fFitFunction->FixParameter(3+i*6+0,Ampli);
}
r = fHistogram->Fit(fFitFunction,"R0S");
}
//Extract Background
if(fBackFunction) delete fBackFunction;
fBackFunction = new TF1("Background", this, &CXFit::StepedBackground, fBackgd[0], fBackgd[1], NPars, "CXFit", "StepedBackground");
fBackFunction->SetParameters(fFitFunction->GetParameters());
fBackFunction->SetNpx(1000);
fBackFunction->SetLineColor(kBlue);
fBackFunction->Draw("same");
//Extract Residue
if(fResidue) delete fResidue;
fResidue = new TF1("Residue", this, &CXFit::Residue, fBackgd[0], fBackgd[1], NPars, "CXFit", "Residue");
fResidue->SetParameters(fFitFunction->GetParameters());
fResidue->SetNpx(1000);
fResidue->SetLineWidth(1);
fResidue->SetLineColor(kBlack);
fResidue->Draw("same");
ostringstream text;
text << "Fit results :";
cout<<text.str()<<endl;fPlayer->PrintInListBox(text.str(),kPrint);text.str("");
text << "Chi2 = "<< r->Chi2();
cout<<text.str()<<endl;fPlayer->PrintInListBox(text.str(),kInfo);text.str("");
text << "Ndf = "<< r->Ndf();
cout<<text.str()<<endl;fPlayer->PrintInListBox(text.str(),kInfo);text.str("");
text << "P val = "<< r->Prob();
cout<<text.str()<<endl;fPlayer->PrintInListBox(text.str(),kInfo);text.str("");
fListOfPeaks->Clear();
for(uint i=0 ; i< fEnergies.size() ; i++) {
text<<"Peak "<<i<<":";
cout<<text.str()<<endl;fPlayer->PrintInListBox(text.str(),kPrint);text.str("");
TF1 *peak = new TF1(Form("Peak%d",i), this, &CXFit::PeakFunction, fBackgd[0], fBackgd[1], NPars, "CXFit", "PeakFunction");
peak->SetParameters(fFitFunction->GetParameters());
peak->SetParErrors(fFitFunction->GetParErrors());
peak->SetParameter(1,1);//with backgroud
peak->SetParameter(0,i);
peak->SetNpx(1000);
peak->SetLineColor(kGreen);
peak->SetLineStyle(kDashed);
peak->Draw("same");
fListOfPeaks->Add(peak);
Float_t Area = (peak->Integral(fBackgd[0],fBackgd[1],1e-6)-fBackFunction->Integral(fBackgd[0],fBackgd[1],1e-6))/fHistogram->GetBinWidth(1);
Float_t AreaErr = 2*sqrt(Area);
Float_t Mean = peak->GetParameter(3+i*6+1);
Float_t MeanErr = peak->GetParError(3+i*6+1);
Float_t FWHM = peak->GetParameter(3+i*6+2);
Float_t FWHMErr = peak->GetParError(3+i*6+2);
Float_t LeftT = TMath::Abs(peak->GetParameter(3+i*6+3));
Float_t LeftTErr = peak->GetParError(3+i*6+3);
Float_t Right = peak->GetParameter(3+i*6+4);
Float_t RightErr = peak->GetParError(3+i*6+4);
peak->SetParameter(1,0);//without backgroud
Float_t Max = peak->GetParameter(3+i*6+0);
Float_t MaxErr = peak->GetParError(3+i*6+0);
Float_t FWHM_L = peak->GetX(Max/2,fBackgd[0],Mean,1e-6);
Float_t FWHM_L_err = peak->GetX((Max-MaxErr)/2,fBackgd[0],Mean,1e-6);
Float_t FWHM_R = peak->GetX(Max/2,Mean,fBackgd[1],1e-6);
Float_t FWHM_R_err = peak->GetX((Max-MaxErr)/2,Mean, fBackgd[1],1e-6);
Float_t FWHM_Real = FWHM_R-FWHM_L;
Float_t FWHM_Real_err = (FWHM_R_err-FWHM_L_err)-FWHM_Real;
peak->SetParameter(1,1);//with backgroud
text<<left<<setw(11)<<"Mean"<<": "<<setprecision(7)<<setw(10)<<Mean<<" ("<<setprecision(7)<<setw(10)<<MeanErr<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"Amplitude"<<": "<<setprecision(7)<<setw(10)<<Max<<" ("<<setprecision(7)<<setw(10)<<MaxErr<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"FWHM (gaus)"<<": "<<setprecision(7)<<setw(10)<<FWHM<<" ("<<setprecision(7)<<setw(10)<<FWHMErr<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"FWHM (real)"<<": "<<setprecision(7)<<setw(10)<<FWHM_Real<<" ("<<setprecision(7)<<setw(10)<<FWHM_Real_err<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"L Tail"<<": "<<setprecision(7)<<setw(10)<<LeftT<<" ("<<setprecision(7)<<setw(10)<<LeftTErr<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"R Tail"<<": "<<setprecision(7)<<setw(10)<<Right<<" ("<<setprecision(7)<<setw(10)<<RightErr<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
text<<left<<setw(11)<<"Area"<<": "<<setprecision(7)<<setw(10)<<Area<<" ("<<setprecision(7)<<setw(10)<<AreaErr<<")";
cout<<text.str()<<endl;
fPlayer->PrintInListBox(text.str(),kInfo);
text.str("");
}
fFitFunction->Draw("same");
fPlayer->GetMainWindow()->RefreshPads();
}
double CXFit::DoubleTailedStepedGaussian(double*xx,double*pp)
{