//dEdx-Fit
//Author: P. Gadow
//status: version 0.8

#include "TROOT.h"
#include "TCanvas.h"
#include "TGraphErrors.h"
#include "TAxis.h"
#include "TFrame.h"
#include "TF1.h"
#include "TLegend.h"
#include "TArrow.h"
#include "TLatex.h"
#include "TFile.h"
#include "TStyle.h"
#include "TString.h"
#include "TH1F.h"
#include "TMath.h"
#include "TPaveText.h"

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


std::string inttostring(Int_t input)
{
    std::string s;
    std::stringstream out;
    out << input;
    s = out.str();
    return s;
}

std::string doubletostring(double input)
{
    std::string s;
    std::stringstream out;
    out << input;
    s = out.str();
    return s;
}

std::string round1(double input)
{
    std::string s;
    std::stringstream out;
    out << std::fixed << std::setprecision(1) << input;
    s = out.str();
    return s;
}
std::string round3(double input)
{
    std::string s;
    std::stringstream out;
    out << std::fixed << std::setprecision(3) << input;
    s = out.str();
    return s;
}


//###################################################################
//$$$$$$$$$$$$$$$$$$$   Fitting Functions   $$$$$$$$$$$$$$$$$$$$$$$$$
//###################################################################

Double_t Gauss(Double_t *x, Double_t *par)
{
    // A single gaussian

    Double_t arg = 0;
    if (par[2]) arg = (x[0] - par[1])/par[2];

    Double_t sig = par[0]*TMath::Exp(-0.5*arg*arg);
    return sig;
}


Double_t Signal(Double_t *x, Double_t *par)
{
    // The signal function: two gaussians

    Double_t arg1 = 0;
    if (par[2]) arg1 = (x[0] - par[1])/par[2];

    Double_t arg2 = 0;
    if (par[5]) arg2 = (x[0] - par[4])/par[5];

    Double_t sig =
            par[0]*TMath::Exp(-0.5*arg1*arg1) + par[3]*TMath::Exp(-0.5*arg2*arg2);
    return sig;
}

Double_t Background(Double_t *x, Double_t *par)
{
    // The background function: p3

    Double_t back =
            par[0] + par[1]*x[0] + par[2]*x[0]*x[0] + par[3]*x[0]*x[0]*x[0];
    return back;
}

Double_t Signal_Background(Double_t *x, Double_t *par)
{
    // The whole function: Signal + Background

    Double_t arg1 = 0;
    if (par[2]) arg1 = (x[0] - par[1])/par[2];

    Double_t arg2 = 0;
    if (par[5]) arg2 = (x[0] - par[4])/par[5];

    Double_t sig =
            par[0]*TMath::Exp(-0.5*arg1*arg1) + par[3]*TMath::Exp(-0.5*arg2*arg2)
            + par[6] + par[7]*x[0] + par[8]*x[0]*x[0] + par[9]*x[0]*x[0]*x[0] ;
    return sig;
}

Double_t Gauss_Background(Double_t *x, Double_t *par)
{
    // Two functions: Gauss + Background

    Double_t arg1 = 0;
    if (par[2]) arg1 = (x[0] - par[1])/par[2];


    Double_t sig =
            par[0]*TMath::Exp(-0.5*arg1*arg1)
            + par[3] + par[4]*x[0] + par[5]*x[0]*x[0] + par[6]*x[0]*x[0]*x[0] ;
    return sig;
}

//###################################################################
//$$$$$$$$$$$$$$$$$$$$$$$$$$$   Fitting  $$$$$$$$$$$$$$$$$$$$$$$$$$$$
//###################################################################


