add in tools the utility of calibration adapted from RecalEnergy2 of Dino

src/root/tools/CMakeLists.txt

@@ -34,6 +34,8 @@ set( with_dictionnaries
 	LevelEditor
 	MatPlayer
 	GSPlayerTUI
+	GwFitSpek
+	GwRecalEnergy
 )
 set( without_dictionnaries
 	RADConverter

src/root/tools/GwFitSpek.h

+#ifndef GwCFitSpek_h
+#define GwCFitSpek_h
+
+// FitSpek.h: interface for the GwCFitSpek class.
+//
+//////////////////////////////////////////////////////////////////////
+
+#if _MSC_VER > 1000
+#pragma once
+#endif // _MSC_VER > 1000
+
+#include <vector>
+
+// Numero di parametri del fit gaussiano
+#define NUMFITPAR  8
+
+///////////////////////////////////////////////////////////////////////////////////
+// Macro per l'indirizzamento C-style [i1][i2] di matrice dimensionata [..][l2]
+// e allocata come array in un blocco contiguo di memoria. 
+// In C++ si puo' fare in modo piu' elegante, ma cosi' e' facile
+#define M12(i1, i2, l2) (i1*l2+i2)
+
+const int p_Free      = 0;
+const int p_Linked    = 1;
+const int p_Fixed     = 2;
+
+const int p_BNone = 0;
+const int p_BLow  = 1;
+const int p_BHig  = 2;
+const int p_BBoth = 3;
+
+///////////////////////////////////////////////////////////////////////////////////
+// In questa classe si usano ( in Curfit() ) puntatori a funzioni non statiche.
+// La sintassi e' complicata 
+// double  (GwCFitSpek::*pfun)(double); // nella dichiarazione si deve specificare classe e segnatura
+// pfun = ExpFunc;                    // niente da segnalare nella assegnazione
+// double aa = ExpFunc(10.);          // la funzione usata di per se
+// double bb = (*this.*pfun)(10.);    // usata come puntatore, ci si deve accedere attraverso un oggetto !!
+// double bb = (this->*pfun)(10.);    //  ""
+
+const double  m_pi    = 4.*atan(1.);
+const double  m_sqrt2 = sqrt(2.);
+const double  sqrt2pi = sqrt(8.*atan(1.));
+const double  s2fwhm  = sqrt(8.*log(2.));
+const double  s2fwtm  = sqrt(8.*log(10.));
+
+class GwCFitSpek
+{
+public:
+  GwCFitSpek();
+  virtual ~GwCFitSpek();
+
+public:
+	void    Reset();
+  bool    Valid()      {return  m_bValid;}
+  bool    Function()   {return  m_bValid && m_bFunct;}
+  bool    Background() {return  m_bValid && m_bBckgr;}
+
+  double  Chi2()      {return (m_bValid) ? m_dChi2   : -1.;}
+
+  int     From()      {return (m_bValid) ? m_nChanA    : 0;}
+  int     To()        {return (m_bValid) ? m_nChanB    : 0;}
+
+  int     BgFrom()    {return (m_bValid) ? m_nBgChanA  : 0;}
+  int     BgTo()      {return (m_bValid) ? m_nBgChanB  : 0;}
+
+  int     Npeaks()    {return (m_bValid) ? m_nNpeaks   : 0;}
+
+  bool    BackFlat()  {return m_bFlat;}   
+  bool    BackSlope() {return m_bSlope;}   
+
+  double  Back0   ()  {return (m_bBckgr) ? m_dBack0    : 0.;}
+  double  Back0Err()  {return (m_bBckgr) ? m_dBack0Err : 0.;}
+
+  double  Back1   ()  {return (m_bBckgr) ? m_dBack1    : 0.;}
+  double  Back1Err()  {return (m_bBckgr) ? m_dBack1Err : 0.;}
+
+  double  BackL1  ()  {return (m_bBckgr) ? m_dBackL1   : 0.;}
+  double  BackL2  ()  {return (m_bBckgr) ? m_dBackL2   : 0.;}
+
+  double  Area   (int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Area  : 0.;}
+  double  AreaErr(int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pErr[np].Area  : 0.;}
+
+  double  Amplitude   (int np) {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Ampli  : 0.;}
+  double  AmplitudeErr(int np) {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pErr[np].Ampli  : 0.;}
+
+  double  Position   (int np)  {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Posi  : 0.;}
+  double  PositionErr(int np)  {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pErr[np].Posi  : 0.;}
+
+  double  Sigma   (int np) {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Sigma  : 0.;}
+  double  SigmaErr(int np) {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pErr[np].Sigma  : 0.;}
+
+  double  Fwhm   (int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Fwhm  : 0.;}
+  double  FwhmErr(int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pErr[np].Fwhm  : 0.;}
+
+  double  Fw05   (int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Fw05  : 0.;}
+  double  Fw01   (int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].Fw01  : 0.;}
+
+  double  W01L   (int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].W01L  : 0.;}
+  double  W01R   (int np)      {return (m_bGauss && np >=0 && np < m_nNpeaks) ? m_pRes[np].W01R  : 0.;}
+
+  double  Step   (int np)      {return (m_bStep && np >=0 && np < m_nNpeaks) ? m_pRes[np].Step  : 0.;}
+  double  StepErr(int np)      {return (m_bStep && np >=0 && np < m_nNpeaks) ? m_pErr[np].Step  : 0.;}
+
+  double  TailLeft    (int np) {return (m_bTailL && np >=0 && np < m_nNpeaks) ? m_pRes[np].TailL : 0.;}
+  double  TailLeftErr (int np) {return (m_bTailL && np >=0 && np < m_nNpeaks) ? m_pErr[np].TailL : 0.;}
+
+  double  TailRight   (int np) {return (m_bTailR && np >=0 && np < m_nNpeaks) ? m_pRes[np].TailR : 0.;}
+  double  TailRightErr(int np) {return (m_bTailR && np >=0 && np < m_nNpeaks) ? m_pErr[np].TailR : 0.;}
+
+  double  ExpAmpli()           {return (m_bExp) ? m_dExpAmpli : 0.;}
+  double  ExpDecay()           {return (m_bExp) ? m_dExpDecay : 0.;}
+
+  double  SinOmega   ()        {return (m_bSin) ? m_pRes[0].Omega : 0.;}
+  double  SinOmegaErr()        {return (m_bSin) ? m_pErr[0].Omega : 0.;}
+
+  double  SinPhase   ()        {return (m_bSin) ? m_pRes[0].Phase : 0.;}
+  double  SinPhaseErr()        {return (m_bSin) ? m_pErr[0].Phase : 0.;}
+
+  double  SinAmpli   ()        {return (m_bSin) ? m_pRes[0].Ampli : 0.;}
+  double  SinAmpliErr()        {return (m_bSin) ? m_pErr[0].Ampli : 0.;}
+
+  double  FitFunc (double x)         {return (m_bValid) ? (this->*m_pfFitFunc)(x)    : 0;}
+  double  FitFuncS(double x) {
+    if(m_bGauss) {
+      bool tTailL = m_bTailL; bool tTailR = m_bTailR;
+      m_bTailL = m_bTailR = false;
+      double res = (m_bValid) ? (this->*m_pfFitPeak)(x,0) + FitBack(x) : 0;
+      m_bTailL = tTailL; m_bTailR = tTailR;
+      return res;
+    }
+    else {
+      return (m_bValid) ? (this->*m_pfFitFunc)(x)    : 0;
+    }
+  }
+  double  FitBack(double x)          {return (m_bBckgr) ? (this->*m_pfFitBack)(x)    : 0;}
+  double  FitPeak(double x, int n=0) {return (m_bFunct) ? (this->*m_pfFitPeak)(x, n) : 0;}
+
+  bool    CalcStrightLine(float *pSpek, int chanA, int chanB, int doSlope);
+  bool    CalcPeakMoments(float *pSpek, int chanA, int chanB, int nChanBack);
+  int     CalcGaussianFit(float *pSpek, int chanA, int chanB, std::vector<double>&vPeaks);
+  int     CalcGaussianFit(float *pSpek, int chanA, int chanB);  // one peak autolocated,
+  bool    CalcExponential(float *pSpek, int chanA, int chanB, int bgA = -1, int bgB = -1);
+  int     CalcSinusoid   (float *pSpek, int chanA, int chanB, double period, bool bFlat, bool bSlope);
+
+public:
+  int     GFitParNumber() {return NUMFITPAR;}
+  int     GFitParB0 ()     {return 0;}
+  int     GFitParB1 ()     {return 1;}
+  int     GFitParAMP()     {return 2;}
+  int     GFitParPOS()     {return 3;}
+  int     GFitParSIG()     {return 4;}
+  int     GFitParSTP()     {return 5;}
+  int     GFitParTL ()     {return 6;}
+  int     GFitParTR ()     {return 7;}
+  void    GFitParGet(int *pP) {for(int ii = 0; ii < NUMFITPAR; ii++) pP[ii] = m_pGFitPar[ii];}
+  void    GFitParSet(int *pP) {for(int ii = 0; ii < NUMFITPAR; ii++) m_pGFitPar[ii] = pP[ii] & 3;}
+  void    GFitParDlg();
+
+public:
+  
+  class CFitRes
+  {
+  public:
+    double  Ampli;
+    double  Area;
+    double  Posi;
+    double  Sigma;
+    double  Fwhm;   // simply Sigma*s2fwhm
+    double  Fw05;   // FW at 1/2  considering tails
+    double  Fw01;   // FW at 1/10 considering tails
+    double  W01L;   // left  width at 1/10 normalized to Fwhm/2
+    double  W01R;   // right width at 1/10 normalized to Fwhm/2
+    double  Step;
+    double  TailL;
+    double  TailR;
+    double  Omega;
+    double  Phase;
+    CFitRes() : Ampli(0), Area(0),  Posi(0),  Sigma(0), Fwhm(0), Fw05(0), Fw01(0), W01L(0), W01R(0),
+                Step(0),  TailL(0), TailR(0), Omega(0), Phase(0)
+    {
+    }
+  };
+  
+  class CFitPar
+  {
+  public:
+    double  Value;
+    double  Error;
+    double  Valmin;
+    double  Valmax;
+    int     Status;
+    int     Bounds;
+    int     LinkTo;
+    int     PackId;
+    CFitPar() : Value(0),  Error(0),  Valmin(0), Valmax(0),
+                Status(0), Bounds(0), LinkTo(0), PackId(0)
+    {
+    }
+  };
+  
+  // i membri puntatori alle funzioni del risultato del fit
+  double (GwCFitSpek::*m_pfFitFunc)(double);
+  double (GwCFitSpek::*m_pfFitBack)(double);
+  double (GwCFitSpek::*m_pfFitPeak)(double, int);
+  
+  // i membri puntatori alle funzioni del fit per Curfit()
+  double (GwCFitSpek::*m_pfFunct)(double, double *);
+  void   (GwCFitSpek::*m_pfDeriv)(int, double *, double *, double &, double &);
+
+private:
+  double  LineFunc(double xx);
+
+  double  ExpFunc(double xx);
+  double  ExpPeak(double xx, int);
+
+  double  GFitFunc(double xx);
+  double  GFitBack(double xx);
+  double  GFitPeak(double xx, int nn);
+  double  GFfunF(double xx, double *par);
+  void    GFderF(int id, double *fpar, double *deriv, double &deltay, double &weight);
+  void    GF_res(void);
+  double  GF_erf(double);
+  double  GF_y2w(int np, double yval, int side = 0);  // -1=Left, 0=full, 1=Right
+
+  double  SFitFunc(double xx);
+  double  SFitWave(double xx, int nn);
+  double  SFfunF(double xx, double *par);
+  void    SFderF(int id, double *fpar, double *deriv, double &deltay, double &weight);
+
+  double  Curfit(int ndat, int npar, int mpar);
+  double  Chisqr(double *par);
+  void    Manage(double *par, int type);
+  double  Matinv(double *array, int nord, int mpar);
+
+  bool    m_bValid;
+  bool    m_bFunct;
+  bool    m_bBckgr;
+
+  bool    m_bLine;
+  bool    m_bFlat;
+  bool    m_bSlope;
+
+  bool    m_bGauss;
+  bool    m_bStep;
+  bool    m_bTailL;
+  bool    m_bTailR;
+
+  int     m_nChanA;
+  int     m_nChanB;
+  int     m_nBgChanA;
+  int     m_nBgChanB;
+  int     m_nNpeaks;
+
+  // control of Gauss Fit parameters
+  // [Bit1 enabled] [Bit2 all equal]
+  int     m_pGFitPar[NUMFITPAR];
+
+  int     m_nNdat;
+  int     m_nNpar;
+  int     m_nFpar;
+  int     m_nOffB;
+  int     m_nOffN;
+
+  double  m_dChi2;
+
+  double  m_dBack0;
+  double  m_dBack0Err;
+  double  m_dBack1;
+  double  m_dBack1Err;
+  double  m_dBackL1;
+  double  m_dBackL2;
+
+  CFitRes *m_pRes;
+  CFitRes *m_pErr;
+
+  bool    m_bExp;
+  double  m_dExpAmpli;
+  double  m_dExpDecay;
+
+  bool    m_bSin;
+
+  CFitPar *m_pPar;        // the parameters of the fit
+  double  *m_pDatVal;     // the data to fit
+  double  *m_pDatErr;     // their error (1/sig2)
+
+};
+
+#endif // !defined(AFX_FITSPEK_H__F2A2E380_9934_11D5_B3CF_000000000000__INCLUDED_)

