FitSpectra.cpp

Go to the documentation of this file.

//typedef char __signed;

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>
#include <map>

#include "FitSpectra.h"

using namespace std;

//ClassImp(FitSpectra);

FitSpectra::FitSpectra(TString type, TString source)
{
    fHistToFit = 0x0;
    fHistName = "";
    fSourceName = source;
    fFitType = type;
    fFitFctName = "";
    fBackGroundName = "";
    fDefinePersoBackground = false;
    fUseOffset = false;

    fWidth = 0.;
    fThreshold = 0.;

    InitSourceProperties();

    fxmin = 0.;
    fxmax = 0.;
}

FitSpectra::FitSpectra(TH1 *hist, TString type, TString source)
{
    fHistToFit = hist;
    fSourceName = source;
    fFitType = type;
    fHistName = fHistToFit->GetName();

    fDefinePersoBackground = false;
    fUseOffset = false;

    fFitFctName = "FitFct_" + (TString)fHistToFit->GetName();
    fBackGroundName = "BackGround_" + (TString)fHistToFit->GetName();

    InitSourceProperties();
}

FitSpectra::~FitSpectra()
{
    delete BackGroundHist;
    delete fFitFct;
}

void FitSpectra::InitSourceProperties()
{
    if(fSourceName == "60Co")
    {
        fEPeaks.push_back(1173.238);
        fEPeaks.push_back(1332.513);
    }
    else if(fSourceName == "22Na") {
        fEPeaks.push_back( 511.006);
        fEPeaks.push_back(1274.545);
    }
    else if(fSourceName == "152Eu") {
        fEPeaks.push_back( 121.7793);
        fEPeaks.push_back( 244.6927);
        fEPeaks.push_back( 344.2724);
        fEPeaks.push_back( 411.111);
        fEPeaks.push_back( 443.979);
        fEPeaks.push_back( 778.890);
        fEPeaks.push_back( 867.388);
        fEPeaks.push_back( 964.014);
        fEPeaks.push_back(1085.793);
        //fEPeaks.push_back(1089.700);
        fEPeaks.push_back(1112.070);
        fEPeaks.push_back(1212.95);
        fEPeaks.push_back(1299.1);
        fEPeaks.push_back(1407.993);
    }
    else if(fSourceName == "88Y") {
        fEPeaks.push_back( 898.045);
        fEPeaks.push_back(1836.062);
        //fEPeaks.push_back(2734.087);
    }
    else if(fSourceName ==  "133Ba") {
        fEPeaks.push_back(  53.156);
        fEPeaks.push_back(  79.623);
        fEPeaks.push_back(  80.999);
        fEPeaks.push_back( 160.609);
        fEPeaks.push_back( 223.116);
        fEPeaks.push_back( 276.404);
        fEPeaks.push_back( 302.858);
        fEPeaks.push_back( 356.014);
        fEPeaks.push_back( 383.859);
    }
    else if(fSourceName == "baseline") {
        fEPeaks.push_back(4000.);
    }

    fNPeaksToFit = fEPeaks.size();
}

void FitSpectra::SetSource(TString source)
{
    fSourceName = source;
    InitSourceProperties();
}

void FitSpectra::SetHistogram(TH1 *hist)
{
    fHistToFit = hist;
    fHistName = fHistToFit->GetName();
}

void FitSpectra::SetHistogram(TString HistName)
{
    fHistToFit = (TH1*)gROOT->FindObject(HistName);
    if(fHistToFit == 0x0)
        Warning("FitSpectra::SetHistogram","Histogram with name %s not found",HistName.Data());
    else fHistName = fHistToFit->GetName();
}

