RecalEnergy.h

Go to the documentation of this file.

#ifndef GWRecal_h
#define GWRecal_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 "ReadSpek.h"
#include "FitSpek.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), eref(-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    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 GWRecal
{

public:

    GWRecal();
    ~GWRecal(){}

    void StartCalib();


public:

    void SetFileName(TString FileName){specName = FileName;}

    TString fSpectraFolderName;
    void SetSpectraFolder(TString name){fSpectraFolderName = name;}

    void SetNSpec(int N){specNN = N;} // number of subspectra to analyze [1]
    void SetHistNumber(int HistNbr){specN0 = HistNbr;} //analysis starts from subspectrum [0]
    void SetStep(int Step){specNs = Step;} //analysis proceeds incrementing subspectrum by [1]

    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 SetHistChannelLimits(int SpecNbr, float ChFrom, float ChTo);//limit the search to this range in channels for spectrum nn

    TString fSourceName;
    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)
    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;nwa.clear();} //default fwhm and minmum amplitude for the peaksearch [10 5]
    void SetHistPeaksLimits(int SpecNbr, float FWHM, float Ampli); //fwhm and minmum amplitude for spectrum nn

    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 UnsetLeftTail(){useTL = false;} // disable using the Left  Tail in peak fit
    void UnsetRightTail(){useTR = false;} // disable using the Right Tail in peak fit
    void UseLinearCalib(){numZero++;} // add a fake peak at (0,0) to push calibration through origin
    void UseAffineCalib(){doPoly1 = true;} // 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 PrintResForXTalk(int VerboseMod){modXtalk = VerboseMod;} //print results as when calculating xTalk coefficients mod==36/37/38
    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




public:

    vector<Fitted> Peaks;
    vector<Fitted> GetFitResults(){return Peaks;}

    string  specName;
    string  specFormat;
    int     specLength;
    float  *specData;
    int     specOffset;           // subtracted immediately after peak fit

    ReadSpek Spek;    // reading managed by this

    int     specNN;     // number of
    int     specN0;     // from
    int     specNs;     // step  (todo)

    int     specFromDef, specToDef;
    int     fCurrentspecFromDef, fCurrentspecToDef;
    int    *specFrom;
    int    *specTo;

    float   specFWHMdef;
    float   specAMPLdef;
    float   fCurrentspecFWHMdef;
    float   fCurrentspecAMPLdef;
    float  *specFWHM;
    float  *specAMPL;

    FILE   *testFP;
    float  *testData;
    float  testGain;
    string testFileName;
    string testSuff;     // suffix for the output spectra
    bool   testResults;

    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)
    int    modXtalk;

    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;

    CFitSpek   *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

    int  gterminate;

    // callback for SIGINT
    void terminate(int /*sig*/) {
        gterminate++;
        if(gterminate > 3)
            exit(EXIT_FAILURE);
    }

    // functions of ReadATCA
    void    numparams (char *, int);
    void    initialize();
    bool    TestResults();
    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(int nspe, 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;
};

#endif