void FindParametersDoubleGauss(TH1* hist1, Double_t* para, int meanpion, int meanelectron, int sigmapion, int sigmaelectron, Double_t scale)
{
    // Find start values for fits, scale is #ADC channels per bin


    //Pions
    Double_t maxpi = 0;
    int meanpi = 0;
    int hwhmpi = 0;
    int backstart = 0;

    Int_t nbins = hist1->GetNbinsX();

    // Take Maximum
    meanpi = meanpion;
    maxpi = hist1->GetBinContent(meanpi);


//    // Find HWHM Pions
//    for (int i = meanpi; i < nbins; i++){
//        if (hist1->GetBinContent(i) < maxpi/2){
//            hwhmpi = i - meanpi;
//            break;
//        }
//    }

    hwhmpi = sigmapion;

    // Find beginning of Spectrum
    for (int i = 0; i < nbins; i++){
        if (hist1->GetBinContent(i) > float(maxpi/50.0)){
            backstart = i;
            break;
        }
    }

    //Electrons
    Double_t maxel = 0;
    int meanel = 0;
    int hwhmel = 0;


    // Take Maximum
    meanel = meanelectron;
    maxel = hist1->GetBinContent(meanel);


//    // Find HWHM Electrons
//    for (int i = meanel; i < nbins; i++){
//        if (hist1->GetBinContent(i) < maxel/2){
//            hwhmel = i - meanel;
//            break;
//        }
//    }

    hwhmel = sigmaelectron;

    std::cout << "\nFindParameters done... \nnBins: "<<nbins;
    std::cout << "\nFit Parameters for Pion Peak: Max: "<< maxpi
            <<" Mean: " << meanpi*scale << " HWHM: "<<hwhmpi*scale
            << " Signal: "<<backstart*scale<< "\n";
    std::cout << "\nFit Parameters for Electron Peak: Max: "<< maxel
            <<" Mean: " << meanel*scale << " HWHM: "<<hwhmel*scale
            << " Signal: "<<backstart*scale<< "\n";

    // Parameters for DoubleGauss Fit
    para[0] = maxpi;
    para[1] = meanpi*scale;
    para[2] = hwhmpi*scale;
    para[3] = maxel;
    para[4] = meanel*scale;
    para[5] = hwhmel*scale;
    para[6] = backstart*scale;

}

//------------------------------------------------------------------------

void FindParametersSingleGauss(TH1* hist1, Double_t* para, int *exportmean,int *exporthwhm, Double_t scale)
{
    // Find start values for fits, scale is #ADC channels per bin

    Double_t max = 0;
    int mean = 0;
    int hwhm = 0;
    int backstart = 0;

    Int_t nbins = hist1->GetNbinsX();

    // Find Maximum
    for (int i = int(nbins/10); i < nbins; i++){
        if (hist1->GetBinContent(i) > max) {
            max = hist1->GetBinContent(i);
            mean = i;
        }
    }

    // Find HWHM
    for (int i = mean; i < nbins; i++){
        if (hist1->GetBinContent(i) < max/2){
            hwhm = i - mean;
            break;
        }
    }

    // Find beginning of Spectrum
    for (int i = 0; i < nbins; i++){
        if (hist1->GetBinContent(i) > float(max/50.0)){
            backstart = i;
            break;
        }
    }

    std::cout << "\nFindParameters done... \nnBins: "<<nbins;
    std::cout << "\nFit Parameters: Max: "<< max
            <<" Mean: " << mean*scale << " HWHM: "<<hwhm*scale
            << " Signal: "<<backstart*scale<< "\n";

    // Parameters for SingleGauss Fit
    para[0] = max;
    para[1] = mean*scale;
    para[2] = hwhm*scale;
    para[3] = backstart*scale;

    *exportmean = mean;
    *exporthwhm = hwhm;
}

//------------------------------------------------------------------------