void FitSpectra::GetCurrentHistogram()
{
    CurrentPad = gPad;

    TList *ListOfPrimitives = (TList*) CurrentPad->GetListOfPrimitives();

    if(ListOfPrimitives->GetSize())
    {
        TIter iter(ListOfPrimitives);
        TObject *o;
        TH1 *h = 0x0;

        while( (o = iter()) )
        {
            if(o->InheritsFrom("TH1"))
            {
                h = (TH1*)o;
                break;
            }
        }
        fHistToFit = h;
        fHistName = fHistToFit->GetName();

        fFitFctName = "FitFct_" + (TString)fHistToFit->GetName();
        fBackGroundName = "BackGround_" + (TString)fHistToFit->GetName();
    }
}

void FitSpectra::Fit()
{

    FitStatus = 0;

    if(fHistToFit==0x0 || fHistToFit->Integral(2,fHistToFit->GetNbinsX()) < 100)
    {
        fHistToFit=0x0;
        return;
    }

    CurrentPad->SetEditable(true);

    gErrorIgnoreLevel = kError;

    gPad->SetLogy();

    if(fxmin==0. && fxmax == 0.) fHistToFit->SetAxisRange(0.,fHistToFit->GetXaxis()->GetXmax());
    else fHistToFit->SetAxisRange(fxmin,fHistToFit->GetXaxis()->GetXmax());

    int peakfactor = 10;
    if(fSourceName=="baseline") peakfactor = 1;

    TSpectrum *s = new TSpectrum(peakfactor*fNPeaksToFit);

//    if(fWidth!=0. && fThreshold!=0.)
//    {
//        fNFound = s->Search(fHistToFit,fWidth,"nodraw",fThreshold);
//    }
//    else
//    {
        fNFound = 100;
        int test=0;
        double FirstStep = fThreshold;
        while(fNFound>(int)fNPeaksToFit*3)
        {
            fNFound = s->Search(fHistToFit,fWidth,"nodraw",FirstStep);
            FirstStep += fThreshold;
            test++;
            if(test>100)break;
        }
//    }

    gPad->Update();
    gPad->Modified();

//    float *xpeak = s->GetPositionX();
//    float *ypeak = s->GetPositionY();

    for(int i=0 ; i<fNFound ; i++)
    {
        fXPeaks.push_back(/*xpeak[i]*/s->GetPositionX()[i]);
        fYPeaks.push_back(/*ypeak[i]*/s->GetPositionY()[i]);
    }

    for(int t=1 ; t<fNFound ; t++)
        for(int i=0 ; i<fNFound-1 ; i++)
            if(fXPeaks[i]>fXPeaks[i+1])
            {
                swap(fXPeaks[i],fXPeaks[i+1]);
                swap(fYPeaks[i],fYPeaks[i+1]);
            }

    if(fNFound!=fNPeaksToFit)
    {
        int min=0;
        int max = fNFound-1;


        double rapmaxmin = fEPeaks[fEPeaks.size()-1]/fEPeaks[0];
        double rapmaxmin2 = fEPeaks[fEPeaks.size()-1]/964.014;

        double bestrap = 1;

        for(int i=0 ; i<fNFound ; i++)
        {
            for(int j=i+1 ; j<fNFound ; j++)
            {
                double rap = TMath::Abs((fXPeaks[j]/fXPeaks[i])/rapmaxmin-1);

                if(rap<0.005 && fYPeaks[i]>0.9*fYPeaks[j] && rap<bestrap)
                {
                    if(fSourceName.Contains("Eu"))
                    {
                        for(int k=0 ; k<fNFound ; k++)
                        {
                            double rap2 = TMath::Abs((fXPeaks[j]/fXPeaks[k])/rapmaxmin2-1);
                            if(rap2<0.005)
                            {
                                min= i;
                                max= j;
                                bestrap=rap;
                            }
                        }
                    }
                    else
                    {
                        min= i;
                        max= j;
                        bestrap=rap;
                    }
                }
            }
        }
        if(bestrap==1)
        {
            Error("FitSpectra::Fit()","First and last peak not found ==> Calibration not feasible, exit Fit() method.");
            return;
        }
        else
        {
            //            cout<<" Peak at "<<fEPeaks[0]<<" kev found at canal "<<fXPeaks[min]<<endl;
            //            cout<<" Peak at "<<fEPeaks[fNPeaksToFit-1]<<" kev found at canal "<<fXPeaks[max]<<endl;
        }

        fXPeaks.erase(fXPeaks.begin()+max+1,fXPeaks.end());
        fYPeaks.erase(fYPeaks.begin()+max+1,fYPeaks.end());

        fXPeaks.erase(fXPeaks.begin(),fXPeaks.begin()+min);
        fYPeaks.erase(fYPeaks.begin(),fYPeaks.begin()+min);

        fNFound = max-min+1;

        //        cout<<fXPeaks.size()<<" "<<fNFound<<endl;
        //        cout<<fXPeaks[0]<<" "<<fXPeaks[fXPeaks.size()-1]<<endl;

        vector < int > BadPeaks;
        vector < int > PeakNotFound;

        int lastgoodpeak=0;

        for(int i=1 ; i<fNPeaksToFit-1 ; i++)
        {
            bestrap = 1;
            bool peakfound = false;

            double rapE = fEPeaks[i]/fEPeaks[0];

            for(int j=lastgoodpeak+1 ; j<fNFound-1 ; j++)
            {
                double rap = TMath::Abs((fXPeaks[j]/fXPeaks[0])/rapE-1);
                if(rap<0.007&&rap<bestrap)
                {
                    for(int k=lastgoodpeak+1 ; k<j ; k++) BadPeaks.push_back(k);
                    lastgoodpeak=j;
                    peakfound = true;
                    break;
                }
            }
            if(!peakfound)
            {
                //                cout<<" Peak at "<<fEPeaks[i]<<" kev not found, peak ignored"<<endl;
                PeakNotFound.push_back(i);
            }
            else
            {
                //                cout<<" Peak at "<<fEPeaks[i]<<" kev found at canal "<<fXPeaks[lastgoodpeak]<<endl;
            }
        }

        fXPeaks.erase(fXPeaks.begin()+lastgoodpeak+1,fXPeaks.end()-1);
        fYPeaks.erase(fYPeaks.begin()+lastgoodpeak+1,fYPeaks.end()-1);

        for(int i=BadPeaks.size()-1 ; i>=0 ; i--)
        {
            fXPeaks.erase(fXPeaks.begin()+BadPeaks[i]);
            fYPeaks.erase(fYPeaks.begin()+BadPeaks[i]);
            fNFound--;
        }
        for(int i=PeakNotFound.size()-1 ; i>=0 ; i--)
        {
            fEPeaks.erase(fEPeaks.begin()+PeakNotFound[i]);
            fNPeaksToFit--;
        }
    }

    double xmin = fXPeaks[0] - fXPeaks[fNPeaksToFit-1]*0.05;
    double xmax = fXPeaks[fNPeaksToFit-1] + fXPeaks[fNPeaksToFit-1]*0.05;
    fHistToFit->SetAxisRange(xmin,xmax);

    //    cout<<fNPeaksToFit<<" "<<fNFound<<" "<<xmin<<" "<<xmax<<endl;

    gPad->Update();
    gPad->Modified();

    if(fSourceName=="baseline")
    {
        xmin =  fXPeaks[0]*(1-0.1);
        xmax =  fXPeaks[fNFound-1]*(1+0.2);
    }

    int nbins = fHistToFit->FindBin(xmax)-fHistToFit->FindBin(xmin);
    CurrentPad->SetLogy();

    // through histograms since TSpectrum API has been modified at some point ... I am not sure the version is exactly the right one
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,0,0)
    double * source = new double[nbins];
