RecalEnergy.cpp

Go to the documentation of this file.

//
// cpp version of recal_cob
// improved search of the peaks
// using user-defined sources
//
// Dino Bazzacco, August 2014
//

#include "TSystem.h"

#include "RecalEnergy.h"


GWRecal::GWRecal() :
    fSpectraFolderName(""),
    fSourceName("60Co"),
    Peaks(),
    specName(""),
    specFormat(""),
    specLength(0),
    specData(NULL),
    specOffset(0),
    Spek(),
    specNN(1),
    specN0(0),
    specNs(1),
    specFromDef(0),
    specToDef(0),
    specFrom(NULL),
    specTo(NULL),
    specFWHMdef(15),
    specAMPLdef(5),
    fCurrentspecFWHMdef(0),
    fCurrentspecAMPLdef(0),
    specFWHM(NULL),
    specAMPL(NULL),
    testFP(NULL),
    testData(NULL),
    testGain(1),
    testFileName(""),
    testSuff(""),
    testResults(false),
    specPeaks(),
    Energies(),
    Delendae(),
    eBhead(0.0),
    eBstep(0.0),
    eBnumb(0),
    refEner(0),
    useTL(true),
    useTR(true),
    fFlatBg(1),
    fAffineBg(0),
    nwa(),
    nlim(),
    bTwoLines(false),
    bOneLine(false),
    doSlope(true),
    doPoly1(false),
    doPoly2(false),
    numZero(0),
    modXtalk(0),
    useErrors(false),
    errE(0.01),
    Dmode(1),
    xP_0(0), xP_1(0), xP_2(0), xP_3(0), xP_4(0), xP_5(0), xP_6(0),
    xP_minAmpl(0.0),
    xP_fThreshCFA(0.0),
    xP_nMinWidthP(0),
    xP_nMinWidthN(0),
    xP_nMinWidthP1(0),
    xP_nMinWidthP2(0),
    m_pCFit(NULL),
    slope(0),
    tshift(0),
    offset1(0),
    slope1(0),
    offset2(0),
    slope2(0),
    curv2(0),
    hGain(1),
    Verbosity(0),
    gterminate(0)
{
    AnalyseSources();
}