src/root/tools/GwRecalEnergy.h

+#ifndef GwRecalEnergy_h
+#define GwRecalEnergy_h
+
+#include <ctime>
+#include <csignal>
+#include <fstream>
+#include <sstream>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <stdlib.h>
+#include <cmath>
+#include <memory.h>
+
+#include <vector>
+#include <algorithm>
+#include <functional>
+
+#include "TString.h"
+#include "GwFitSpek.h"
+
+using namespace std;
+
+struct Fitted
+{
+    Fitted() : NSubPeaks(0), BgdOff(0), BgdSlope(0), BgFrom(0), BgTo(0), area(0), ampli(0), posi(0), fw05(0), fw01(0), fwhm(0), tailL(0), tailR(0), Lambda(0), Rho(0), S(0), erefindex(-1), good(false) {;}
+    int NSubPeaks;
+    double BgdOff;
+    double BgdSlope;
+    double BgFrom;
+    double BgTo;
+    double area;
+    double ampli;
+    double posi;
+    double fw05;
+    double fw01;
+    double fwhm;
+    double tailL;
+    double tailR;
+    double Lambda;
+    double Rho;
+    double S;
+    int    erefindex; // the corresponding calibration line, if good==true
+    double eref; // the corresponding calibration line, if good==true
+    bool   good;
+};
+
+struct NWA_t
+{
+    NWA_t() : index(-1), width(0), ampli(0) {;}
+    int   index;
+    float width;
+    float ampli;
+};
+
+struct NLimits_t
+{
+    NLimits_t() : index(-1), from(0), to(0) {;}
+    int   index;
+    int from;
+    int to;
+};
+
+struct smallerPosi
+{
+    bool operator()( Fitted elem1, Fitted elem2 ) {
+        return elem1.posi < elem2.posi;
+    }
+};
+struct largerArea
+{
+    bool operator()( Fitted elem1, Fitted elem2 ) {
+        return elem1.area > elem2.area;
+    }
+};
+struct largerAmplitude
+{
+    bool operator()( Fitted elem1, Fitted elem2 ) {
+        return elem1.ampli > elem2.ampli;
+    }
+};
+
+class TH1;
+
+class GwRecalEnergy
+{
+
+public:
+
+    GwRecalEnergy();
+    virtual ~GwRecalEnergy(){}
+
+    void StartCalib();
+
+    TString fSpectraFolderName;
+    TString fSourceName;
+
+    Int_t fId;
+
+public:
+
+    void SetDataFromHistTH1(TH1 *hist, Int_t Id=0);
+    void SetFileName(TString FileName){specName = FileName;}
+
+    void SetChannelOffset(int Off){specOffset = Off;} //channel offset to subtract to the position of the peaks [0]
+    void SetGlobalChannelLimits(int ChFrom, int ChTo){specFromDef = ChFrom; specToDef = ChTo; nlim.clear();} //limit the search to this range in channels
+
+    void SetSource(TString SourceName){fSourceName = SourceName; AnalyseSources();}
+    void AnalyseSources();
+    void AddPeak(double EPeak){Energies.push_back(EPeak);} //add this energy to the list of lines (can be given more than once)
+    Int_t GetNEnergies(){return Energies.size();}
+    vector< double > GetEnergies(){return Energies;}
+    void RemoveClosestPeakTo(double EPeak){Delendae.push_back(EPeak);} //remove the line closest to this energy from the list of lines
+    void SetRefPeak(double EPeak){refEner = EPeak;} //energy (keV) of the reference peak for extended printouts
+    void SetGlobalPeaksLimits(float DefFWHM, float DefAmpli){specFWHMdef = DefFWHM; specAMPLdef = DefAmpli;} //default fwhm and minmum amplitude for the peaksearch [10 5]
+
+    void UseFlatBackGround(){fFlatBg = true; fAffineBg = false;} // Use a flat function background to fit the peaks
+    void UseAffineBackGround(){fAffineBg = true; fFlatBg = false;} // Use a affine function background to fit the peaks
+    void NoBackGround(){fAffineBg = false; fFlatBg = false;}    // No background estimation to fit the peaks
+
+    void UseFirstDerivativeSearch(){Dmode = 1;} // use the 1st-derivative search (default, always for 2-line sources)
+    void UseSecondDerivativeSearch(){Dmode = 2;} // use the 2nd-derivative search
+
+    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 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]
+    void SetFullVerbosity(){Verbosity = 0xFFFF;} //Set full verbosity
+    void Reset();
+
+    Float_t GetOffset(){return offset1;}
+    Float_t GetSlope(){return slope1*hGain;}
+
+public:
+
+    vector<Fitted> Peaks;
+    vector<Fitted> GetFitResults(){return Peaks;}
+
+    string  specName;
+    int     specLength;
+    float  *specData = nullptr;
+    int     specOffset;           // subtracted immediately after peak fit
+
+    int     specFromDef, specToDef;
+
+    float   specFWHMdef;
+    float   specAMPLdef;
+
+    vector<double> specPeaks;
+    vector<double> Energies;
+    vector<double> Delendae;
+
+    double eBhead;      // band head
+    double eBstep;      // delta
+    int    eBnumb;  // numbero of peaks
+
+    double refEner;
+    bool  useTL;
+    bool  useTR;
+
+    bool fFlatBg;
+    bool fAffineBg;
+
+    vector<NWA_t> nwa;
+    vector<NLimits_t> nlim;
+
+    bool   bTwoLines;
+    bool   bOneLine;
+    bool   doSlope;
+    bool   doPoly1;   // linear
+    bool   doPoly2;   // cubic
+    int    numZero;       // number of fake peaks at (0,0)
+
+    bool   useErrors ;   // fit using position and energy errors
+    double errE;
+
+    int    Dmode;   // peak search using first derivative
+    int    xP_0, xP_1, xP_2, xP_3, xP_4, xP_5, xP_6;
+    float  xP_minAmpl;
+    float  xP_fThreshCFA;
+    int    xP_nMinWidthP;
+    int    xP_nMinWidthN;
+    int    xP_nMinWidthP1;
+    int    xP_nMinWidthP2;
+
+    GwCFitSpek   *m_pCFit;
+
+    double slope, tshift;
+    double offset1, slope1;               // should be generalized to polynomials
+    double offset2, slope2, curv2;
+
+    float hGain;
+
+    int  Verbosity;  // &1 PeakFit  &2 CalibPeaks   &4 CalibSort
+
+    // functions of ReadATCA
+    void    initialize();
+    int     PeakSearch1(float * data, int chA, int chB, float fwhm, float minAmpl, std::vector<double>&vPeaks);
+    int     PeakSearch2(float * data, int chA, int chB, float fwhm, float minAmpl, std::vector<double>&vPeaks);
+    int     LargePeaks1(float * data, int chA, int chB, float fwhm, float minAmpl, std::vector<double>&vPeaks, int maxNum);
+    void    SmoothSpek(float *spek, int nchan, double sigma, int ndiff, float *& tSpek, int &start);
+    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)
+    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);
+
+    clock_t startTime, stopTime;
+
+    ClassDef(GwRecalEnergy,0);
+};
+
+#endif

src/root/tools/LinkDef.h

@@ -44,6 +44,6 @@
 #pragma link C++ class Gw::LevelEditor;
 #pragma link C++ class Gw::MatPlayer;
 #pragma link C++ class Gw::GSPlayerTUI+;
-
+#pragma link C++ class GwRecalEnergy;
 
 #endif