#else
    float * source = new float[nbins];
#endif

    for(int i=0 ; i< nbins ; i++) source[i] = fHistToFit->GetBinContent(fHistToFit->FindBin(xmin)+i+1);

    if(fDefinePersoBackground)
    {
        s->Background(source,nbins,fNumberIterations,fDirection,fFilterOrder,fSmoothing,fsmoothingWindow,fCompton);
    }
    else if(fSourceName=="152Eu")
    {
        s->Background(source,nbins,30,TSpectrum::kBackDecreasingWindow,TSpectrum::kBackOrder2,true,TSpectrum::kBackSmoothing3,false);
    }
    else if(fSourceName=="baseline")
    {
        s->Background(source,nbins,20,TSpectrum::kBackDecreasingWindow,TSpectrum::kBackOrder2,true,TSpectrum::kBackSmoothing3,false);
    }
    else if(fHistName.Contains("36") || fHistName.Contains("37"))
    {
        s->Background(source,nbins,90,TSpectrum::kBackDecreasingWindow,TSpectrum::kBackOrder2,true,TSpectrum::kBackSmoothing3,false);
    }
    else
    {
        s->Background(source,nbins,300,TSpectrum::kBackDecreasingWindow,TSpectrum::kBackOrder8,true,TSpectrum::kBackSmoothing3,true);
    }

    if(gROOT->FindObject(fBackGroundName) != 0x0)
    {
        TH1 *h = (TH1*) gROOT->FindObject(fBackGroundName);
        delete h;
    }

    BackGroundHist = new TH1F(fBackGroundName,fBackGroundName,nbins,xmin,xmax);

    for (int i = 0; i < nbins; i++) BackGroundHist->SetBinContent(i+1,source[i]);

    if(gROOT->FindObject(fFitFctName) != 0x0)
    {
        TH1 *f = (TH1*) gROOT->FindObject(fFitFctName);
        delete f;
    }

    int NPar = 0;

    if(fFitType =="Gaus")
    {
        NPar = 3;

        fFitFct = new TF1(fFitFctName,this,&FitSpectra::BGPlusGaus,xmin,xmax,fNPeaksToFit*NPar,"FitSpectra","BGPlusGaus");
        fFitFct->SetNpx(nbins);

        for(int i=0 ; i<fNPeaksToFit ; i++)
        {
            fFitFct->SetParameter(0+NPar*i,fYPeaks[i]-BackGroundHist->GetBinContent(BackGroundHist->FindBin(fXPeaks[i])));
            fFitFct->SetParLimits(0+NPar*i,fFitFct->GetParameter(0+NPar*i)/2.,fFitFct->GetParameter(0+NPar*i)*2.);

            fFitFct->SetParameter(1+NPar*i,fXPeaks[i]);
            fFitFct->SetParLimits(1+NPar*i,fXPeaks[i]*(1-2e-3),fXPeaks[i]*(1+2e-3));

            fFitFct->SetParameter(2+NPar*i,fXPeaks[i]*5e-4);
            fFitFct->SetParLimits(2+NPar*i,0.,20.);
        }
    }
    else if(fFitType =="CBDS")
    {
        NPar = 7;

        fFitFct = new TF1(fFitFctName,this,&FitSpectra::BGPlusCBDS,xmin,xmax,fNPeaksToFit*NPar,"FitSpectra","BGPlusCBDS");
        fFitFct->SetNpx(nbins);

        for(int i=0 ; i<fNPeaksToFit ; i++)
        {
            fFitFct->SetParameter(0+NPar*i,fYPeaks[i]-BackGroundHist->GetBinContent(BackGroundHist->FindBin(fXPeaks[i])));
            fFitFct->SetParLimits(0+NPar*i,fFitFct->GetParameter(0+NPar*i)*0.95,fFitFct->GetParameter(0+NPar*i)*1.05);

            fFitFct->SetParameter(1+NPar*i,fXPeaks[i]);
            fFitFct->SetParLimits(1+NPar*i,fXPeaks[i]*(1-0.001),fXPeaks[i]*(1+0.001));

            fFitFct->SetParameter(2+NPar*i,fXPeaks[i]*5e-4);
            //fFitFct->SetParLimits(2+NPar*i,0.,fXPeaks[i]*1e-2);

            fFitFct->SetParameter(3+NPar*i,2.);
            fFitFct->SetParLimits(3+NPar*i,0.,5.);

            fFitFct->SetParameter(4+NPar*i,5.);
            fFitFct->SetParLimits(4+NPar*i,0.,25.);

            fFitFct->SetParameter(5+NPar*i,2.);
            fFitFct->SetParLimits(5+NPar*i,0.,5.);

            fFitFct->SetParameter(6+NPar*i,5.);
            fFitFct->SetParLimits(6+NPar*i,0.,25.);
        }
    }
    else if(fFitType =="Dino")
    {
        NPar = 3;
        double par[fNPeaksToFit*NPar];

        for(int i=0 ; i<fNPeaksToFit ; i++)
        {
            TF1 *Gaus = new TF1("Gaus",this,&FitSpectra::BGPlusGaus,fXPeaks[i]-5.,fXPeaks[i]+10.,NPar,"FitSpectra","BGPlusGaus");

            Gaus->SetParameter(0,fYPeaks[i]-BackGroundHist->GetBinContent(BackGroundHist->FindBin(fXPeaks[i])));

            Gaus->SetParameter(1,fXPeaks[i]);
            Gaus->SetParLimits(1,fXPeaks[i]-1.,fXPeaks[i]+1.);

            Gaus->SetParameter(2,fXPeaks[i]*1e-3);
            Gaus->SetParLimits(2,0.,fXPeaks[i]*2e-3);

            fHistToFit->Fit(Gaus,"QRS0");

            par[0+NPar*i] = Gaus->GetParameter(0);
            par[1+NPar*i] = Gaus->GetParameter(1);
            par[2+NPar*i] = Gaus->GetParameter(2);

            delete Gaus;
        }
        NPar = 6;

        fFitFct = new TF1(fFitFctName,this,&FitSpectra::BGPlusDinoFct,xmin,xmax,fNPeaksToFit*NPar,"FitSpectra","BGPlusDinoFct");
        fFitFct->SetNpx(nbins);

        for(int i=0 ; i<fNPeaksToFit ; i++)
        {
            fFitFct->FixParameter(0+NPar*i,par[0+3*i]);

            fFitFct->FixParameter(1+NPar*i,par[1+3*i]);

            fFitFct->SetParameter(2+NPar*i,par[2+3*i]);
            fFitFct->SetParLimits(2+NPar*i,1./sqrt(12.),par[2+3*i]*2.35);

            fFitFct->SetParameter(3+NPar*i,-2.);
            fFitFct->SetParLimits(3+NPar*i,-50,-0.1);

            fFitFct->SetParameter(4+NPar*i,2.);
            fFitFct->SetParLimits(4+NPar*i,0.1,50.);

            fFitFct->FixParameter(5+NPar*i,0.);
        }
    }
    fHistToFit->Fit(fFitFct,"QRS0+");
    TF1 *FitFctClone = (TF1*)fFitFct->Clone();


    for(int i=0 ; i<fNPeaksToFit ; i++)
    {
        TF1 *FNoBG = 0x0;
        if(fFitType =="Gaus") FNoBG = new TF1("FNoBG",this,&FitSpectra::MyGaus,xmin,xmax,NPar,"FitSpectra","MyGaus");
        else if(fFitType =="CBDS") FNoBG = new TF1("FNoBG",this,&FitSpectra::fnc_dscb,xmin,xmax,NPar,"FitSpectra","fnc_dscb");
        else if(fFitType =="Dino") FNoBG = new TF1("FNoBG",this,&FitSpectra::DinoFct,xmin,xmax,NPar,"FitSpectra","DinoFct");

        for(int j=0 ; j<NPar ; j++) FNoBG->SetParameter(j,fFitFct->GetParameter(j+NPar*i));

        double mean = FNoBG->GetParameter(1);
        double sigma = FNoBG->GetParameter(2);

        fMean.push_back(mean);
        fSigma.push_back(sigma);

        double N = FNoBG->GetParameter(0);

        double fwhm = FNoBG->GetX(N/2.,mean,mean+5*sigma)-FNoBG->GetX(N/2.,mean-5*sigma,mean);

        fFWHM.push_back(fwhm);

        double fwhm_10 = FNoBG->GetX(N/10.,mean,mean+5*sigma)-FNoBG->GetX(N/10.,mean-5*sigma,mean);

        fFWHM_10.push_back(fwhm_10);

        double LTProp = (mean-FNoBG->GetX(N/10.,mean-5*sigma,mean))/(mean-FNoBG->GetX(N/2.,mean-5*sigma,mean));

        fLeftTailProp.push_back(LTProp);

        double RTProp = (mean-FNoBG->GetX(N/10.,mean,mean+5*sigma))/(mean-FNoBG->GetX(N/2.,mean,mean+5*sigma));

        fRightTailProp.push_back(RTProp);

        double area = FNoBG->Integral(mean-5*sigma,mean+5*sigma);
        fArea.push_back(area);

        delete FNoBG;
    }

    if(fNPeaksToFit>1)
    {
        FitStatus = 1;
        Calib();
    }


    FitFctClone->Draw("same");

    TH1F *BackGroundHistClone = (TH1F*)BackGroundHist->Clone();
    BackGroundHistClone->SetLineColor(kGreen);
    BackGroundHistClone->SetLineStyle(7);
    BackGroundHistClone->Draw("same");

    fHistToFit->SetAxisRange(xmin,xmax);

    CurrentPad->SetEditable(true);

    gPad->Update();
    gPad->Modified();

    gErrorIgnoreLevel = kPrint;

    delete [] source;
    delete s;
}

