//Macro to test the prediction of Allison and Cobb
//Author: P. Gadow
//status: version 1.0

#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 "TTree.h"
#include "TString.h"
#include "TH1F.h"
#include "TMath.h"
#include "TVector3.h"
#include "TClonesArray.h"
#include "TStyle.h"
#include "TPaveText.h"



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


#include <TpcSample.h>
#include <TpcEvent.h>
#include <TpcCluster.h>
#include <TpcDigi.h>
#include <TpcdEdx.h>
#include <GEMEvent.h>
#include <GFTrack.h>



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

string doubletostring(double input)
{
    string s;
    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,500);
    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,500);
    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 ExecutedEdxfit(TH1F *histpion, TH1F *histelectron, double& resolutionpion, double& resolutionpionerror, double& resolutionelectron, double& resolutionelectronerror){ //TString histofile


    // create canvas
    TCanvas *cdEdx = new TCanvas("cdEdxFit","dEdxFit",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->Divide(2,1);
    cdEdx->SetGrid();

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

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

    cdEdx->cd(1);

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

    FitHistogramSingleGauss(histpion, parpionfit, &parpionfiterrors[0], scalepi);


    // Create Lables with Fit results

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

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

    resolutionpion = sigmapi;

    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;

    resolutionpionerror = sigmapierr;


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


    ptX3->AddText(txtX31.Data());
    ptX3->AddText(txtX32.Data());
    ptX3->AddText(txtX33.Data());
    ptX3->Draw();


    //------------------------------------------------------------
    //Electron fit
    cdEdx->cd(2);

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

    // Create Lables with Fit results

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


    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;

    resolutionelectron = sigmael;

    resolutionelectronerror = sigmaelerr;


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


    ptX4->AddText(txtX41.Data());
    ptX4->AddText(txtX42.Data());
    ptX4->AddText(txtX43.Data());
    ptX4->Draw();

    cdEdx->Write();

    return;
}








//##############################################################################################################
//Fitting Task End #############################################################################################
//##############################################################################################################





bool dedxSort(std::pair<double,double> p1,std::pair<double,double> p2){
  if(fabs(p1.second)<1.E-10) return false;
  if(fabs(p2.second)<1.E-10) return true;
  return p1.first/p1.second < p2.first/p2.second;
}


double truncatedmeanNsamples(unsigned int nsamples, double trunclow, double trunchigh, TpcdEdx* tpcdedx)
{
    //check if track has less than nsamples, then abort
    if (tpcdedx->nEntries()<nsamples){
        return -1;
    }

    std::vector < std::pair<double,double> > dEdxsamples;

    //the samples from TpcdEdxTaskAlice are already sorted for rows, first sample is from first row
    for (unsigned int i = 0;  i<nsamples; i++){
        std::pair<double, double> dedxsample = tpcdedx->getdEdx(i);
        dEdxsamples.push_back(dedxsample);
    }

    sort(dEdxsamples.begin(),dEdxsamples.end(),dedxSort);


    double cutLow = trunclow;
    double cutHigh = trunchigh;  //throw away this percentage of entries (high tail)

    int N = dEdxsamples.size();

    int intCutLow = (int)(cutLow*N);
    int intCutHigh = (int)(cutHigh*N);

    double sum=0.;
    int NinSum=0;
    for(int i=intCutLow;i<(N-intCutHigh);++i){
          if(fabs(dEdxsamples[i].second)>1.E-10){
            sum+=dEdxsamples[i].first/dEdxsamples[i].second;
            ++NinSum;
          }
    }
    if(NinSum>0)return sum/NinSum;
    return -1.;
}




void ExecuteACPrediction(bool verbose, int runnumber, int datasetnumber, int numevents, unsigned int nsamples, double cutlow=0., double cuthigh=.4){


    int lowtrunc = int ((cutlow+0.005) *100);
    int hightrunc = int ((cuthigh+0.005) *100);

    TString indirname = "/nfs/mds/data/tpc/alice/ps_test_beam/RecoOut/";
    TString infilename = "merge_" + inttostring(runnumber) + "_"+inttostring(datasetnumber)+"_corrected_hough_smoothed.reco.root";
    TString inpathname = indirname + infilename;

    //define output filenames
    TString dirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdx/ACP/";
    TString filename = "merge_" + inttostring(runnumber) + "_"+ inttostring(datasetnumber) + "_dEdxACP";
    filename+= "_truncmean_"+inttostring(lowtrunc)+"_"+inttostring(hightrunc);
    TString savefilename = dirname + filename + ".root";
    TFile *savefile = new TFile(savefilename.Data(),"UPDATE");

    ofstream outputfile;
    outputfile.open (dirname + filename + "_ResolutionACP.dat", std::ofstream::app);

    // open file(s)
    TFile *reco = new TFile(inpathname);

    // create canvas
    TCanvas *cdEdx = new TCanvas("cdEdx","dEdx",1000,700);
    TCanvas *cTracklength = new TCanvas("cTracklength","underlying distribution for cut",1000,700);
    cdEdx->SetFillColor(00);
    cdEdx->Divide(2,2);
    cdEdx->SetGrid();


    TString histname = "dEdx";
    histname+=" truncated mean ("+doubletostring(cutlow)+", "+doubletostring(cuthigh)+") -";
    histname+= " single track events";

    TString histnamePID = histname+ ": e and Pion cut";
    TString histnameE1 = histname + ": e-cher cut(>400)";
    TString histnameE2 = histname + ": e- double cut(>400)";
    TString histnamePi = histname + ": Pionen";

    // create histogram
    TH1F *histall = new TH1F ("histall", histname, 600,0,600);
    TH1F *histPID = new TH1F ("histPID", histnamePID, 600,0,600);
    TH1F *histelEsinglecut = new TH1F ("histEsinglecut", histnameE1, 600,0,600);
    TH1F *histEledoublecut = new TH1F ("histEdoublecut", histnameE2, 600,0,600);
    TH1F *histpion = new TH1F ("histpion", histnamePi, 400,0,400);
    TH1F *histtracklength = new TH1F ("histtracklength", "TrackLength", 160,0,80);
    // create Tree
    TTree *tpcTree = (TTree*)reco->Get("cbmsim");


    //save time and load only necessary branches
    tpcTree->SetBranchStatus("*",0);
    tpcTree->SetBranchStatus("TrackPostFit.*",1);
    tpcTree->SetBranchStatus("GEMEvent.*",1);
    tpcTree->SetBranchStatus("dEdx.*",1);

    // create arrays
    TClonesArray *dedx = new TClonesArray("TpcdEdx");
    TClonesArray *tracks = new TClonesArray("GFTrack");
    TClonesArray *gemevent = new TClonesArray("GEMEvent");


    tpcTree->SetBranchAddress("dEdx",&dedx);
    tpcTree->SetBranchAddress("TrackPostFit",&tracks);
    tpcTree->SetBranchAddress("GEMEvent",&gemevent);

    savefile->cd();

    // loop over dedx

    //events
    if (numevents == 0){
        numevents = tpcTree->GetEntries();
    }

    std::cout<<std::endl;
    for( Int_t ev = 0; ev <numevents; ev++){

        if(ev%1000==0){
            std::cout<<"#################### Event done: "<<ev<<"##############"<<std::endl;
        }

        dedx->Clear();
        tracks->Clear();
        gemevent->Clear();

        tpcTree->GetEvent(ev);

        GEMEvent* gemEvent = (GEMEvent*)gemevent->At(0); //there is only one element in gemevent
        Int_t NumTpcdEdx = dedx->GetEntries();
        if (NumTpcdEdx == 0) continue;


        int numberoftracks = tracks->GetEntries();
        if(verbose){
            std::cout << "number of tracks per event:" << numberoftracks << std::endl;
        }
        //select only single track events
        if (numberoftracks!=1){
                continue;
            }

        for (Int_t idedex = 0; idedex < NumTpcdEdx; ++idedex){
            TpcdEdx* tpcdedx = (TpcdEdx*) dedx->At(idedex);


            //tracklength
            //detector length in x = 46.5 cm
            double tracklength = 0.;
            for(unsigned int idx=0;idx<tpcdedx->nEntries();++idx){
                tracklength += tpcdedx->getdx(idx);
            }
            histtracklength->Fill(tracklength);

            // select only tracks that go all over the detetor
            if (tracklength < 35){ ///
                continue;
            }
            TVector3 position = ((GFTrack*) tracks->At(0))->getCardinalRep()->getPos();

            //select only tracks with z-position above 2.713 cm

            if (position.Z() < 2.7){
                continue;
            }

            if (position.Z()>10.6){
	      continue;
	    }

            double mean = 0;

            //calculate the reduced truncmean and put it in a map with the number of samples as the first index
            mean = truncatedmeanNsamples(nsamples, cutlow, cuthigh, tpcdedx);

            if (mean<0){
                continue;
            }

            if(verbose){
                std::cout<< "num of dedx entries:" << NumTpcdEdx <<std::endl;
                std::cout<< "num of dedx entries in single dedx:" << tpcdedx->nEntries() <<std::endl;
                for(unsigned int i = 0; i < tpcdedx->nEntries(); i++){
                    std::cout<<"dedx" << idedex<<"       de: " << tpcdedx->getdx(i)<< ", dx: "<<tpcdedx->getdE(i)<<std::endl;
                }
                std::cout<<mean<<std::endl;
            }
            histall->Fill(mean);
            // is it an electron?
            if (gemEvent->GetCherenkovValue()>400){
                histelEsinglecut->Fill(mean);
                histPID->Fill(mean);
            }
            if (gemEvent->GetCherenkovValue()>400 && gemEvent->GetPbGlassValue()>400){
                histEledoublecut->Fill(mean);
            }

            // is it a pion?
            if (gemEvent->GetCherenkovValue() < 150){
                histpion->Fill(mean);
                histPID->Fill(mean);
            }

        }//end ofdedxloop
    }//end of eventloop

    double respion = 0;
    double respionerr = 0;
    double reselectron = 0;
    double reselectronerr = 0;

    ExecutedEdxfit(histpion, histelEsinglecut, respion,respionerr,reselectron,reselectronerr);


    outputfile << inttostring(runnumber) << " " << inttostring(datasetnumber)<< " "
            << inttostring(lowtrunc)<< " " << inttostring(hightrunc)<< " "
            << inttostring(nsamples) << " "
            << doubletostring(respion) << " " << doubletostring(respionerr) << " "
            << doubletostring(reselectron)<< " " << doubletostring(reselectronerr) << "\n";

    outputfile.close();

    //  Draw the histogram
    cdEdx->cd(1);
    histall->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histall->GetYaxis()->SetTitle("# of entries");
    histall->Write();
    cdEdx->cd(2);
    histpion->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histpion->GetYaxis()->SetTitle("# of entries");
    histpion->Write();
    cdEdx->cd(3);
    histPID->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histPID->GetYaxis()->SetTitle("# of entries");
    histPID->Write();

    cdEdx->cd(4);
    histelEsinglecut->GetXaxis()->SetTitle("dE/dx (a.u.)");
    histelEsinglecut->GetYaxis()->SetTitle("# of entries");
    histelEsinglecut->Write();

    histEledoublecut->GetYaxis()->SetTitle("# of entries");
    histEledoublecut->Write();


    cdEdx->Write();

    cTracklength->cd();
    histtracklength->Write();
    std::cout<<std::endl;


    savefile->Close();
  }