void FitHistogramSingleGauss (TH1* hist1, Double_t* para, Double_t* paraerr, Double_t scale)
{
    Double_t min = 0;
    Double_t max = 0;
    Double_t startfit = 0;
    Double_t endfit = 0;
    Double_t backstart = 0;

    gROOT->Reset();

    max = hist1->GetNbinsX();
    backstart = para[3];


    // Fit Signal in subrange
    //=======================

    // Get Fit Range
    startfit = para[1] - 2*para[2];
    endfit = para[1] + 2*para[2];

    std::cout << startfit << " single " << endfit << std::endl;
    // Create Fitting Function
    TF1 *func1 = new TF1("Gauss",Gauss,startfit,endfit,3);
    func1->SetNpx(1000);      // Draw 1000 Points
    func1->SetLineColor(4);   // Green
    func1->SetLineWidth(1);
    func1->SetParameters(para[0],para[1],para[2]);
    func1->SetParLimits(1,0,500000);
    func1->SetParLimits(1,0,1000);
    func1->SetParLimits(2,0,300);

    func1->SetParNames("Amp","Mean","#sigma");

    hist1->Fit("Gauss", "RN");  // Fit without drawing


    // Fit Background
    //===============


    // Create Histogram
    TH1F *hist2 = (TH1F*)hist1->Clone();
    hist2->Reset();
    for (int i = 1;i<=max;i++) {
        Float_t x = hist1->GetBinCenter(i);
        Double_t fval = func1->Eval(x);
        Double_t diff = hist1->GetBinContent(i)-fval;
        if (diff < 0) diff = 0;
        hist2->Fill(x,diff);
    }

    // Create Fitting Function
    TF1 *func2 = new TF1("Background",Background,backstart,max*scale,4);
    func2->SetNpx(1000);      // Draw 1000 Points
    func2->SetLineColor(6);   // Magenta
    func2->SetLineWidth(1);

    hist2->Fit("Background", "RN");



    // Fit Signal+Background
    //=======================


    // Read back Fit values
    func1->GetParameters(&para[0]);
    func2->GetParameters(&para[3]);


    // Create Fitting Function
    TF1 *func3 = new TF1("Gauss_Background",Gauss_Background,
                         backstart,max*scale,7);
    func3->SetNpx(1000);      // Draw 1000 Points
    func3->SetLineColor(2);   // Red
    func3->SetLineWidth(2);
    func3->SetParameters(para[0],para[1],para[2],
                         para[3],para[4],para[5],
                         para[6]);
    func3->SetParNames("Amp","Mean","#sigma",
                       "b0", "b1", "b2", "b3");
    func3->SetParLimits(1,0,1000);
    func3->SetParLimits(2,0,300);

    hist1->Fit("Gauss_Background", "RN");

    func3->SetRange(backstart, max*scale);
    func3->GetParameters(&para[0]);

    //retrieve parameter errors
    paraerr[0] = func3->GetParError(0); paraerr[1] = func3->GetParError(1); paraerr[2] = func3->GetParError(2); paraerr[3] = func3->GetParError(3); paraerr[4] = func3->GetParError(4); paraerr[5] = func3->GetParError(5); paraerr[6] = func3->GetParError(6);

    TF1 *func2_final = new TF1("Background_final",
                               Background,backstart,max*scale,4);
    func2_final->SetLineColor(6);   // Magenta
    func2_final->SetLineWidth(1);
    func2_final->SetParameters(func3->GetParameter(3),
                               func3->GetParameter(4),
                               func3->GetParameter(5),
                               func3->GetParameter(6));

    // Drawing Histogram + Functions, draw all into the same Frame
    hist1->Draw();
    //func1->DrawCopy("lsame");   // use DrawCopy instead of Draw because
    //funcg1->DrawCopy("lsame");
    //func2->DrawCopy("lsame");   // drawing was disabled for Fit
    //funcg2->DrawCopy("lsame");
    func2_final->DrawCopy("lsame");
    func3->DrawCopy("lsame");
}

//------------------------------------------------------------------------