void FitSpectra::Calib()
{
    const int n = fNPeaksToFit;

    double Raw[n];
    double E[n];

    for(int i=0 ; i<fNPeaksToFit ; i++)
    {
        Raw[i] = fMean[i];
        E[i] = fEPeaks[i];
    }

    TGraph *g = new TGraph(fNPeaksToFit,Raw,E);

    TF1 *fpol = new TF1("f","[0]+[1]*x",0.,Raw[fNPeaksToFit-1]*1.2);
    if(!fUseOffset) fpol->FixParameter(0,0.);

    g->Fit(fpol,"RQ0");

    fCalibOffset = fpol->GetParameter(0);
    fCalibSlope = fpol->GetParameter(1);

    for(int i=0 ; i<fNPeaksToFit ; i++)
    {
        fEcal.push_back(fpol->Eval(fMean[i]));
        fFWHMCal.push_back(fCalibSlope*fFWHM[i]);
        fFWHM_10Cal.push_back(fCalibSlope*fFWHM_10[i]);
    }

    delete fpol;
}

double FitSpectra::fnc_dscb(double*xx,double*pp)
{
    double x   = xx[0];
    // gaussian core
    double N   = pp[0];//norm
    double mu  = pp[1];//mean
    double sig = pp[2];//variance
    // transition parameters
    double a1  = pp[3];
    double p1  = pp[4];
    double a2  = pp[5];
    double p2  = pp[6];

    double u   = (x-mu)/sig;
    double A1  = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2);
    double A2  = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2);
    double B1  = p1/TMath::Abs(a1) - TMath::Abs(a1);
    double B2  = p2/TMath::Abs(a2) - TMath::Abs(a2);

    double result(N);
    if      (u<-a1) result *= A1*TMath::Power(B1-u,-p1);
    else if (u<a2)  result *= TMath::Exp(-u*u/2);
    else            result *= A2*TMath::Power(B2+u,-p2);
    return result;
}

