#include "TGraph.h"
#include "TGraphErrors.h"
#include "TString.h"
#include "TMultiGraph.h"
#include "TLegend.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "TAxis.h"
#include "TROOT.h"
#include "TFile.h"
#include "TMath.h"
#include <map>
#include <iostream>
#include <sstream>
#include <fstream>


//Definition of ACP for Resolution

Double_t ACPred(Double_t *N, Double_t  *par)
{   Double_t n = N[0];
    Double_t prefactor = 96./2.35*TMath::Power((0.75 * 0.938),-0.32); //prediction for pure argon
    prefactor *= TMath::Power(((11.2)/(35.2)),-0.32) / TMath::Power(((18)/(26.4)),-0.32); //scaling for NeCo2: (mean # electrons gas)/(mean ionization potential gas)^(-0.32)/(mean # electrons Ar)/(mean ionization potential Ar)^(-0.32)
    prefactor *= TMath::Power(par[0],2);
    Double_t resolution = prefactor * TMath::Power(n,-0.43);
    return resolution;
    
}



struct FitResults
{
    unsigned int run;
    unsigned int data;
    unsigned int lowtrunc;
    unsigned int hightrunc;
    unsigned int nsamples;
    double sigmapi;
    double sigmapierr;
    double sigmael;
    double sigmaelerr;
};

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


//temporary solution: a map - better with root trees
std::map<unsigned int, FitResults> fitresultsmap;


void LoadData(TString datafile, unsigned int datasetnumber, unsigned int setlowtrunc, unsigned int sethightrunc){
    ifstream instr;
    instr.open(datafile.Data());

    if (instr.is_open())
        std::cout<<"opened file: "<<datafile.Data()<<std::endl;
    else{
        std::cerr << "Error: file not found!!! Aborting";
        return;
    }

    while(1) {
        unsigned int run = 0;
        unsigned int data=0;
        unsigned int lowtrunc=0;
        unsigned int hightrunc=0;
        unsigned int nsamples=0;
        double sigmapi=0;
        double sigmapierr=0;
        double sigmael=0;
        double sigmaelerr=0;



        // input from file: 703 2 5 25 20  0.191379 0.005231 0.184228 0.00544.
        instr>>run>>data>>lowtrunc>>hightrunc>>nsamples>>sigmapi>>sigmapierr>>sigmael>>sigmaelerr;
        if(instr.eof())
            break;
        cout << "reading data: " << run <<" "<< data << " "<< lowtrunc <<" "<< hightrunc<<std::endl;


        // select correct data
        if (data!=datasetnumber){
            continue;
        }

        if (lowtrunc!=setlowtrunc){
            continue;
        }

        if (hightrunc!=sethightrunc){
            continue;
        }


        if( sigmael < 0||   sigmapi< 0||   sigmaelerr < 0||  sigmapierr< 0)  {
            std::cerr << "dEdxACPMakePictures:: Error in Fit, negative Results in run " << run << " with data " << data << std::endl;
            continue;
        }
        if(fitresultsmap.find(run) != fitresultsmap.end() && fitresultsmap.find(data) != fitresultsmap.end()) {
            std::cerr << "Error: double run number in data!!! Aborting";
            return;
        }

        FitResults results;
        results.run = run;
        results.lowtrunc = lowtrunc;
        results.hightrunc = hightrunc;
        results.nsamples = nsamples;
        results.sigmael = sigmael;
        results.sigmapi = sigmapi;
        results.sigmaelerr = sigmaelerr;
        results.sigmapierr = sigmapierr;

        cout << "data read! run: "<< run <<std::endl;

        
        unsigned int fitresultsentry = run*100;
        fitresultsentry +=nsamples;
        
        fitresultsmap[fitresultsentry]= results;
        cout << "data written in map" <<std::endl;

    }

    std::cout<<"Data successfully read"<<std::endl;

    return;

}

//##############################################################################################################