void FitHistogramDoubleGauss (TH1* hist1, Double_t* para, Double_t* paraerr, Double_t scale)
{
  Double_t min = 0;
  Double_t max = 0;
  Double_t startfit = 0;
  Double_t endfit = 0;
  Double_t backstart = 0;

  gROOT->Reset();

  max = hist1->GetNbinsX();
  backstart = para[6];


  // Fit Signal in subrange
  //=======================

  // Get Fit Range
  startfit = para[1] - 2*para[2];
  endfit = para[4] + 2*para[5];

  std::cout << startfit << " double " << endfit << std::endl;

  // Create Fitting Function
  TF1 *func1 = new TF1("Signal",Signal,startfit,endfit,6);
  func1->SetNpx(1000);      // Draw 1000 Points
  func1->SetLineColor(4);   // Green
  func1->SetLineWidth(1);
  func1->SetParameters(para[0],para[1],para[2],para[3],para[4],para[5]);
  func1->SetParNames("Amp1","Mean1","#sigma1","Amp2","Mean2","#sigma2");
  func1->SetParLimits(1,0,1000);
  func1->SetParLimits(2,0,300);
  func1->SetParLimits(4,0,1000);
  func1->SetParLimits(5,0,300);
  hist1->Fit("Signal", "RN");  // Fit without drawing

  TF1 *funcg1 = new TF1("g1",Gauss,startfit,endfit,3);
  funcg1->SetLineColor(3);   // Green
  funcg1->SetLineWidth(1);
  funcg1->SetParameters(func1->GetParameter(0),
                        func1->GetParameter(1),
                        func1->GetParameter(2));
  funcg1->SetParLimits(1,0,1000);
  funcg1->SetParLimits(2,0,300);

  TF1 *funcg2 = new TF1("g2",Gauss,startfit,endfit,3);
  funcg2->SetLineColor(3);   // Green
  funcg2->SetLineWidth(1);
  funcg2->SetParameters(func1->GetParameter(3),
                        func1->GetParameter(4),
                        func1->GetParameter(5));
  funcg2->SetParLimits(1,0,1000);
  funcg2->SetParLimits(2,0,300);

  // Fit Background
  //===============


  // Create Histogram
  TH1F *hist2 = (TH1F*)hist1->Clone();
  hist2->Reset();
  for (int i = 1;i<=max;i++) {
    Float_t x = hist1->GetBinCenter(i);
    Double_t fval = func1->Eval(x);
    Double_t diff = hist1->GetBinContent(i)-fval;
    if (diff < 0) diff = 0;
    hist2->Fill(x,diff);
  }

  // Create Fitting Function
  TF1 *func2 = new TF1("Background",Background,backstart,max*scale,4);
  func2->SetNpx(1000);      // Draw 1000 Points
  func2->SetLineColor(6);   // Magenta
  func2->SetLineWidth(1);

  hist2->Fit("Background", "RN");


  // Fit Signal+Background
  //=======================


  // Read back Fit values
  func1->GetParameters(&para[0]);
  func2->GetParameters(&para[6]);


  // Create Fitting Function
  TF1 *func3 = new TF1("Signal_Background",Signal_Background,
                       backstart,max*scale,10);
  func3->SetNpx(1000);      // Draw 1000 Points
  func3->SetLineColor(2);   // Red
  func3->SetLineWidth(2);
  func3->SetParameters(para[0],para[1],para[2],
                       para[3],para[4],para[5],
                       para[6],para[7], para[8], para[9]);
//  func3->FixParameter(1,para[1]);
//  func3->FixParameter(2,para[2]);
//  func3->FixParameter(4,para[4]);
//  func3->FixParameter(5,para[5]);
  func3->SetParNames("Amp1","Mean1","#sigma1",
                     "Amp2","Mean2","#sigma2",
                     "b0", "b1", "b2", "b3");
  func3->SetParLimits(1,0,1000);
  func3->SetParLimits(2,0,300);
  func3->SetParLimits(4,0,1000);
  func3->SetParLimits(5,0,300);
  hist1->Fit("Signal_Background", "RN");

  func3->SetRange(backstart, max*scale);
  func3->GetParameters(&para[0]);

  //retrieve parameter errors
  paraerr[0] = func3->GetParError(0); paraerr[1] = func3->GetParError(1); paraerr[2] = func3->GetParError(2); paraerr[3] = func3->GetParError(3); paraerr[4] = func3->GetParError(4); paraerr[5] = func3->GetParError(5); paraerr[6] = func3->GetParError(6); paraerr[7] = func3->GetParError(7); paraerr[8] = func3->GetParError(8); paraerr[9] = func3->GetParError(9);


  TF1 *funcg1_final = new TF1("g1_final",Gauss,backstart,max*scale,3);
  funcg1_final->SetLineColor(3);   // Green
  funcg1_final->SetLineWidth(1);
  funcg1_final->SetParameters(func3->GetParameter(0),
                              func3->GetParameter(1),
                              func3->GetParameter(2));

  TF1 *funcg2_final = new TF1("g2_final",Gauss,backstart,max*scale,3);
  funcg2_final->SetLineColor(3);   // Green
  funcg2_final->SetLineWidth(1);
  funcg2_final->SetParameters(func3->GetParameter(3),
                              func3->GetParameter(4),
                              func3->GetParameter(5));

  TF1 *func2_final = new TF1("Background_final",
                             Background,backstart,max*scale,4);
  func2_final->SetLineColor(6);   // Magenta
  func2_final->SetLineWidth(1);
  func2_final->SetParameters(func3->GetParameter(6),
                              func3->GetParameter(7),
                              func3->GetParameter(8),
                              func3->GetParameter(9));

  // Drawing Histogram + Functions, draw all into the same Frame
  hist1->Draw();
  //func1->DrawCopy("lsame");   // use DrawCopy instead of Draw because
  //funcg1->DrawCopy("lsame");
  //func2->DrawCopy("lsame");   // drawing was disabled for Fit
  //funcg2->DrawCopy("lsame");
  funcg1_final->DrawCopy("lsame");
  funcg2_final->DrawCopy("lsame");
  func2_final->DrawCopy("lsame");
  func3->DrawCopy("lsame");
  }