double FitSpectra::DinoFct(double*xx,double*pp)
{
    double x   = xx[0];
    // gaussian core
    double A   = pp[0];//norm
    double P  = pp[1];//mean
    double sigma = pp[2];//variance
    // transition parameters
    double lambda  = pp[3];
    double rho  = pp[4];
    double S  = pp[5];

    double u   = (x-P)/sigma;

    double f_g = A*TMath::Exp(-u*u*0.5);

    double f_lambda = A*TMath::Exp(-0.5*lambda*(2.*u-lambda));

    double f_rho = A*TMath::Exp(-0.5*rho*(2.*u-rho));

    double f_S = A*S*1./((1+TMath::Exp(u))*(1+TMath::Exp(u)));

    double f_tot = 0.;

    if(u<lambda) f_tot += f_lambda;
    else if(u>rho) f_tot += f_rho;
    else f_tot += f_g;

    f_tot += f_S;

    return f_tot;
}

double FitSpectra::MyGaus(double *xx, double *par)
{
    double x = xx[0];
    double N = par[0];
    double mean = par[1];
    double sigma = par[2];

    double f = N*exp(-0.5*((x-mean)/sigma)*((x-mean)/sigma));

    return f;
}

double FitSpectra::BGPlusGaus(double *xx, double *par)
{
    double x = xx[0];

    double f = BackGroundHist->GetBinContent(BackGroundHist->FindBin(x));

    for(int i=0 ; i<fNFound ; i++)
    {
        double *p = new double[3];

        p[0] = par[3*i+0];
        p[1] = par[3*i+1];
        p[2] = par[3*i+2];

        f += MyGaus(xx,p);
        delete p;
    }

    return f;
}