void GWRecal::StartCalib()
{
    initialize();                 // initialize parameters related to digitisers and detectors

    //gterminate = 0;
    //signal(SIGINT, terminate);
    //cout << "\nAnalysis can be stopped by typing CTRL_C\n" << endl << flush;

//    int  evnum = 0;
//    bool evok  = true;
    startTime  = clock();

    for(int nspe = 0; nspe < specNN; nspe++) {

        int ispec = specN0 + nspe*specNs;
        if( !Spek.Read(specLength, ispec) )
            break;
        specData = Spek.Data();

        fCurrentspecFWHMdef = specFWHM[nspe];
        fCurrentspecAMPLdef = specAMPL[nspe];

        fCurrentspecFromDef = specFrom[nspe];
        fCurrentspecToDef   = specTo[nspe];

        int np    = 0;
        int ngood = 0;

        if(bOneLine)
            np = LargePeaks1(specData, fCurrentspecFromDef, fCurrentspecToDef, fCurrentspecFWHMdef, fCurrentspecAMPLdef, specPeaks, 1);
        else if(bTwoLines)
            np = LargePeaks1(specData, fCurrentspecFromDef, fCurrentspecToDef, fCurrentspecFWHMdef, fCurrentspecAMPLdef, specPeaks, 6);
        else if(Dmode == 1)
            np = PeakSearch1(specData, fCurrentspecFromDef, fCurrentspecToDef, fCurrentspecFWHMdef, fCurrentspecAMPLdef, specPeaks);
        else
            np = PeakSearch2(specData, fCurrentspecFromDef, fCurrentspecToDef, fCurrentspecFWHMdef, fCurrentspecAMPLdef, specPeaks);

        //cout << specName << " # " << setw(2) /*<< setfill('0')*/ << nspe /*<< setfill(' ')*/ << setw(5) <<  np << "   ";
        //if(Verbosity) cout << endl;
        //cout<<Form("%5d  %s # %5d  ( %3d", nspe, specName.c_str(), ispec, np);
        //cout<<Form("%5d # %5d  %3d", nspe, ispec, np);

        if(Verbosity & 1) {
            cout<<Form("#1 %5d %6d %4d", nspe, ispec, np);
        }

        slope=0, tshift = 0;
        offset1=0, slope1=0;               // should be generalized to polynomials
        offset2=0, slope2=0, curv2=0;

        if(np) {
            int fp = FitPeaks( Verbosity );
            if(fp) {
                if(specOffset) {               // correct peak positions for specOffs
                    for(size_t np_ = 0; np_ < Peaks.size(); np_++) {
                        Peaks[np_].posi -= specOffset;
                    }
                }
                if(doSlope){
                    slope = ECalibration(nspe, offset1, slope1, offset2, slope2, curv2, Verbosity );
                    if(slope) {
                        for(size_t np_ = 0; np_ < Peaks.size(); np_++) {
                            if( Peaks[np_].good )
                                ngood++;
                        }
                    }
                }
                else {
                    tshift = refEner - Peaks[0].posi;
                    ngood  = 1;
                }
            }
        }
        else {
            Peaks.clear();
            if(Verbosity & 1)
                cout<<Form("\n");
        }

        if(doSlope) {
            double refPeakPosi=0, refPeakEner=0;
            int iref = -1;;
            if(refEner > 0 && slope) {
                refPeakPosi = refEner/slope;                            // position of the reference peak derived from calibration
                for(size_t irp = 0; irp < Peaks.size(); irp++) {
                    if( (abs(Peaks[irp].posi-refPeakPosi) < fCurrentspecFWHMdef)/* && Peaks[irp].good*/ ) { // not restricted to the good peaks
                        iref = irp;
                        refPeakPosi = Peaks[irp].posi;                      // actual peak position
                        refPeakEner = slope*refPeakPosi;
                    }
                }
            }
            if(modXtalk==0)
                cout<<Form("%6d %6d %4d %3d", nspe, ispec, np, ngood);
            else
                cout<<Form("%6d %6d %4d %3d", nspe%modXtalk, ispec/modXtalk, np, ngood);
            double chi2 = 0;
            int nused = 0;
            for(size_t np_ = 0; np_ < Peaks.size(); np_++) {
                if(Peaks[np_].good) {
                    nused++;
                    double ee = slope*(Peaks[np_].posi);
                    chi2 += pow(Energies[Peaks[np_].eref]-ee, 2);
                }
            }
            chi2 = (nused-1) ? chi2/(nused-1) : 0;
            if(iref>=0) {
                if(modXtalk==0) {
                    cout<<Form(" %9.2f %8.3f %8.3f %9.0f %9.2f %6.1f %7.0f %7.3f %7.3f",
                               refPeakEner, Peaks[iref].fw05*slope, Peaks[iref].fw01*slope,
                               Peaks[iref].area, Peaks[iref].posi, Peaks[iref].fwhm, Peaks[iref].ampli, Peaks[iref].tailL, Peaks[iref].tailR);
                }
                else {
                    cout<<Form(" %9.2f %8.3f %8.3f %9.0f %9.2f %6.1f %7.0f %7.3f %7.3f",
                               refPeakEner, Peaks[iref].fw05, Peaks[iref].fw01,
                               Peaks[iref].area, Peaks[iref].posi, Peaks[iref].fwhm, Peaks[iref].ampli, Peaks[iref].tailL, Peaks[iref].tailR);
                }
            }
            else {
                cout<<Form(" %9.2f %8.3f %8.3f %9.0f %9.2f %6.1f %7.0f %7.3f %7.3f", 0., 0., 0., 0., 0., 0., 0., 0., 0.);
            }
            if(modXtalk==0)
                cout<<Form("  %10.6f %6.2f", slope*hGain, min(999.99, chi2*100));
            else {
                if(slope)
                    cout<<Form("  %12.7f", 1./(slope*hGain));
                else
                    cout<<Form("  %12.7f", 0.);
            }
            if(doPoly1) {
                chi2 = 0;
                nused = 0;
                for(size_t np_ = 0; np_ < Peaks.size(); np_++) {
                    if(Peaks[np_].good) {
                        nused++;
                        double ee = offset1 + slope1*(Peaks[np_].posi);
                        chi2 += pow(Energies[Peaks[np_].eref]-ee, 2);
                    }
                }
                chi2 = (nused-2)>0 ? chi2/(nused-2) : 0;
                cout<<Form(" %9.3f %10.6f %6.2f", offset1, slope1*hGain, min(999.99, chi2*100));
            }
            if(doPoly1 && doPoly2) {
                chi2 = 0;
                nused = 0;
                for(size_t np_ = 0; np_ < Peaks.size(); np_++) {
                    if(Peaks[np_].good) {
                        nused++;
                        double ee = offset2 + slope2*Peaks[np_].posi + curv2*Peaks[np_].posi*Peaks[np_].posi;
                        chi2 += pow(Energies[Peaks[np_].eref]-ee, 2);
                    }
                }
                chi2 = (nused-3)>0 ? chi2/(nused-3) : 0;
                cout<<Form(" %9.3f %10.6f % 10.4e %6.2f", offset2, slope2*hGain, curv2*hGain*hGain, min(999.99, chi2*100));
            }
        }
        else {
            cout<<Form("%6d %6d %4d %3d", nspe, ispec, np, ngood);
            if(Peaks.size())
                cout<<Form(" %9.2f %8.3f %8.3f %9.0f %9.2f %6.1f %7.0f %7.3f %7.3f %11.3f",
                           Peaks[0].posi, Peaks[0].fw05, Peaks[0].fw01,
                        Peaks[0].area, Peaks[0].posi, Peaks[0].fwhm, Peaks[0].ampli, Peaks[0].tailL, Peaks[0].tailR, tshift*hGain);
            else
                cout<<Form(" %9.2f %8.3f %8.3f %9.0f %9.2f %6.1f %7.0f %7.3f %7.3f %11.3f", 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
        }
        cout<<Form("\n");

        if(testResults) {
            if(!TestResults()) {
                cout << "Error writing " << testFileName << endl;
                break;
            }
        }
    }

    // TODO : write the calibration files in the needed formats

    stopTime = clock();
//    double realTime = (double)(stopTime - startTime) / CLOCKS_PER_SEC;

    Spek.Close();
}


void GWRecal::initialize()
{
    m_pCFit   = new CFitSpek();
    int  npars = m_pCFit->GFitParNumber();
    int* vpars = (int *) new int [npars];
    m_pCFit->GFitParGet(vpars);

    if(specFormat.empty()) {
//        int spSize = 0;
        if( !SpecNameDecode(specName, specFormat, specLength) ) {
            specFormat = "UI";    // hopefully
            specLength = 16384;
        }
    }

    if( !Spek.Open(specName, specFormat) ) {
        cerr << "Error opening " << specName << endl;
    }

    // exclude the very first and last part of the spectrum
    if(specFromDef == specToDef) {
        specFromDef = int(5*specFWHMdef);
        specToDef   = specLength-1 - specFromDef;
    }

    specFrom = new int[specNN];
    specTo   = new int[specNN];
    for(int nn = 0; nn < specNN; nn++) {
        specFrom[nn] = specFromDef;
        specTo[nn]   = specToDef;
    }
    for(size_t nn = 0; nn < nlim.size(); nn++) {
        int ii = nlim[nn].index;
        if(specNN==1 && specN0==nlim[nn].index)
        {
            specFrom[0] = nlim[nn].from;
            specTo[0]   = nlim[nn].to;
        }
        else if(ii >= 0 && ii < specNN) {
            specFrom[ii] = nlim[nn].from;
            specTo[ii]   = nlim[nn].to;
        }
    }

    cout << endl;
    cout << "# ***********************************"<<endl;
    cout << "# *****                         *****"<<endl;
    cout << "# ***       GWRecal Fit results   ***"<<endl;
    cout << "# *****                         *****"<<endl;
    cout << "# ***********************************"<<endl;
    cout << "# " << endl;

    cout << "# Data format " << specFormat << " " << specLength <<endl;

    cout<< "# Global channel range : [ "<<specFromDef<<" ; "<<specToDef<<" ]"<<endl;
    if(nlim.size()) cout<< "#     Specific channel ranges:"<<endl;
    for(size_t i=0 ; i<nlim.size() ; i++)
    {
        NLimits_t lim = nlim[i];
        int LibNumber = lim.index/38;
        int HistNbr = lim.index%38;
        cout<< "#       ==> Seg nbr "<<Form("%2.0d in library %d",HistNbr,LibNumber)<<": [ "<<lim.from<<" ; "<<lim.to<<" ]"<<endl;
    }
    cout << "# " << endl;

    cout<< "# Global Search peak limits : FWHM < "<< specFWHMdef<<" ; Amplitude > "<<specAMPLdef<<endl;
    if(nwa.size()) cout<< "#     Specific peak limits:"<<endl;
    for(size_t i=0 ; i<nwa.size() ; i++)
    {
        NWA_t tnwa = nwa[i];
        int LibNumber = tnwa.index/38;
        int HistNbr = tnwa.index%38;
        cout<< "#       ==> Seg nbr "<<Form("%2.0d in library %d",HistNbr,LibNumber)<<": FWHM < "<<tnwa.width<<" ; Amplitude > "<<tnwa.ampli<<endl;
    }
    cout << "# " << endl;

    if(fFlatBg)
    {
        vpars[0] = 1;
        vpars[1] = 0;
        fAffineBg = false;
        cout << "# Flat Background estimation " << endl;
    }
    else if(fAffineBg)
    {
        vpars[0] = 1;
        vpars[1] = 1;
        fFlatBg = false;
        cout << "# Affine Background estimation " << endl;
    }
    else if(!fAffineBg && ! fFlatBg)
    {
        vpars[0] = 0;
        vpars[1] = 0;
        fFlatBg = false;
        fAffineBg = false;
        cout << "# No Background estimation " << endl;
    }
    cout << "#  "<< endl;

    if(Dmode==1) cout << "# Using 1st-derivative search"<<endl;
    else if(Dmode==2) cout << "# Using 2nd-derivative search"<<endl;
    cout << "# " << endl;

    if(!useTL) cout << "# Left Tail in peak fit disabled"<<endl<< "# " << endl;
    if(!useTR) cout << "# Right Tail in peak fit disabled"<<endl<< "# " << endl;

    if(numZero) cout << "# Linear calibration enabled (no offset)"<<endl;
    else if(doPoly1 && doPoly2)cout << "#  Parabolic calibration enabled (y = offset + slope*x + curv*x*x)"<<endl;
    else if(doPoly1) cout << "#  Affine calibration enabled (y = offset + slope*x)"<<endl;
    cout<< "# Scaling factor for the slope = "<<hGain<<endl;
    cout<< "# " << endl;

    cout << "#  Selected source : " << fSourceName<<endl;
    cout << "# " << endl;

    // enabling the tails should be decided from the command line
    useTL ? vpars[6] |= 1 : vpars[6] &= ~1;   // control of the Tail to the Left
    useTR ? vpars[7] |= 1 : vpars[7] &= ~1;   // control of the Tail to the Right
    if(!doSlope)
        vpars[5] = 0;  // no step for time shifts
    m_pCFit->GFitParSet(vpars);

    specFWHM = new float[specNN];
    specAMPL = new float[specNN];
    for(int nn = 0; nn < specNN; nn++) {
        specFWHM[nn] = specFWHMdef;
        specAMPL[nn] = specAMPLdef;
    }
    for(size_t nn = 0; nn < nwa.size(); nn++) {
        int ii = nwa[nn].index;
        if(specNN==1 && specN0==nwa[nn].index)
        {
            specFWHM[0] = nwa[nn].width;
            specAMPL[0] = nwa[nn].ampli;
        }
        else if(ii >= 0 && ii < specNN) {
            specFWHM[ii] = nwa[nn].width;
            specAMPL[ii] = nwa[nn].ampli;
        }
    }

    for(int nn = 0; nn < eBnumb; nn++)
        Energies.push_back(eBhead + eBstep*nn);

    if(Energies.size() < 1) {
        Energies.push_back(1173.238);
        Energies.push_back(1332.513);
    }
    if(Delendae.size()) {
        for(size_t dd = 0; dd < Delendae.size(); dd++) {
            size_t closest = 1;
            double delta   = fabs(Delendae[dd]-Energies[0]);
            for (size_t nn = 1; nn < Energies.size(); nn++) {
                if(fabs(Delendae[dd]-Energies[nn]) < delta) {
                    closest = nn;
                    delta   = fabs(Delendae[dd]-Energies[nn]);
                }
            }
            Energies.erase(Energies.begin()+closest);
        }
    }
    sort( Energies.begin( ), Energies.end( ) );

    if(Energies.size()==1) {
        bOneLine  = true;
        bTwoLines = false;
        refEner = Energies[0];
    }
    else if(Energies.size()==2) {
        bOneLine  = false;
        bTwoLines = true;
        if(!refEner) {
            refEner = Energies[1];
        }
        else {
            if(abs(refEner-Energies[0])>1. && abs(refEner-Energies[1])>1.)
                refEner = Energies[1];
        }
    }
    else {
        bOneLine  = false;
        bTwoLines = false;
        if(!refEner)
            refEner = Energies[Energies.size()-1];
    }

    if(doSlope)
        cout << "#  " << "Energies used for the calibration " << endl;
    else
        cout << "#  " << "Position used for the calibration " << endl;
    for (size_t nn = 0; nn < Energies.size(); nn++)
    {
        cout << "#  " << setw(3) << nn << fixed << setprecision(3) << setw(10) << Energies[nn];
        if(refEner==Energies[nn]) cout<<" ==> Reference peak for extended printouts";
        cout<<endl;
    }

    if(testResults) {
        gSystem->mkdir(fSpectraFolderName);
        ostringstream name;
        name << fSpectraFolderName << "/Recal__" << specNN << "-" << specLength << "-F__" << testSuff << ".dat" ;
        testFileName = name.str();
        testFP = fopen(testFileName.c_str(), "wb");
        if(!testFP) {
            cerr << "Error opening " << testFileName << endl;
        }
        testData = new float[specLength];
        memset(testData, 0, sizeof(float)*specLength);
    }

    cout << "#  " << endl;
    cout << "#  " << "Spectra taken from file : " << specName << endl;
    if(testResults)
        cout << "#  Calibrated spectra written to : " << testFileName << endl;


    cout << "#  " << endl;
    if(doSlope) {
        if(modXtalk==0) {
            cout << "# indx  #spec #pks #ok   rEnergy     FW05     FW01      Area  Position  Width   Ampli    WTML    WTMR  slope*gain rChi2%";
            if(doPoly1)
                cout << "   offs1*g   slope1*g rChi2%";
            if(doPoly2)
                cout << "   offs2*g   slope2*g    coeff2*g  rChi2%";
        }
        else {
            cout << "# indx  #spec #pks #ok   rEnergy     FW05     FW01      Area  Position  Width   Ampli    WTML    WTMR    cross-talk";
        }

    }
    else {
        cout << "# indx  #spec #pks #ok   rEnergy     FW05     FW01      Area  Position  Width   Ampli    WTML    WTMR  shift*gain";
    }
    cout << "\n# " << endl;

}

void GWRecal::AnalyseSources()
{
    Energies.clear();

    if(fSourceName.Contains("22Na",TString::kIgnoreCase))
    {
        fSourceName = "22Na";
        Energies.push_back( 511.006);
        Energies.push_back(1274.545);
    }
    if(fSourceName.Contains("56Co",TString::kIgnoreCase))
    {
        fSourceName = "56Co";
        Energies.push_back( 846.772);
        Energies.push_back(1037.840);
        Energies.push_back(1175.102);
        Energies.push_back(1238.282);
        Energies.push_back(1360.250);
        Energies.push_back(1771.351);
        Energies.push_back(2015.181);
        Energies.push_back(2034.755);
        Energies.push_back(2598.458);
        Energies.push_back(3201.962);
        Energies.push_back(3253.416);
        Energies.push_back(3272.990);
        Energies.push_back(3451.152);
        Energies.push_back(3547.925);
    }
    if(fSourceName.Contains("40K",TString::kIgnoreCase))
    {
        fSourceName = "40K";
        Energies.push_back( 122.0614);
        Energies.push_back( 136.4743);
    }
    if(fSourceName.Contains("57Co",TString::kIgnoreCase))
    {
        fSourceName = "57Co";
        //Energies.push_back(  14.4130);
        Energies.push_back( 122.0614);
        Energies.push_back( 136.4743);
    }
    if(fSourceName.Contains("60Co",TString::kIgnoreCase))
    {
        fSourceName = "60Co";
        Energies.push_back(1173.238);
        Energies.push_back(1332.513);
    }
    if(fSourceName.Contains("88Y",TString::kIgnoreCase))
    {
        fSourceName = "88Y";
        Energies.push_back( 898.045);
        Energies.push_back(1836.062);
        //Energies.push_back(2734.087);
    }
    if(fSourceName.Contains("133Ba",TString::kIgnoreCase))
    {
        fSourceName = "133Ba";
        Energies.push_back(  53.156);
        Energies.push_back(  79.623);
        Energies.push_back(  80.999);
        Energies.push_back( 160.609);
        Energies.push_back( 223.116);
        Energies.push_back( 276.404);
        Energies.push_back( 302.858);
        Energies.push_back( 356.014);
        Energies.push_back( 383.859);
    }
    if(fSourceName.Contains("134Cs",TString::kIgnoreCase))
    {
        fSourceName = "134Cs";
        Energies.push_back( 475.36);
        Energies.push_back( 563.27);
        Energies.push_back( 569.30);
        Energies.push_back( 604.68);
        Energies.push_back( 795.78);
        Energies.push_back( 801.86);
        Energies.push_back(1038.53);
        Energies.push_back(1167.89);
        Energies.push_back(1365.17);
    }
    if(fSourceName.Contains("137Cs",TString::kIgnoreCase))
    {
        fSourceName = "137Cs";
        Energies.push_back( 661.661);
    }
    if(fSourceName.Contains("152Eu",TString::kIgnoreCase))
    {
        fSourceName = "152Eu";
        Energies.push_back( 121.7793);
        Energies.push_back( 244.6927);
        Energies.push_back( 344.2724);
        Energies.push_back( 411.111);
        Energies.push_back( 443.979);
        Energies.push_back( 778.890);
        Energies.push_back( 964.014);
        Energies.push_back(1085.793);
        //Energies.push_back(1089.700);
        Energies.push_back(1112.070);
        Energies.push_back(1407.993);
    }
    if(fSourceName.Contains("208Pb",TString::kIgnoreCase))
    {
        fSourceName = "208Pb";
        Energies.push_back(2614.5);   // increase precision
    }
    if(fSourceName.Contains("226Ra",TString::kIgnoreCase))
    {
        fSourceName = "226Ra";
        Energies.push_back( 186.211);
        Energies.push_back( 241.981);
        Energies.push_back( 295.213);
        Energies.push_back( 351.921);
        Energies.push_back( 609.312);
        Energies.push_back( 768.356);
        Energies.push_back( 934.061);
        Energies.push_back(1120.287);
        Energies.push_back(1238.110);
        Energies.push_back(1377.669);
        Energies.push_back(1509.228);
        Energies.push_back(1729.595);
        Energies.push_back(1764.494);
        Energies.push_back(1847.420);
        Energies.push_back(2118.551);
        Energies.push_back(2204.215);
        Energies.push_back(2447.810);
    }
    if(fSourceName.Contains("241Am",TString::kIgnoreCase))
    {
        fSourceName = "241Am";
        Energies.push_back(  26.345);
        Energies.push_back(  59.537);
    }
}


void GWRecal::numparams(char * cmd, int nn)
{
    if(nn >= 0) return;
    nn = -nn;
    cout<<Form("error decoding command %s ==> %d parameter(s) missing \n", cmd, nn);
    exit(EXIT_FAILURE);
}


bool GWRecal::TestResults()
{
    // delete previous
    memset(testData, 0, sizeof(float)*specLength);

    // fill with calibrated spectrum
    double e0(0), e1(0);  // calibrated beginning and end of channel
    e1 = Calibrated(0)*testGain;
    for(int nn = 0; nn < specLength; nn++) {
        e0 = e1;
        e1 = Calibrated(nn)*testGain;
        double val = specData[nn];
        int n0 = int(e0);
        int n1 = int(e1);
        if(n0 == n1) {        // all in one channel
            if(n0 >=0 && n0 < specLength)
                testData[n0] += float(val);
        }
        else if (n0 < n1) {   // spread linearly in the target range
            val /= (e1-e0);
            if(n0 >=0 && n0 < specLength)
                testData[n0] += float(val*(n0+1-e0));
            for(int ii = n0+1 ; ii < n1; ii++) {
                if(ii >=0 && ii < specLength)
                    testData[ii] += float(val);
            }
            if(n1 >=0 && n1 < specLength)
                testData[n1] += float(val*(e1-n1));
        }
        else {                //  spread linearly in the target range, in reverse order
            val /= (e0-e1);
            {int nt = n0; n0 = n1; n1 = nt;}
            if(n0 >=0 && n0 < specLength)
                testData[n0] += float(val*(n0+1-e1));
            for(int ii = n0+1; ii < n1; ii++) {
                if(ii >=0 && ii < specLength)
                    testData[ii] += float(val);
            }
            if(n1 >=0 && n1 < specLength)
                testData[n1] += float(val*(e0-n1));
        }
    }

    // write result
    int nn = fwrite(testData, sizeof(float), specLength, testFP);

    return (nn == specLength);
}


// - smooth the spectrum with a gaussian kernel
// - differentiate ONCE
// Loop (maxNum times)
//   find the minimum of the derivative and the associated +++--- pattern
//   accept it with conditions with the cross-over as peak position
//   erase the region
// 
int GWRecal::LargePeaks1(float *m_pSpek, int chA, int chB, float fwhm, float minAmpl, std::vector<double>&vPeaks, int maxNum)
{
    //CheckLimits(chA, chB);
    int nChan = chB-chA+1;

    double  sigma = fwhm/s2fwhm;
    float  *tSpek = NULL;     // the smoothed and differentiated spectrum
    int     xMin  = 0;        // position of chA in tSpek
    SmoothSpek(m_pSpek+chA, nChan, sigma, 1, tSpek, xMin);
    int xMax = xMin + nChan;  // position of chB+1 in tSpek

    // ready to search peaks from xMin to xMax
    vPeaks.clear();

    // put here sensible values
    xP_minAmpl     = float(0.5*abs(minAmpl*exp(-.5)/sigma)); // max(g') = |g'(sigma)| = g(0)*exp(-0.5)/sigma
    xP_fThreshCFA  = 0;
    xP_nMinWidthP  = max(2, int(fwhm)-2);
    xP_nMinWidthN  = max(2, int(fwhm)-2);

    while(vPeaks.size() < size_t(maxNum)) {
        // variables named as in xP_Next1
        xP_0 = xP_1 = xP_2 = -1; // positive lobe
        xP_3 = xP_4 = xP_5 = -1; // negative lobe
        // find the largest negative peak
        float yMin = 0;
        for(int nn = xMin+1; nn < xMax; nn++) {
            if(tSpek[nn] < yMin) {
                yMin = tSpek[nn];
                xP_4 = nn;
            }
        }
        if(xP_4 < 0)
            break;      // empty?
        // find the zero crossing at the left of the minimum
        for(int nn = xP_4; nn > xMin; nn--) {
            if(tSpek[nn] > 0) {
                xP_2 = nn;
                break;
            }
        }
        if(xP_2 < 0)
            break;
        xP_3 = xP_2 + 1;
        // find the end of the negative lobe
        for(int nn = xP_3; nn < xMax; nn++) {
            if(tSpek[nn] >= 0) {
                xP_5 = nn;
                break;
            }
        }
        if(xP_5 < 0)
            break;
        // find the beginning of the positive lobe, and the maximum
        float yMax = 0;
        for(int nn = xP_2; nn > xMin; nn--) {
            if(tSpek[nn] <= 0) {
                xP_0 = nn;
                break;
            }
            if(tSpek[nn] > yMax) {
                yMax = tSpek[nn];
                xP_1 = nn;
            }
        }
        if(xP_0 < 0)
            break;
        int widthP = xP_2-xP_0;
        int widthN = xP_5-xP_3;
        bool good = true;
        if(widthP < xP_nMinWidthP || widthN < xP_nMinWidthN)
            good = false;
        if(yMax < xP_minAmpl || yMin > -xP_minAmpl)
            good = false;
        float cmp = abs(yMax/yMin);
        if(cmp < 0.2f || cmp > 2.f)
            good = false;
        if(good) {
            // peak position is transition from positive to negative
            double xPos = 0.5f*(xP_2+xP_3);     // position of peak in tSpek coordinates
            vPeaks.push_back(xPos-xMin+chA+1);  //     ""            m_pSpek   ""
        }
        // zero the region
        for(int nn = xP_0; nn <xP_5; nn++)
            tSpek[nn] = 0;
    }

#if 0
    // Copy back (the remaining part of) the working spectrum
    memcpy(m_pSpek+chA, tSpek+xMin, sizeof(float)*nChan);
    float scale = float(sigma/exp(-.5));
    for(int nn = chA; nn < chA+nChan; nn++)
        m_pSpek[nn] *= scale;
#endif

    delete [] tSpek;

    return vPeaks.size();
}


// imported from TkT
// - smooth the spectrum with a gaussian kernel
// - differentiate ONCE
// Loop
//   look for the +++--- pattern with conditions on its shape and amplitude
//   test significance of peak before accepting it (nothing yet)
// 
int  GWRecal::PeakSearch1(float *m_pSpek, int chA, int chB, float fwhm, float minAmpl, std::vector<double>&vPeaks)
{
    //CheckLimits(chA, chB);
    int nChan = chB-chA+1;

    double  sigma = fwhm/s2fwhm;
    float  *tSpek = NULL;     // the smoothed and differentiated spectrum
    int     xMin  = 0;        // position of chA in tSpek
    SmoothSpek(m_pSpek+chA, nChan, sigma, 1, tSpek, xMin);
    int xMax = xMin + nChan;  // position of chB+1 in tSpek

    // ready to search peaks from xMin to xMax
    vPeaks.clear();

    // put here sensible values
    xP_minAmpl     = float(0.5*abs(minAmpl*exp(-.5)/sigma)); // max(g') = |g'(sigma)| = g(0)*exp(-0.5)/sigma
    xP_fThreshCFA  = 0;
    xP_nMinWidthP  = max(2, int(fwhm)-2);
    xP_nMinWidthN  = max(2, int(fwhm)-2);

    int xx = xMin;
    while(xx < xMax) {
        if( xP_Next1(tSpek+xx, nChan-xx) <= 0 )
            break;
        // peak position is transition from positive to negative
        double xPos = xx + 0.5f*(xP_2+xP_3);      // position of peak in tSpek coordinates
        vPeaks.push_back(xPos + chA + 1 - xMin);  //     ""            m_pSpek   ""
        xx += xP_5;
    }

#if 0
    memcpy(m_pSpek+chA, tSpek+xMin, sizeof(float)*nChan);
    float scale = float(sigma/exp(-.5));
    for(int nn = chA; nn < chA+nChan; nn++)
        m_pSpek[nn] *= scale;
#endif

    delete [] tSpek;

    return vPeaks.size();
}


// imported from TkT
// - smooth the spectrum with a gaussian kernel
// - differentiate TWICE
// Loop
//   look for the next mexican-hat pattern with conditions on its shape and amplitude
//   test significance of peak before accepting it (nothing yet)
// 
int GWRecal::PeakSearch2(float *m_pSpek, int chA, int chB, float fwhm, float minAmpl, std::vector<double>&vPeaks)
{
    //CheckLimits(chA, chB);
    int nChan = chB-chA+1;

    double  sigma = fwhm/s2fwhm;
    float  *tSpek = NULL;     // the smoothed and differentiated spectrum
    int     xMin  = 0;        // position of chA in tSpek
    SmoothSpek(m_pSpek+chA, nChan, sigma, 2, tSpek, xMin);
    int xMax = xMin + nChan;  // position of chB+1 in tSpek

    // ready to search peaks from xMin to xMax
    vPeaks.clear();

    // put here sensible values
    xP_minAmpl     = float(abs(minAmpl/(sigma*sigma))); // max(g'') =  |g''(0)| = |g(0)|/sigma**2
    xP_fThreshCFA  = 0;
    xP_nMinWidthN  = max(2, int(fwhm)-2);
    xP_nMinWidthP1 = max(2, xP_nMinWidthN-1);
    xP_nMinWidthP2 = max(2, (xP_nMinWidthN+1)/2);

    int xx = xMin;
    while(xx < xMax) {
        if( xP_Next2(tSpek+xx, xMax-xx) <= 0 )
            break;
        // first moments of the negative lobe to get the peak position
        float xP0 = 0, xP1 = 0;
        int nn0 = xx + xP_2 + 1;
        int nn1 = xx + xP_4 + 1;
        for(int nn = nn0; nn < nn1; nn++) {
            float x = float(nn-nn0)+0.5f;
            float y = tSpek[nn];
            xP0 += y;
            xP1 += y*x;
        }
        float xPos = nn0 + xP1/xP0;           // position of peak in tSpek coordinates
        vPeaks.push_back(xPos + chA - xMin);  //     ""            m_pSpek   ""
        xx += xP_4;
    }

#if 0
    memcpy(m_pSpek+chA, tSpek+xMin, sizeof(float)*nChan);
    float scale = float(sigma*sigma);
    for(int nn = chA; nn < chA+nChan; nn++)
        m_pSpek[nn] *= scale;
#endif

    delete [] tSpek;

    return vPeaks.size();
}


// generates a smoothed and differentiated version of the input data
// returns the  pointer and the position of the first valid channel
// the caller has to free tSpek
void GWRecal::SmoothSpek(float *spek, int nchan, double sigma, int ndiff, float *&tSpek, int &start)
{
    // generate the gaussian kernel for the smoothing operation
    int     nKern  = 2*int(5.5*sigma)+1;    // < 1ppm
    int     zKern  =  nKern/2;
    double  gfact1 = -1./(2.*sigma*sigma);
    double  gfact2 = 1./(sigma*sqrt(2.*m_pi));

    float  *gKern  = new float [nKern];
    for(int nn = 0; nn <= zKern; nn++) {
        double gvalue = gfact2*exp( ((nn-zKern)*(nn-zKern))*gfact1 );
        gKern[nn] = gKern[nKern-1-nn] = float(gvalue);
    }

    // the temporary space for the data
    int tLen = nKern + nchan + nKern;

    start = nKern;
    tSpek = new float [tLen];

    // copy the original spectrum to the temporary storage,
    // leaving nKern empty channels on both sides
    memset(tSpek, 0, sizeof(float)*nKern);
    memcpy(tSpek+nKern, spek, nchan*sizeof(float));
    memset(tSpek+nKern+nchan, 0, sizeof(float)*nKern);

    // convolution --> the resulting spectrum is right shifted by zKern channels
    for(int ii = tLen-1; ii >= nKern; ii--) {
        float csum = 0;
        for(int jj = 0; jj < nKern; jj++)
            csum += gKern[jj]*tSpek[ii-jj];
        tSpek[ii] = csum;
    }
    start += zKern;

    // differentiate as requested
    for(int nd = 0; nd < ndiff; nd++) {
        for(int ii = tLen-1; ii >0; ii--)
            tSpek[ii] -= tSpek[ii-1];
    }
    //start += ndiff/2;

    delete [] gKern;
}


// trigger on a +++--- sequence, first derivative of a gaussian
// with conditions on the width and amplitude of the lobes
int GWRecal::xP_Next1(float *yDat, int yLen)
{
    enum {
        stateUnknown,   // initial or negative before starting a new sequence
        stateP,         // positive: xP_0(first positive chan),xP_1(max), xP_2(last positive)
        stateN         // negative: xP_3(first negative), xP_4(min), xP_5(last negative)
    };

    int    state = stateUnknown;            // start in an indefinite state
    int    widthP = 0, widthN = 0;
    float  yMax=0;
    float  yMin=0;

    for(int nn = 0; nn < yLen; nn++) {
        float yInp = yDat[nn];
        switch (state) {
        case stateUnknown:                  // was below CF threshold in an indefinite state
            if(yInp > xP_fThreshCFA) {        // and goes above
                state = stateP;                 // initiate a possible trigger sequence
                xP_0  = xP_1 = xP_2 = nn;
                yMax = yInp;
                widthP = 1;                     // start checking width
            }                                 // if it stays below there is nothing to do
            break;
        case stateP:                        // was above threshold
            if(yInp >= xP_fThreshCFA) {       // and stays above
                if(yInp > yMax) {
                    xP_1 = nn;                    // record maximum of first positive lobe
                    yMax = yInp;
                }
                xP_2 = nn;
                widthP++;                       // incr positive width and stay in this state
            }
            else {                            // transition to negative
                if(widthP >=  xP_nMinWidthP &&  // positive lobe wide enough
                        //widthP < 4*xP_nMinWidthP &&  // but not too wide
                        xP_1 > xP_0 &&               // and has a maximum
                        yMax > xP_minAmpl) {         // which is large enough
                    state = stateN;             // initiate a negative lobe sequence
                    xP_3  = xP_4 = xP_5 = nn;
                    yMin  = yInp;
                    widthN = 1;                // start counting width of negative lobe
                }
                else
                    state = stateUnknown;         // go back to indefinite below threshold
            }
            break;
        case stateN:                        // below threshold in a valid-sequence state
            if(yInp < xP_fThreshCFA) {        // and stays below
                if(yInp < yMin) {
                    xP_4 = nn;                    // record minimum
                    yMin = yInp;
                }
                xP_5 = nn;
                widthN++;                       // incr negative width and stay in this state
                if(widthN > xP_nMinWidthN &&    // second positive lobe wide enough for a valid sequence
                        yInp > yMin &&                // and signal is beyond minimum
                        yMin < -xP_minAmpl) {         // and minimum large enough
                    float cmp = abs(yMax/yMin);
                    if(cmp < 0.2f || cmp > 2.f)
                        state = stateUnknown;
                    else
                        return nn;                  // this is a good trigger
                }
                break;
            }
            else {                            // goes positive without having triggered
                state = stateP;                 // initiate a new possible trigger sequence
                xP_0 = xP_1 = xP_2 = nn;
                yMax = yInp;
                widthP = 1;                     // start checking width
                break;
            }
        }
    }
    return -1;
}


// trigger on a +++---++++ sequence
// with conditions on the width and amplitude of the lobes
int GWRecal::xP_Next2(float *yDat, int yLen)
{
    enum {
        stateUnknown,   // initial or negative before starting a new sequence
        stateP1,        // positive: xP_0(first positive chan),xP_1(max), xP_2(last positive)
        stateN,         // negative: xP_3(min), xP_4(last negative)
        stateP2         // positive: xP_5(max), xP_6(last)
    };

    int    state = stateUnknown;            // start in an indefinite state
    int    widthP1 = 0, widthN = 0, widthP2 = 0;
    float  yMax1=0;
    float  yMin=0;
    float  yMax2=0;

    for(int nn = 0; nn < yLen; nn++) {
        float yInp = yDat[nn];
        switch (state) {
        case stateUnknown:                  // was below CF threshold in an indefinite state
            if(yInp > xP_fThreshCFA) {        // and goes above
                state = stateP1;                // initiate a possible trigger sequence
                xP_0  = xP_1 = xP_2 = nn;
                yMax1 = yMin = yInp;
                widthP1 = 1;                    // start checking width
            }                                 // if it stays below there is nothing to do
            break;
        case stateP1:                       // was above threshold
            if(yInp >= xP_fThreshCFA) {       // and stays above
                if(yInp > yMax1) {
                    xP_1 = nn;                    // record maximum of first positive lobe
                    yMax1 = yInp;
                }
                xP_2 = nn;
                widthP1++;                      // incr positive width and stay in this state
            }
            else {                            // transition to negative
                if(widthP1 >= xP_nMinWidthP1 && // positive lobe wide enough
                        xP_1 > xP_0) {              // and has a maximum
                    state = stateN;             // initiate a negative lobe sequence
                    xP_3  = xP_4 = nn;
                    yMin    = yInp;
                    widthN  = 1;                // start counting width of negative lobe
                }
                else
                    state = stateUnknown;         // go back to indefinite below threshold
            }
            break;
        case stateN:                        // below threshold in a valid-sequence state
            if(yInp < xP_fThreshCFA) {        // and stays below
                if(yInp < yMin) {
                    xP_3 = nn;                    // record minimum
                    yMin = yInp;
                }
                xP_4 = nn;
                widthN++;                       // incr positive width and stay in this state
            }
            else {                            // transition to positive
                if(widthN >= xP_nMinWidthN &&   // negative lobe wide enough
                        yMin < -xP_minAmpl &&         // its amplitude large enough
                        yMin < -1.1*yMax1 &&          // and minimum-to-maximum is acceptable
                        yMin > -4.5*yMax1) {          // (exact value is -0.5*exp(1.5)) = -2.2408)
                    state = stateP2;            // initiate the second positive lobe sequence
                    xP_5  = xP_6 = nn;
                    yMax2 = yInp;
                    widthP2 = 1;                // start counting width of second positive lobe
                }
                else {
                    state = stateP1;              // restart a new sequence
                    xP_0  = xP_1 = xP_2 = nn;
                    yMax1 = yMin =  yInp;
                    widthP1 = 1;                  // start checking width
                }
            }
            break;
        case stateP2:                       // above threshold in a valid-sequence state
            if(yInp > xP_fThreshCFA) {        // and is still above
                if(yInp > yMax2) {
                    xP_5  = nn;                   // record max of second positive lobe
                    yMax2 = yInp;
                }
                xP_6 = nn;
                widthP2++;                      // incr positive width and stay in this state
                if(widthP2 > xP_nMinWidthP2 &&  // second positive lobe wide enough for a valid sequence
                        yInp < yMax2) {              // and signal is beyond maximum
                    return nn;                 // this is a good trigger
                }
            }
            else {                            // goes positive (with or without trigger)
                state = stateUnknown;           // ready  to start a possible new trigger sequence
            }
            break;
        }
    }
    return -1;
}


int GWRecal::FitPeaks(int verbose)
{
    if(specPeaks.size() < 1)
        return 0;

    // order peaks (needed to merge correctly close neighbours)
    sort( specPeaks.begin( ), specPeaks.end( ) );

    const float nwLe = 4.f; //  define left border
    const float nwOv = 2.f; //  check region overlap
    const float nwRi = 3.f; //  define right border

    Peaks.clear();

    // instead of specFWHMdef, we should use realistic width to define regions and overlaps

    int nmult = 0;
    int numpk = int(specPeaks.size());
    for(int nn = 0; nn < numpk; nn += 1 + nmult) {
        int chanA = int(specPeaks[nn]-nwLe*fCurrentspecFWHMdef);
        int chanB = int(specPeaks[nn]+nwRi*fCurrentspecFWHMdef);
        chanA = max(1,chanA);
        chanB = min(chanB, specLength-2);
        nmult = 0;
        if(nn < numpk-2) {
            // check if the next peaks are in the same region
            for(size_t jj = nn+1; jj < specPeaks.size(); jj++) {
                if(specPeaks[jj]-nwOv*fCurrentspecFWHMdef >= chanB)
                    break;
                nmult++;
                chanB = int(specPeaks[jj]+nwRi*fCurrentspecFWHMdef);
                chanB = min(chanB, specLength-2);
            }
        }

        int np = m_pCFit->CalcGaussianFit(specData, chanA, chanB, specPeaks);

        for(int jj = 0; jj < np; jj++) {
            Fitted res;

            res.NSubPeaks = np;
            res.BgFrom= m_pCFit->BgFrom();
            res.BgTo  = m_pCFit->BgTo();
            res.BgdOff= m_pCFit->Back0();
            res.BgdSlope= m_pCFit->Back1();

            //      cout<<"[ "<<res.BgFrom<<" ; "<<res.BgTo<<" ] ==> Bg = "<<res.BgdSlope<<"*x + "<<res.BgdOff<<endl;

            res.area   = m_pCFit->Area(jj);
            res.ampli  = m_pCFit->Amplitude(jj);
            res.posi   = m_pCFit->Position(jj);
            res.fw05   = m_pCFit->Fw05(jj);
            res.fw01   = m_pCFit->Fw01(jj);
            res.fwhm   = m_pCFit->Fwhm(jj);
            res.tailL  = m_pCFit->W01L(jj);
            res.tailR  = m_pCFit->W01R(jj);
            res.Lambda = m_pCFit->TailLeft(jj);
            res.Rho    = m_pCFit->TailRight(jj);
            res.S      = m_pCFit->Step(jj);

            bool bad = (res.area < 10) || (res.fwhm >= 5*fCurrentspecFWHMdef) || (res.tailR > 5.f*res.tailL);
            if(verbose&1) {
                if(nn || jj) cout<<Form("#1                  ");
                if(jj == 0) cout<<Form(" %3d ", nn+jj);
                else        cout<<Form("  %3d", nn+jj);
                cout<<Form(" %10.1f %10.2f  %10.2f %8.3f %8.3f %8.3f %8.3f %8.3f", res.area, res.ampli, res.posi, res.fw05, res.fw01, res.fwhm, res.tailL, res.tailR);
                if(bad)
                    cout<<Form("  ***");
                cout<<Form("\n");
            }
            if(!bad) {
                Peaks.push_back(res);
            }
        }
    }
    return Peaks.size();
}


// the energies should be already sorted in ascending order
// the peaks will be sorted in largest area order and the 
// first Energies.size() of them will be connected to all energies, checking
// how many other peaks match the specific gain
double GWRecal::ECalibration(int nspe, double& offset1_, double& slope1_, double& offset2_, double& slope2_, double& curv2_, int verbose_)
{
    vector <Fitted>::iterator Iter;

    if(Energies.size() == 1) {
        // select the peak with the largest area
        int pn_max = 0;
        for(size_t nn = 1; nn < Peaks.size(); nn++) {
            if(Peaks[nn].area > Peaks[pn_max].area) {
                pn_max = nn;
            }
        }
        offset1_ = slope1_ = 0;
        double slope_  = Energies[0] / Peaks[pn_max].posi;
        Peaks[pn_max].good = true;
        Peaks[pn_max].eref = 0;
        if(verbose_&2)
            cout<<Form("#2 %5d  Slope = %10.6f\n", nspe, slope_);
        return slope_;
    }

    if( Peaks.size() < 2 || Energies.size() < 2) {
        if(verbose_&2)
            cout << "#   Number of peaks (" << Peaks.size() << ") or number of energies (" << Energies.size() << ")  too small" << endl;
        return 0;
    }

#if 0
    // all peaks used as pivot
    int npmax = Peaks.size();
    for(int np = 0; np < npmax; np++)
        Peaks[np].good = true;
#else
    // only the peak with largest area used as pivot
    sort( Peaks.begin( ), Peaks.end( ), largerAmplitude() );
    int npmax = min( min(Peaks.size(), size_t(4)), Energies.size() );  // limit the number of peaks to min(np,ne,4)
    for(int np = 0; np < npmax; np++)
        Peaks[np].good = true;
    for(size_t np = npmax; np < Peaks.size(); np++)
        Peaks[np].good = false;
#endif

    // (re)order the peaks
    sort( Peaks.begin( ), Peaks.end( ), smallerPosi() );
    if(verbose_&4) {
        cout<<Form("#  Sorted Peak Positions = (");
        for ( Iter = Peaks.begin( ) ; Iter != Peaks.end( ) ; Iter++ ) cout<<Form( " %d", int((*Iter).posi) );
        cout<<Form(" )\n");
    }

    // Connect every (large) peak with every energy and check how many other (peak,energy) pairs agree with this slope
    // The combination with the largest number of agreeing pairs provides the reference gain for the final average
    // Equal number of matches are ordered by chi2
    int    np_max = 0;      // number of matches
    int    np_ind = 0;      // reference peak
    int    np_ref = 0;      // reference energy
    double np_chi = 1.e40;  // to distinguish among equal matchers
    double dpos2max = fCurrentspecFWHMdef*fCurrentspecFWHMdef;
    for(size_t np = 0; np < Peaks.size(); np++) {
        if(!Peaks[np].good)
            continue;
        int    ec_max = 0;      // number of matches for peak np
        int    ec_ref = 0;      // reference energy
        double ec_chi = 1.e30;  // its chi2
        for(size_t ne = 0; ne < Energies.size(); ne++) {            // connecting peak np with energy ne
            double lgain = (Peaks[np].posi)/Energies[ne];             // gives this gain
            int match = 0;
            double chi2 = 0;
            for(size_t nee = 0; nee < Energies.size(); nee++) {       // count the number of matches for this combination
                double epos = Energies[nee]*lgain;                      // expected position
                for(size_t ip = 0; ip < Peaks.size(); ip++) {
                    double dpos  = Peaks[ip].posi-epos;
                    double dpos2 = dpos*dpos;
                    if(dpos2 < dpos2max) {
                        match++;
                        chi2 += dpos2;
                        break;
                    }
                }
            }
            if( (match > ec_max)  ||                        // more matches
                    ((match == ec_max) && (chi2 < ec_chi)) ) {  // or better chi2
                ec_max = match;       // record the best combination so far
                ec_ref = ne;
                ec_chi = chi2;
            }
        }
        if( (ec_max > np_max) ||                           // more matches
                ((ec_max == np_max) && (ec_chi < np_chi)) ) {  // or better chi2
            np_max = ec_max;       // record the best combination so far
            np_ind = np;
            np_ref = ec_ref;
            np_chi = ec_chi;
        }
        //if(np_max == Energies.size())
        //  break;    // cannot do better than this ?
    }

    double bestSlope = Energies[np_ref]/(Peaks[np_ind].posi);
    if(verbose_&4) {
        cout<<Form("#  Best-match slope %g [p=%d e=%d] with %d values\n", bestSlope, np_ind, np_ref, np_max);
    }

    // find the good peaks
    for(size_t np = 0; np < Peaks.size(); np++) {
        Peaks[np].eref = -1;
        Peaks[np].good = false;
    }
    int match = 0;
    for(size_t ne = 0; ne < Energies.size(); ne++) {
        double epos = Energies[ne]/bestSlope;
        for(size_t np = 0; np < Peaks.size(); np++) {
            if(abs(Peaks[np].posi-epos) < fCurrentspecFWHMdef && !Peaks[np].good) {
                match++;
                Peaks[np].eref = ne;
                Peaks[np].good = true;
                break;
            }
        }
    }

    // fit slope using the good peaks, starting without errors
    double val_xx = 0;
    double val_xy = 0;
    int used = 0;
    for(size_t np = 0; np < Peaks.size(); np++) {
        if(Peaks[np].good) {
            ++used;
            double x  = Peaks[np].posi;
            double y  = Energies[Peaks[np].eref];
            val_xx += x*x;
            val_xy += x*y;
        }
    }
    double slope0 = val_xy / val_xx;
    double slope_  = slope0;

    // if enabled, error**2 for chi1 is: errPos**2 + errEner**2
    // errP = slope * sigma/sqrt(Area)  Approx for a gaussian peak. Scaled to keV
    // errE taken from command line -errors errE (should be specific of each calibration line)

    // Not really worthwhile in actual practice.

    if(useErrors) {
        // repete, using errors
        val_xx = 0;
        val_xy = 0;
        used = 0;
        for(size_t np = 0; np < Peaks.size(); np++) {
            if(Peaks[np].good) {
                ++used;
                double x  = Peaks[np].posi;
                double y  = Energies[Peaks[np].eref];
                double errP = slope0* pow(Peaks[np].fwhm/s2fwhm, 2.)/Peaks[np].area;
                double w2 = 1. / (errP*errP+errE*errE);
                val_xx += w2*x*x;
                val_xy += w2*x*y;
            }
        }
        slope_ = val_xy / val_xx;
    }

    bool fit1Done = false;
    offset1_ = slope1_ = 0;
    if(match>=2 && doPoly1) {
        // fit a stright line through the good peaks
        double sx0 = 0., sx1 = 0., sx2 = 0., syx0 = 0., syx1 = 0.;
        double w2s = 0; // to record the average weight
        used = 0;
        double w2 = 1.;
        for(size_t np = 0; np < Peaks.size(); np++) {
            if(Peaks[np].good) {
                ++used;
                double x  = Peaks[np].posi;
                double y  = Energies[Peaks[np].eref];
                if(useErrors) {
                    double errP = slope_* pow(Peaks[np].fwhm/s2fwhm, 2.)/Peaks[np].area;
                    w2 = 1. / (errP*errP+errE*errE);
                }
                sx0  += w2;
                sx1  += w2*x;
                sx2  += w2*x*x;
                syx0 += w2*y;
                syx1 += w2*y*x;
                w2s  += w2;
            }
        }
        if(numZero > 0 ) {
            // multiple artificial peaks at (0,0) to help system pass through origin
            //double x=0, y=0;      // --> only sx0 affected
            w2 = w2s/used;   // ??
            sx0  += w2*numZero;     // with multiplicity
        }
        double deter = sx0*sx2-sx1*sx1;
        if(deter > 0.) {
            offset1_ = ( sx2*syx0 - sx1*syx1)/deter;
            slope1_  = (-sx1*syx0 + sx0*syx1)/deter;
//            double err_offset1 = sqrt(sx2/deter);
//            double err_slope1  = sqrt(sx0/deter);
            fit1Done = true;
        }
    }

    bool fit2Done = false;
    offset2_ = slope2_ = curv2_ = 0;
    if(match>=3 && doPoly2) {
        // fit a parabola line through the good peaks
        double sx0 = 0., sx1 = 0., sx2 = 0., sx3 = 0., sx4 = 0., syx0 = 0., syx1 = 0., syx2 = 0.;
        double w2s = 0; // to record the average weight
        used = 0;
        double w2 = 1.;
        for(size_t np = 0; np < Peaks.size(); np++) {
            if(Peaks[np].good) {
                ++used;
                double x  = Peaks[np].posi;
                double y  = Energies[Peaks[np].eref];
                if(useErrors) {
                    double errP = slope_* pow(Peaks[np].fwhm/s2fwhm, 2.)/Peaks[np].area;
                    w2 = 1. / (errP*errP+errE*errE);
                }
                sx0  += w2;
                sx1  += w2*x;
                sx2  += w2*x*x;
                sx3  += w2*x*x*x;
                sx4  += w2*x*x*x*x;
                syx0 += w2*y;
                syx1 += w2*y*x;
                syx2 += w2*y*x*x;
                w2s  += w2;
            }
        }
        if(numZero > 0 ) {
            // multiple artificial peaks at (0,0) to help system pass through origin
            //double x=0, y=0;      // --> only sx0 affected
            w2 = w2s/used;   // ??
            sx0  += w2*numZero;     // with multiplicity
        }
        double dd[9]; memset(dd, 0, sizeof(dd));    // row-wise
        double xd[3]; memset(xd, 0, sizeof(xd));
        // [sx0]  [sx1]  [sx2]   [b0]    [ex0]
        // [sx1]  [sx2]  [sx3] * [b1] =  [ex1]
        // [sx2]  [sx3]  [sx4]   [b2]    [ex3]

        dd[0] = sx0; dd[1] = sx1; dd[2] = sx2; xd[0] = syx0;
        dd[3] = sx1; dd[4] = sx2; dd[5] = sx3; xd[1] = syx1;
        dd[6] = sx2; dd[7] = sx3; dd[8] = sx4; xd[2] = syx2;

        double DD[9];
        if(InvertMatrix3(dd, DD)) {
            offset2_ = DD[0]*xd[0] + DD[1]*xd[1] + DD[2]*xd[2];
            slope2_  = DD[3]*xd[0] + DD[4]*xd[1] + DD[5]*xd[2];
            curv2_   = DD[6]*xd[0] + DD[7]*xd[1] + DD[8]*xd[2];
            fit2Done = true;
        }
    }

    if(verbose_&2) {
        cout<<Form("#2 %5d  Slope = %1.6f", nspe, slope_);
        if(fit1Done)
            cout<<Form("    Cal1=[ %7.4f  %1.6f ]", offset1_, slope1_);
        if(fit1Done && fit2Done)
            cout<<Form("    Cal2=[ %7.4f  %1.6f % .5e ]", offset2_, slope2_, curv2_);
        cout<<Form("\n");
        size_t nn = 0;
        double chi2[5] = {0};
        int nused = 0;
        for(size_t np = 0; np < Peaks.size(); np++) {
            if(Peaks[np].good) {
                nused++;
                double lineEner = Energies[Peaks[np].eref];
                double seenPosi = Peaks[np].posi;
                double seenEner = seenPosi*slope_;
                double diffEner = seenEner-lineEner;
                chi2[0] += diffEner*diffEner;
                cout<<Form("#2 %5d  %4d %4d %15.1f %10.2f %9.3f %9.3f   %11.3f %7.3f",
                           nspe, int(nn++), int(np), Peaks[np].area, Peaks[np].posi, Peaks[np].fwhm, Peaks[np].fwhm*slope_, lineEner, diffEner);
                if(fit1Done) {
                    seenEner = offset1_ + slope1_*seenPosi;
                    diffEner = seenEner-lineEner;
                    cout<<Form(" %7.3f", diffEner);
                    chi2[1] += diffEner*diffEner;
                }
                if(fit1Done && fit2Done) {
                    seenEner = offset2_ + slope2_*seenPosi + curv2_*seenPosi*seenPosi;
                    diffEner = seenEner-lineEner;
                    cout<<Form(" %7.3f", diffEner);
                    chi2[2] += diffEner*diffEner;
                }
                cout<<Form("\n");
            }
        }
        chi2[0] = (nused-1) ? chi2[0]/(nused-1) : 0;
        cout<<Form("#2 %5d  Chi2 %73.3f", nspe, chi2[0]);
        if(fit1Done) {
            chi2[1] = (nused-2) ? chi2[1]/(nused-2) : 0;
            cout<<Form(" %7.3f", chi2[1]);
        }
        if(fit1Done && fit2Done) {
            chi2[2] = (nused-3) ? chi2[2]/(nused-3) : 0;
            cout<<Form(" %7.3f", chi2[2]);
        }
        cout<<Form("\n");
    }

    return slope_;
}


double GWRecal::TCalibration()
{
    return 0;
}


bool GWRecal::InvertMatrix3(const double m[9], double invOut[9])
{
    double d00 = m[0];
    double d01 = m[1];
    double d02 = m[2];
    double d10 = m[3];
    double d11 = m[4];
    double d12 = m[5];
    double d20 = m[6];
    double d21 = m[7];
    double d22 = m[8];

    double Det =  d00*(d11*d22 - d12*d21) + d01*(d12*d20 - d10*d22) + d02*(d10*d21 - d11*d20);
    if (Det == 0)
        return false;

    double D00 =  d11*d22 - d12*d21; double D01 =  d12*d20 - d10*d22; double D02 =  d10*d21 - d11*d20;
    double D10 =  d21*d02 - d22*d01; double D11 =  d22*d00 - d20*d02; double D12 =  d20*d01 - d21*d00;
    double D20 =  d01*d12 - d02*d11; double D21 =  d02*d10 - d00*d12; double D22 =  d00*d11 - d01*d10;

    Det = 1/Det;

    invOut[0] = D00 * Det;
    invOut[1] = D01 * Det;
    invOut[2] = D02 * Det;
    invOut[3] = D10 * Det;
    invOut[4] = D11 * Det;
    invOut[5] = D12 * Det;
    invOut[6] = D20 * Det;
    invOut[7] = D21 * Det;
    invOut[8] = D22 * Det;

#if 0
#define indx(n1, n2) (n1*3+n2)
    double unit[9]; memset(unit, 0, sizeof(unit));
    for(int i1 = 0; i1 < 3; i1++) {
        for(int i2 = 0; i2 < 3; i2++) {
            double accu = 0;
            for(int ii = 0; ii < 3; ii++) {
                accu += m[indx(i1,ii)]*invOut[indx(ii,i2)];
            }
            unit[indx(i1,i2)] = accu;
        }
    }
#undef indx
#endif

    return true;
}


double GWRecal::Calibrated(double x)
{
    if(!doSlope)     return slope*x;
    else if(doPoly2) return offset1 + slope1*x;
    else if(doPoly1) return offset2 + (slope2 + curv2*x)*x;
    else             return tshift  + x;
}


void GWRecal::SetHistChannelLimits(int SpecNbr, float ChFrom, float ChTo)
{
    bool IndexAlreadyDefined = false;
    for(size_t i=0 ; i<nlim.size() ; i++)
    {
        NLimits_t lim = nlim[i];
        if(lim.index == SpecNbr)
        {
            lim.from = (int)ChFrom;
            lim.to = (int)ChTo;
            IndexAlreadyDefined = true;
        }
    }
    if(!IndexAlreadyDefined)
    {
        NLimits_t lim;
        lim.index = SpecNbr;
        lim.from = (int)ChFrom;
        lim.to = (int)ChTo;
        nlim.push_back(lim);
    }
}


void GWRecal::SetHistPeaksLimits(int SpecNbr, float FWHM, float Ampli)
{
    bool IndexAlreadyDefined = false;
    for(size_t i=0 ; i<nwa.size() ; i++)
    {
        NWA_t tnwa = nwa[i];
        if(tnwa.index == SpecNbr)
        {
            tnwa.width = FWHM;
            tnwa.ampli = Ampli;
            IndexAlreadyDefined = true;
        }
        nwa[i] = tnwa;
    }
    if(!IndexAlreadyDefined)
    {
        NWA_t tnwa;
        tnwa.index = SpecNbr;
        tnwa.width = FWHM;
        tnwa.ampli = Ampli;
        nwa.push_back(tnwa);
    }
}