//------------------------------------------------------------------------

int GetSeparationPower (TH1* hist1, double *separationpower, int meanpion, int meanelectron, int hwhmpion, int hwhmelectron, Double_t scale)
{
    // Declare data variables
    Double_t s[15];
    Double_t serr[15];


    for (int i =0; i<15; i++) s[i] = 0;
    for (int i =0; i<15; i++) serr[i] = 0;


   FindParametersDoubleGauss(hist1, s, meanpion, meanelectron, hwhmpion, hwhmelectron, 1.);
//    para[0] = max;
//    para[1] = mean*scale;
//    para[2] = hwhm*scale;
//    para[3] = max;//0.1 * max;
//    para[4] = mean*scale; // 0.673*mean*scale;
//    para[5] = hwhm*scale;
//    para[6] = backstart*scale;
    FitHistogramDoubleGauss(hist1, s,&serr[0], scale);
//    "Amp1","Mean1","#sigma1",
//    "Amp2","Mean2","#sigma2",
//    "b0", "b1", "b2", "b3"

    // Extract Separation Power

    *separationpower = 2*((s[4]-s[1]))/(s[2]+ s[5]);

    return 0;

}



//###################################################################
//$$$$$$$$$$$$$$$$$$$$$$$$    Start Macro   $$$$$$$$$$$$$$$$$$$$$$$$$
//###################################################################