double FitSpectra::BGPlusCBDS(double *xx, double *par)
{
    double x = xx[0];

    double f = BackGroundHist->GetBinContent(BackGroundHist->FindBin(x));

    for(int i=0 ; i<fNFound ; i++)
    {
        const int size = 7;

        double *p = new double[size];

        p[0] = par[size*i+0];
        p[1] = par[size*i+1];
        p[2] = par[size*i+2];
        p[3] = par[size*i+3];
        p[4] = par[size*i+4];
        p[5] = par[size*i+5];
        p[6] = par[size*i+6];

        f += fnc_dscb(xx,p);
        delete p;
    }

    return f;
}

double FitSpectra::BGPlusDinoFct(double *xx, double *par)
{
    double x = xx[0];

    double f = BackGroundHist->GetBinContent(BackGroundHist->FindBin(x));

    for(int i=0 ; i<fNFound ; i++)
    {
        const int size = 6;

        double *p = new double[size];

        p[0] = par[size*i+0];
        p[1] = par[size*i+1];
        p[2] = par[size*i+2];
        p[3] = par[size*i+3];
        p[4] = par[size*i+4];
        p[5] = par[size*i+5];

        f += DinoFct(xx,p);
        delete p;
    }

    return f;
}

void FitSpectra::SetPersoBackground(int NumberIterations, int Direction,int FilterOrder,bool Smoothing,int smoothingWindow,bool Compton)
{
    fDefinePersoBackground = true;
    fNumberIterations = NumberIterations;
    fDirection = Direction;
    fFilterOrder = FilterOrder;
    fSmoothing = Smoothing;
    fsmoothingWindow = smoothingWindow;
    fCompton = Compton;
}

void FitSpectra::ReduceListOfPeak()
{
    double FirstEToKeep = 344.272;
    int Pos=0;
    for(int i=0 ; i<fNPeaksToFit ; i++)
    {
        if(fEPeaks[i] < FirstEToKeep) Pos = i;
    }

    fEPeaks.erase(fEPeaks.begin(),fEPeaks.begin()+Pos+1);
    fNPeaksToFit = fEPeaks.size();
}