void ExecuteMakeACPPictures(int runnumber, int datasetnumber, unsigned int trunclow, unsigned int trunchigh, bool pretty = false) {

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

    }

    TString dirname = "/nfs/mds/data/tpc/alice/ps_test_beam/dEdx/ACP/";
    
    TString infilename = dirname + "merge_" + inttostring(runnumber) + "_" + inttostring(datasetnumber) + "_dEdxACP_truncmean_" + inttostring(trunclow) + "_" + inttostring(trunchigh) + "_ResolutionACP.dat";
    
    LoadData(infilename, datasetnumber, trunclow, trunchigh);

    

    TString filename = "dEdxACPResults_";
    filename+= inttostring(runnumber) + "_"+inttostring(datasetnumber)+"_"+inttostring(trunclow)+"_"+inttostring(trunchigh);

    TString savefilename = dirname + filename + ".root";
    TFile *savefile = new TFile(savefilename.Data(),"RECREATE");
    savefile->cd();

    // Plot Graphs
    
    //################################################################################################################################################
    //Resolution for different truncations ###########################################################################################################
    //################################################################################################################################################

    TCanvas *cstdsigma = new TCanvas("cstdsigma","dEdx resolution for different truncations",1000,700);
    cstdsigma->SetFillColor(00);
    cstdsigma->SetGrid();
    cstdsigma->cd();
    TMultiGraph *mgstdsigma = new TMultiGraph();


    const int nstdsigma = 10;
    double samplesstdsigma[nstdsigma] = {10., 15., 20., 25., 30., 35., 40., 45., 50., 55.};

    double sigmael[nstdsigma];
    double sigmapi[nstdsigma];
    double dsigmael[nstdsigma];
    double dsigmapi[nstdsigma];
    unsigned int mapentry = runnumber*100;
    sigmael[0] = (fitresultsmap[mapentry+10]).sigmael*100; sigmael[1] = (fitresultsmap[mapentry+15]).sigmael*100; sigmael[2] = (fitresultsmap[mapentry+20]).sigmael*100; sigmael[3] = (fitresultsmap[mapentry+25]).sigmael*100; sigmael[4] = (fitresultsmap[mapentry+30]).sigmael*100; sigmael[5] = (fitresultsmap[mapentry+35]).sigmael*100; sigmael[6] = (fitresultsmap[mapentry+40]).sigmael*100; sigmael[7] = (fitresultsmap[mapentry+45]).sigmael*100; sigmael[8] = (fitresultsmap[mapentry+50]).sigmael*100; sigmael[9] = (fitresultsmap[mapentry+55]).sigmael*100;
    sigmapi[0] = (fitresultsmap[mapentry+10]).sigmapi*100; sigmapi[1] = (fitresultsmap[mapentry+15]).sigmapi*100; sigmapi[2] = (fitresultsmap[mapentry+20]).sigmapi*100; sigmapi[3] = (fitresultsmap[mapentry+25]).sigmapi*100; sigmapi[4] = (fitresultsmap[mapentry+30]).sigmapi*100; sigmapi[5] = (fitresultsmap[mapentry+35]).sigmapi*100; sigmapi[6] = (fitresultsmap[mapentry+40]).sigmapi*100; sigmapi[7] = (fitresultsmap[mapentry+45]).sigmapi*100; sigmapi[8] = (fitresultsmap[mapentry+50]).sigmapi*100; sigmapi[9] = (fitresultsmap[mapentry+55]).sigmapi*100;
    dsigmael[0] = (fitresultsmap[mapentry+10]).sigmaelerr*100; dsigmael[1] = (fitresultsmap[mapentry+15]).sigmaelerr*100; dsigmael[2] = (fitresultsmap[mapentry+20]).sigmaelerr*100; dsigmael[3] = (fitresultsmap[mapentry+25]).sigmaelerr*100; dsigmael[4] = (fitresultsmap[mapentry+30]).sigmaelerr*100; dsigmael[5] = (fitresultsmap[mapentry+35]).sigmaelerr*100; dsigmael[6] = (fitresultsmap[mapentry+40]).sigmaelerr*100; dsigmael[7] = (fitresultsmap[mapentry+45]).sigmaelerr*100; dsigmael[8] = (fitresultsmap[mapentry+50]).sigmaelerr*100; dsigmael[9] = (fitresultsmap[mapentry+55]).sigmaelerr*100;
    dsigmapi[0] = (fitresultsmap[mapentry+10]).sigmapierr*100; dsigmapi[1] = (fitresultsmap[mapentry+15]).sigmapierr*100; dsigmapi[2] = (fitresultsmap[mapentry+20]).sigmapierr*100; dsigmapi[3] = (fitresultsmap[mapentry+25]).sigmapierr*100; dsigmapi[4] = (fitresultsmap[mapentry+30]).sigmapierr*100; dsigmapi[5] = (fitresultsmap[mapentry+35]).sigmapierr*100; dsigmapi[6] = (fitresultsmap[mapentry+40]).sigmapierr*100; dsigmapi[7] = (fitresultsmap[mapentry+45]).sigmapierr*100; dsigmapi[8] = (fitresultsmap[mapentry+50]).sigmapierr*100; dsigmapi[9] = (fitresultsmap[mapentry+55]).sigmapierr*100;

    

    //Plot Allison Cobb prediction
    TF1* grACPredictionEl = new TF1("ACPred",ACPred,0,60,1);
    grACPredictionEl->SetParameter(0,1); //set beta to 1 for electrons and pions (pions: beta = 0.990)
    grACPredictionEl->SetLineColor(kRed);
    grACPredictionEl->SetLineStyle(2);


    //Plot measured resolutions
    
    TGraphErrors* grsigmael = new TGraphErrors(nstdsigma,samplesstdsigma,sigmael,0,dsigmael);
    grsigmael->SetLineColor(1);
    grsigmael->SetLineWidth(1);
    grsigmael->SetMarkerColor(1);
    grsigmael->SetMarkerStyle(20);
    grsigmael->SetTitle("electrons");
    mgstdsigma->Add(grsigmael);

    TGraphErrors* grsigmapi = new TGraphErrors(nstdsigma,samplesstdsigma,sigmapi,0,dsigmapi);
    grsigmapi->SetLineColor(15);
    grsigmapi->SetLineWidth(1);
    grsigmapi->SetMarkerColor(15);
    grsigmapi->SetMarkerStyle(20);
    grsigmapi->SetTitle("pions");
    mgstdsigma->Add(grsigmapi);

    mgstdsigma->Draw("ap");
    mgstdsigma->GetXaxis()->SetTitle("# of dE/dx samples");
    mgstdsigma->GetYaxis()->SetTitle("#sigma(dE/dx)/<dE/dx> [%]");
    mgstdsigma->GetXaxis()->SetLimits(0,60);
    mgstdsigma->GetYaxis()->SetRangeUser(9,30);

    TLegend *legstdsigma=new TLegend(0.6,0.78,0.9,.9);
    legstdsigma->SetBorderSize(1);
    legstdsigma->SetFillColor(10);
    legstdsigma->AddEntry(grsigmael,"electrons","p");
    legstdsigma->AddEntry(grsigmapi,"pions","p");
    legstdsigma->AddEntry(grACPredictionEl,"PAI prediction","l");
    legstdsigma->Draw("same");

    cstdsigma->Update();
    cstdsigma->Modified();
    grACPredictionEl->Draw("Same");
    cstdsigma->Write();
    cstdsigma->SaveAs(dirname+filename+"ACP.pdf");

    //-----------------------------------------------------------------------------------
    savefile->Close();

}