void ExecutedEdxfit(bool verbose, int runnumber, int datasetnumber, int lowtrunc, int hightrunc, bool pretty = false){ //TString histofile

    if (pretty){
        gROOT->ProcessLine(".x ~/rootlogon_Bernhard.C");
        gROOT->SetStyle("col");
        gROOT->ForceStyle();

    }

    TString indirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdx/";
    //TString indirname = "/nfs/hicran/project/panda/SIM/pgadow/fopiROOT_trunk/macro/ALICE/plots/dEdx/truncations/";

    TString histofile = "merge_" + inttostring(runnumber) + "_"+inttostring(datasetnumber)+"_dEdx_truncmean_" +inttostring(lowtrunc)+"_" +inttostring(hightrunc)+"histograms.root";
    TString inpathname = indirname + histofile;


    // create output file
    ofstream outputfile;
    outputfile.open ("/nfs/hicran/project/panda/SIM/pgadow/fopiROOT_trunk/macro/ALICE/plots/dEdx/FITDATA.dat", std::ofstream::app);

    TString dirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdx/fit/";
    TString filename = "merge_" + inttostring(runnumber) + "_"+ inttostring(datasetnumber) + "_dEdx";
    filename+= "_truncmean_"+inttostring(lowtrunc)+"_"+inttostring(hightrunc);


    TString savefilename = dirname + filename + "_fit.root";
//    TFile *savefile = new TFile(savefilename.Data(),"UPDATE");


    // open file
    TFile *histograms = new TFile(inpathname);

    TH1F *histall = (TH1F *)histograms->Get("histall");
    TH1F *histPID = (TH1F *)histograms->Get("histPID");
    TH1F *histpion = (TH1F *)histograms->Get("histpion");
    TH1F *histelectron = (TH1F *)histograms->Get("histEsinglecut");

    histPID->SetTitle("");
    histPID->GetXaxis()->SetRangeUser(0,400);
    // create canvas
    TCanvas *cdEdx = new TCanvas("cdEdxFit","",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->Divide(2,2);
    cdEdx->SetGrid();
    TCanvas *cdEdxbothpeaks = new TCanvas("cdEdxFitSpectrum","",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->SetGrid();
    TCanvas *cdEdxelectron = new TCanvas("cdEdxFitElectron","",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->SetGrid();
    TCanvas *cdEdxpion = new TCanvas("cdEdxFitPion","",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->SetGrid();

    gStyle->SetOptStat(0);      // switch off Statistics Display

    if (verbose){

        cout << "\n\n================================\n";
        cout <<     "Starting Analysis \n";
        cout <<     "================================\n\n";
    }



    // Analyze Spectrum
    if (verbose){
        cout << "Analyzing Spectrum... \n \n";
    }

    //------------------------------------------------------------
    //Pion fit

//    cdEdx->cd(3);
    cdEdxpion->cd();

    Double_t parpionfit[15];
    for (int i =0; i<15; i++) parpionfit[i] = 0;
    Double_t parpionfiterrors[15];
    for (int i =0; i<15; i++) parpionfiterrors[i] = 0;

    int pionmean = 0;
    int pionhwhm = 0;

    double scalepi = 1;

    FindParametersSingleGauss(histpion, parpionfit, &pionmean,&pionhwhm, scalepi);
    if (verbose){
        std::cout<<"before pion fit: "<<parpionfit[0]<<" "<<parpionfit[1]<<" "<<parpionfit[2]<<" "<<std::endl;
    }
    FitHistogramSingleGauss(histpion, parpionfit, &parpionfiterrors[0], scalepi);

    if (verbose){
        std::cout<<"after pion fit: "<<parpionfit[0]<<" "<<parpionfit[1]<<" "<<parpionfit[2]<<" "<<std::endl;
        std::cout<<"errors of pion fit: "<<parpionfiterrors[0]<<" "<<parpionfiterrors[1]<<" "<<parpionfiterrors[2]<<" "<<std::endl;
    }



    // Create Lables with Fit results

    TPaveText *ptX3 = new TPaveText(0.7, 0.71, 0.95, 0.87, "NDC"); //    TPaveText *ptX3 = new TPaveText(0.75, 0.76, 0.95, 0.87, "NDC");

    ptX3->SetFillColor(10);

    Double_t meanpi = parpionfit[1];
    Double_t sigmapi = parpionfit[2]/parpionfit[1];
    Double_t dERespi = sigmapi*2.3548;

    Double_t meanpierr = parpionfiterrors[1];
    Double_t sigmapierr = parpionfit[2]/parpionfit[1] * TMath::Sqrt(TMath::Power((parpionfiterrors[2]/parpionfit[2]),2)+TMath::Power((parpionfiterrors[1]/parpionfit[1]),2));
    Double_t dERespierr = sigmapierr*2.3548;

//    char valX31[10] = "";
//    char errX31[10] = "";
//    char txtX31[25] = "Fitted Mean #mu = ";
//    char valX32[10] = "";
//    char errX32[10] = "";
//    char txtX32[25] = "#sigma/#mu [%] = ";
//    char valX33[10] = "";
//    char errX33[10] = "";
//    char txtX33[25] = "#Delta E/E [%] = ";
//    int dummy3 = 0;

//    dummy3 = sprintf(valX31, "%.3f", double(meanpi)*scalepi);
//    dummy3 = sprintf(valX32, "%.3f", 100*sigmapi);
//    dummy3 = sprintf(valX33, "%.3f", 100*dERespi);
//    strcat(txtX31, valX31);
//    strcat(txtX32, valX32);
//    strcat(txtX33, valX33);
//
//    dummy3 = sprintf(errX31, "#pm %.3f", double(meanpierr)*scalepi);
//    dummy3 = sprintf(errX32, "#pm %.3f", 100* sigmapierr);
//    dummy3 = sprintf(errX33, "#pm %.3f", 100* dERespierr);
//    strcat(txtX31, errX31);
//    strcat(txtX32, errX32);
//    strcat(txtX33, errX33);

    TString txtX31 = "Fitted Mean #mu = ";
    TString txtX32 = "#sigma/#mu [%] = ";
    TString txtX33 = "#Delta E/E [%] = ";

    txtX31+= round1(double(meanpi)*scalepi);
    txtX32+= round3(100*sigmapi);
    txtX33+= round3(100*dERespi);

    txtX31+= " #pm ";
    txtX32+= " #pm ";
    txtX33+= " #pm ";

    txtX31+= round1(double(meanpierr)*scalepi);
    txtX32+= round3(100*sigmapierr);
    txtX33+= round3(100*dERespierr);
    if (verbose){
        std::cout<<"fitted mean: "<<meanpi<<"+/-"<<meanpierr<<", sigma/mean: "<<sigmapi<<"+/-"<<sigmapierr<<", Delta E/E: "<<dERespi<<"=/-"<<dERespierr<<std::endl;
    }

//    ptX3->AddText(txtX31);
//    ptX3->AddText(txtX32);
//    ptX3->AddText(txtX33);

    TString pionstitle = "Pions";
    ptX3->AddText(pionstitle.Data());
    ptX3->AddText(txtX31.Data());
    ptX3->AddText(txtX32.Data());
//    ptX3->AddText(txtX33.Data());


    //------------------------------------------------------------
    //Electron fit
//    cdEdx->cd(4);
    cdEdxelectron->cd();

    Double_t parelectronfit[15];
    for (int i =0; i<15; i++) parelectronfit[i] = 0;
    Double_t parelectronfiterrors[15];
    for (int i =0; i<15; i++) parelectronfiterrors[i] = 0;

    int electronmean = 0;
    int electronhwhm = 0;

    double scaleel = 1.;

    FindParametersSingleGauss(histelectron, parelectronfit, &electronmean, &electronhwhm, scaleel);
    FitHistogramSingleGauss(histelectron, parelectronfit,&parelectronfiterrors[0], scaleel);
    if (verbose){
        cout << "Electron Fit, Electron mean: "<<electronmean<<"  \n \n";
    }

    // Create Lables with Fit results

    TPaveText *ptX4 = new TPaveText(0.7, 0.55, 0.95, 0.71, "NDC");
    ptX4->SetFillColor(10);

//    char valX41[10] = "";
//    char errX41[10] = "";
//    char txtX41[25] = "Fitted Mean #mu = ";
//    char valX42[10] = "";
//    char errX42[10] = "";
//    char txtX42[25] = "#sigma/#mu [%] = ";
//    char valX43[10] = "";
//    char errX43[10] = "";
//    char txtX43[25] = "#Delta E/E [%] = ";
//    int dummy4 = 0;

    Double_t meanel = parelectronfit[1];
    Double_t sigmael = parelectronfit[2]/parelectronfit[1];
    Double_t dEResel = sigmael*2.3548;

    Double_t meanelerr = parelectronfiterrors[1];
    Double_t sigmaelerr = parelectronfit[2]/parelectronfit[1]*TMath::Sqrt(TMath::Power((parelectronfiterrors[2]/parelectronfit[2]),2)+TMath::Power((parelectronfiterrors[1]/parelectronfit[1]),2));
    Double_t dEReselerr = sigmaelerr*2.3548;

//    dummy4 = sprintf(valX41, "%.3f" , double(meanel)*scaleel);
//    dummy4 = sprintf(valX42, "%.3f", 100*sigmael);
//    dummy4 = sprintf(valX43, "%.3f", 100*dEResel);
//    strcat(txtX41, valX41);
//    strcat(txtX42, valX42);
//    strcat(txtX43, valX43);
//
//    dummy4 = sprintf(errX41, "#pm %.3f" , double(meanelerr)*scaleel);
//    dummy4 = sprintf(errX42, "#pm %.3f", 100*sigmaelerr);
//    dummy4 = sprintf(errX43, "#pm %.3f", 100*dEReselerr);
//    strcat(txtX41, errX41);
//    strcat(txtX42, errX42);
//    strcat(txtX43, errX43);

    TString txtX41 = "Fitted Mean #mu = ";
    TString txtX42 = "#sigma/#mu [%] = ";
    TString txtX43 = "#Delta E/E [%] = ";

    txtX41+= round1(double(meanel)*scaleel);
    txtX42+= round3(100*sigmael);
    txtX43+= round3(100*dEResel);

    txtX41+= " #pm ";
    txtX42+= " #pm ";
    txtX43+= " #pm ";



    txtX41+= round1(double(meanelerr)*scalepi);
    txtX42+= round3(100*sigmaelerr);
    txtX43+= round3(100*dEReselerr);

    if (verbose){
        std::cout<<"fitted mean: "<<meanel<<"+/-"<<meanelerr<<", sigma/mean: "<<sigmael<<"+/-"<<sigmaelerr<<", Delta E/E: "<<dEResel<<"=/-"<<dEReselerr<<std::endl;
    }

//    ptX4->AddText(txtX41);
//    ptX4->AddText(txtX42);
//    ptX4->AddText(txtX43);

    TString electronstitle = "Electrons";
    ptX4->AddText(electronstitle.Data());
    ptX4->AddText(txtX41.Data());
    ptX4->AddText(txtX42.Data());
//    ptX4->AddText(txtX43.Data());


    //------------------------------------------------------------
    cdEdx->cd(1);

    // Declare ROOT variables
    TPaveText *ptX1;
    // Others:
    double separationpower1=0.;


    GetSeparationPower(histall, &separationpower1, pionmean, electronmean, pionhwhm, electronhwhm, 1.0);

    if (verbose){
        cout << "\n\n===============================\n\n";
        cout << "Fitting done... \n\n";

        cout << "Separation Power:" << separationpower1 << "\n";
    }

    // Create Lables with Fit results


    ptX1 = new TPaveText(0.75, 0.87, 0.95, 0.98, "NDC");
    ptX1->SetFillColor(10);

    char valX1[10] = "";
    char txtX1[25] = "Separation Power = ";

    int dummy1 = 0;

    dummy1 = sprintf(valX1, "%.1f", separationpower1);
    strcat(txtX1, valX1);
    ptX1->AddText(txtX1);
    ptX1->Draw();


    //------------------------------------------------------------

//    cdEdx->cd(2);
    cdEdxbothpeaks->cd();

    // Declare ROOT variables
    TPaveText *ptX2;
    // Others:
    double separationpower2=0.;


    GetSeparationPower(histPID, &separationpower2, pionmean, electronmean, pionhwhm, electronhwhm, 1.0);

    if (verbose){
        cout << "\n\n===============================\n\n";
        cout << "Fitting done... \n\n";

        cout << "Separation Power:" << separationpower2 << "\n";
    }

    // Create Lables with Fit results


    ptX2 = new TPaveText(0.7, 0.87, 0.95, 0.98, "NDC");
    ptX2->SetFillColor(10);

    char valX2[10] = "";
    char txtX2[25] = "Separation Power = ";

    int dummy2 = 0;

    dummy2 = sprintf(valX2, "%.1f", separationpower2);
    strcat(txtX2, valX2);
    ptX2->AddText(txtX2);
    ptX2->Draw();
    ptX3->Draw();
    ptX4->Draw();

    TString savetemp = "/nfs/mnemosyne/user/pgadow/private/BAThesis/spectra/" + inttostring(runnumber) + ".pdf";
    TString savetemproot = "/nfs/mnemosyne/user/pgadow/private/BAThesis/spectra/" + inttostring(runnumber) + ".root";
    cdEdxbothpeaks->SaveAs(savetemp.Data());
    cdEdxbothpeaks->SaveAs(savetemproot.Data());

    //------------------------------------------------------------

//    savefile->cd();
//    cdEdx->Write();
//    cdEdxbothpeaks->Write();
//    cdEdxelectron->Write();
//    cdEdxpion->Write();
//    savefile->Close();
    if(verbose){
        cout << "Goodbye...\n";
    }

    //runnumber datasetnumber electronmean electronsigma pionmean pionsigma separationpower
    outputfile << inttostring(runnumber) << " " << inttostring(datasetnumber)<< " "
            << inttostring(lowtrunc)<< " " << inttostring(hightrunc)<< " "
            << doubletostring(meanel) << " " << doubletostring(meanelerr) << " "
            << doubletostring(sigmael) << " " << doubletostring(sigmaelerr) << " "
            << doubletostring(meanpi)<< " " << doubletostring(meanpierr)<< " "
            << doubletostring(sigmapi)<< " " << doubletostring(sigmapierr)<< " "
            << doubletostring(separationpower1)<< " " << doubletostring(0)<< " "
            << doubletostring(separationpower2)<< " " << doubletostring(0)<< "\n";

    outputfile.close();
    return;
}