/**
 * \file CbmRichGeoTestOpt.cxx
 * \author Semen Lebedev <s.lebedev@gsi.de>
 * \date 2019
 */
#include "CbmRichGeoTestOpt.h"

#include "CbmDrawHist.h"
#include "utils/CbmUtils.h"

#include "TFile.h"
#include <iostream>
#include <utility>
#include "TLine.h"
#include "TStyle.h"

using namespace std;

CbmRichGeoTestOpt::CbmRichGeoTestOpt():
                        fHM(nullptr)
{

}

CbmRichGeoTestOpt::~CbmRichGeoTestOpt()
{

}

void CbmRichGeoTestOpt::SetFilePathes(vector<vector<string> > geoTestPathes, vector<vector<string> > urqmdTestPathes, vector<vector<string> > recoQaPathes)
{
    fGeoTestPathes = geoTestPathes;
    fUrqmdTestPathes = urqmdTestPathes;
    fRecoQaPathes = recoQaPathes;
}

int CbmRichGeoTestOpt::GetNofFiles()
{
    return fGeoTestPathes[0].size();
}

string CbmRichGeoTestOpt::GetFilePath(int simInd, CbmRichGeoTestOptFileEnum fileEnum, int iFile)
{
    if (fileEnum == kGeoTestFile) return fGeoTestPathes[simInd][iFile];
    if (fileEnum == kUrqmdTestFile) return fUrqmdTestPathes[simInd][iFile];
    if (fileEnum == kRecoQaFile) return fRecoQaPathes[simInd][iFile];

    return "";
}

pair<double, double> CbmRichGeoTestOpt::H1MeanRms(int simInd, CbmRichGeoTestOptFileEnum fileEnum, int iFile, const string& histName)
{
    string path = GetFilePath(simInd, fileEnum, iFile);
    if (path == "") return make_pair(0., 0.);
    TFile* file = new TFile(path.c_str(), "READ");
    if (file == nullptr) return make_pair(0., 0.);
    TH1D* hist = (TH1D*)file->Get(histName.c_str());
    if (hist == nullptr) return make_pair(0., 0.);
    double mean = hist->GetMean();
    double rms = hist->GetRMS();
    file->Close();
    delete file;
    return make_pair(mean, rms);
}

pair<double, double> CbmRichGeoTestOpt::H2ProjYMeanRms(int simInd, CbmRichGeoTestOptFileEnum fileEnum, int iFile, const string& histName)
{
    string path = GetFilePath(simInd, fileEnum, iFile);
    if (path == "") return make_pair(0., 0.);
    TFile* file = new TFile(path.c_str(), "READ");
    if (file == nullptr) return make_pair(0., 0.);
    TH2D* hist = (TH2D*)file->Get(histName.c_str());
    if (hist == nullptr) return make_pair(0., 0.);
    TH1D* py = hist->ProjectionY( (histName + to_string(iFile) + "_py" ).c_str() );
    double mean = py->GetMean();
    double rms = py->GetRMS();
    file->Close();
    delete file;
    return make_pair(mean, rms);
}

double CbmRichGeoTestOpt::HEntries(int simInd, CbmRichGeoTestOptFileEnum fileEnum, int iFile, const string& histName)
{
    string path = GetFilePath(simInd, fileEnum, iFile);
    if (path == "") return 0.;
    TFile* file = new TFile(path.c_str(), "READ");
    if (file == nullptr) return 0.;
    TH1* hist = (TH1*)file->Get(histName.c_str());
    if (hist == nullptr) return 0.;
    double entries = hist->GetEntries();
    file->Close();
    delete file;
    return entries;
}


void CbmRichGeoTestOpt::DrawMeanRms(int simInd, CbmRichGeoTestOptFileEnum fileEnum, const string& histName, CbmRichGeoTestOptHistEnum histEnum,
        const string& titleY, double minY, double maxY, int nofParts, int nofFilesPart)
{
    string canvasName = "richgeoopt_sim" + to_string(simInd) + "_" + histName + to_string(nofParts);
    TCanvas* c = fHM->CreateCanvas(canvasName, canvasName, 1800, 600);
    vector<TH1*> hists;
    vector<string> legend;
    for (int iP = 0; iP < nofParts; iP++) {
        string histNameInd = "sim" + to_string(simInd) + histName + "_" + to_string(nofParts) + "_" +to_string(iP);
        fHM->Create1<TH1D>(histNameInd, histNameInd + ";Geometry;" + titleY, nofFilesPart, .5, nofFilesPart + .5);
        for (int iF = 0; iF < nofFilesPart; iF++) {
            double value = 0.;
            int ind = iP * nofFilesPart + iF;
            if (histEnum == kH1MeanHist) value = H1MeanRms(simInd, fileEnum, ind, histName).first;
            else if (histEnum == kH1RmsHist) value = H1MeanRms(simInd, fileEnum, ind, histName).second;
            else if (histEnum == kH2MeanHist) value = H2ProjYMeanRms(simInd, fileEnum, ind, histName).first;
            else if (histEnum == kH2RmsHist) value = H2ProjYMeanRms(simInd, fileEnum, ind, histName).second;
            else if (histEnum == kHEntriesHist) value = HEntries(simInd, fileEnum, ind, histName) / 1000.; // TODO: get nof events

            fHM->H1(histNameInd)->SetBinContent(iF+1, value);
        }
        hists.push_back(fHM->H1(histNameInd));
        legend.push_back(to_string(iP));
    }

    DrawManyH1(hists,legend, minY, maxY);
}


void CbmRichGeoTestOpt::DrawMeanEff(int simInd, CbmRichGeoTestOptFileEnum fileEnum, const string& histName1, const string& histName2,
        const string& titleY, double minY, double maxY, int nofParts, int nofFilesPart)
{
    string canvasName = "richgeoopt_sim" + to_string(simInd) + "_" + histName1 + "_" + histName2 + to_string(nofParts);
    TCanvas* c = fHM->CreateCanvas(canvasName, canvasName, 1800, 600);
    vector<TH1*> hists;
    vector<string> legend;
    for (int iP = 0; iP < nofParts; iP++) {
        string histEffName = "sim" + to_string(simInd) + histName1 + "_" + histName2 + "_" + to_string(nofParts) + "_" +to_string(iP);
        //int nofFilesPart = (GetNofFiles() / nofParts) + 1;
        fHM->Create1<TH1D>(histEffName, histEffName + ";Geometry;"+titleY, nofFilesPart, .5, nofFilesPart + .5);
        for (int iF = 0; iF < nofFilesPart; iF++)  {
            int ind = iP * nofFilesPart + iF;
            double entries1 = HEntries(simInd, fileEnum, ind, histName1);
            double entries2 = HEntries(simInd, fileEnum, ind, histName2);
            double value = (entries2 != 0)?100. * entries1 / entries2 : 0.;
            fHM->H1(histEffName)->SetBinContent(iF+1, value);
        }
        hists.push_back(fHM->H1(histEffName));
        legend.push_back(to_string(iP));
    }
    DrawManyH1(hists, legend, minY, maxY);
}

void CbmRichGeoTestOpt::DrawMeanRms2D(int simInd, CbmRichGeoTestOptFileEnum fileEnum, const string& histName, CbmRichGeoTestOptHistEnum histEnum,
        const string& titleZ, double minZ, double maxZ)
{
    string canvasName = "richgeoopt_sim" + to_string(simInd) + "_" + histName + "_2d";
    TCanvas* c = fHM->CreateCanvas(canvasName, canvasName, 1800, 1200);
    c->Divide(3,2);
    vector<string> camTilt = {"-15 deg", "-10 deg", "-5 deg", "0 deg", "5 deg", "10 deg"};
    int nofParts = 6;
    int nofFilesPart = (GetNofFiles() / 30.) + 1;
    int nofX = 7;
    int nofY = 7;
    double shift = 0.5 * 300. / 7.;
    for (int iP = 0; iP < nofParts; iP++) {
        string histNameInd = "sim" + to_string(simInd) + histName + "_" + to_string(iP)+ "_2d";
        fHM->Create2<TH2D>(histNameInd, histNameInd + ";Shift CamY [mm];Shift CamZ [mm];" + titleZ, nofX, -150. - shift, 150 + shift, nofY, -150. - shift, 150 + shift);
        for (int iX = 0; iX < nofX; iX++) {
            for (int iY = 0; iY < nofY; iY++) {
                double value = 0.;
                int ind = iP * nofFilesPart + iX * nofX + iY;
                //cout << iP << " " << iX << " " << iY << " " << ind << endl;
                if (histEnum == kH1MeanHist) value = H1MeanRms(simInd, fileEnum, ind, histName).first;
                else if (histEnum == kH1RmsHist) value = H1MeanRms(simInd, fileEnum, ind, histName).second;
                else if (histEnum == kH2MeanHist) value = H2ProjYMeanRms(simInd, fileEnum, ind, histName).first;
                else if (histEnum == kH2RmsHist) value = H2ProjYMeanRms(simInd, fileEnum, ind, histName).second;
                else if (histEnum == kHEntriesHist) value = HEntries(simInd, fileEnum, ind, histName) / 1000.; // TODO: get nof events

                fHM->H2(histNameInd)->SetBinContent(iX + 1, iY + 1, value);
            }
        }
        c->cd(iP+1);
        DrawH2(fHM->H2(histNameInd), kLinear, kLinear, kLinear, "COLZ TEXT");
        fHM->H2(histNameInd)->GetZaxis()->SetRangeUser(minZ, maxZ);
        fHM->H2(histNameInd)->SetMarkerSize(1.2);
        DrawTextOnPad(camTilt[iP], 0.1, 0.9, 0.9, 0.99);
    }

}


void CbmRichGeoTestOpt::DrawMeanEff2D(int simInd, CbmRichGeoTestOptFileEnum fileEnum, const string& histName1, const string& histName2,
        const string& titleZ, double minZ, double maxZ)
{
    string canvasName = "richgeoopt_sim" + to_string(simInd) + "_" + histName1 + "_" + histName2 + "_2d";
    TCanvas* c = fHM->CreateCanvas(canvasName, canvasName, 1800, 1200);
    c->Divide(3,2);
    vector<string> camTilt = {"-15 deg", "-10 deg", "-5 deg", "0 deg", "5 deg", "10 deg"};
    int nofParts = 6;
    int nofFilesPart = (GetNofFiles() / 30.) + 1;
    int nofX = 7;
    int nofY = 7;
    double shift = 0.5 * 300. / 7.;
    for (int iP = 0; iP < nofParts; iP++) {
        string histNameInd = "sim" + to_string(simInd) + histName1 + "_" + histName2 + "_" + to_string(iP)+ "_2d";
        fHM->Create2<TH2D>(histNameInd, histNameInd + ";Shift CamY [mm];Shift CamZ [mm];" + titleZ, nofX, -150. - shift, 150 + shift, nofY, -150. - shift, 150 + shift);
        for (int iX = 0; iX < nofX; iX++) {
            for (int iY = 0; iY < nofY; iY++) {
                int ind = iP * nofFilesPart + iX * nofX + iY;
                //cout << iP << " " << iX << " " << iY << " " << ind << endl;
                double entries1 = HEntries(simInd, fileEnum, ind, histName1);
                double entries2 = HEntries(simInd, fileEnum, ind, histName2);
                double value = (entries2 != 0)?100. * entries1 / entries2 : 0.;

                fHM->H2(histNameInd)->SetBinContent(iX + 1, iY + 1, value);
            }
        }
        c->cd(iP+1);
        DrawH2(fHM->H2(histNameInd), kLinear, kLinear, kLinear, "COLZ TEXT");
        fHM->H2(histNameInd)->SetMarkerSize(1.2);
        fHM->H2(histNameInd)->GetZaxis()->SetRangeUser(minZ, maxZ);
        DrawTextOnPad(camTilt[iP], 0.1, 0.9, 0.9, 0.99);
    }

}

void CbmRichGeoTestOpt::DrawManyH1(const vector<TH1*>& hist, const vector<string>& legend, double minY, double maxY)
{
    vector<string> legCamRadius = {"1550 mm", "1600 mm","1650 mm", "1700 mm", "1750 mm"};
    vector<string> camTilt = {"-15 deg", "-10 deg", "-5 deg", "0 deg", "5 deg", "10 deg"};
    bool drawLegend = (hist.size()==1)?false:true;
    DrawH1(hist, (legend.size() == 5)?legCamRadius:(legend.size()==6)?camTilt:legend, kLinear, kLinear, drawLegend, 0.93, 0.8, 0.99, 0.99);
    gPad->SetLeftMargin(0.06);
    gPad->SetRightMargin(0.01);
    hist[0]->GetYaxis()->SetTitleOffset(0.5);
    hist[0]->GetYaxis()->SetRangeUser(minY, maxY);
    DrawLines(true, true, false, minY, maxY);
}

void CbmRichGeoTestOpt::DrawLines(bool drawCamRadius, bool drawCamTilt, bool drawCamY, double minY, double maxY)
{
    double nMinY = minY + 0.8 * (maxY - minY);
    double nMaxY = maxY;
    // cam Y
    if (drawCamY) {
        for (int i = 1; i <= 210; i++) {
            TLine* line = new TLine(i * 7 + .5, nMinY, i * 7 + .5, nMaxY);
            line->SetLineColor(kGreen+2);
            line->SetLineWidth(2);
            line->Draw();
        }
    }

    // cam Tilt
    if (drawCamTilt) {
        for (int i = 1; i <= 30; i++) {
            TLine* line = new TLine(i * 49 + .5, nMinY, i * 49 + .5, nMaxY);
            line->SetLineColor(kBlue+2);
            line->SetLineWidth(2);
            line->Draw();
        }
    }

    // radius
    if (drawCamRadius) {
        for (int i = 1; i <= 5; i++) {
            TLine* line = new TLine(i * 294 + .5, nMinY, i * 294 + .5, nMaxY);
            line->SetLineColor(kRed + 1);
            line->SetLineWidth(2);
            line->Draw();
        }
    }

}

void CbmRichGeoTestOpt::Draw()
{
    fHM = new CbmHistManager();

    gStyle->SetPaintTextFormat(".2f");

    SetDefaultDrawStyle();
    int nofFilesPartAll = (GetNofFiles() / 1.) + 1;
    int nofFilesPartCamRadius = (GetNofFiles() / 5.) + 1;
    int nofFilesPartCamTilt = (GetNofFiles() / 30.) + 1;

    for (int iS = 0; iS < 2; iS++) {



    // DrawMeanRms(kGeoTestFile, "fhNofHits_hits", kH1MeanHist, "#hits in ring", 23, 32, 1, nofFilesPartAll);
    // DrawMeanRms(kGeoTestFile, "fhNofHits_points", kH1MeanHist, "#MC points in ring", 160, 200, 1);
    //  DrawMeanRms(kGeoTestFile, "fhBoverAVsMom_hits", kH2MeanHist, "B/A (hit fit)", 0.85, 0.95, 1);
    //  DrawMeanRms(kGeoTestFile, "fhBoverAVsMom_points", kH2MeanHist, "B/A (point fit)", 0.85, 0.95, 1, nofFilesPartAll);
    //
    //    DrawMeanRms(kGeoTestFile, "fhRadiusVsNofHits", kH2RmsHist, "RMS.Radius [cm]", 0.3, 0.75, 1);
    //    DrawMeanRms(kGeoTestFile, "fhAaxisVsNofHits", kH2RmsHist, "RMS.Aaxis [cm]", 0.3, 0.75, 1);
    //    DrawMeanRms(kGeoTestFile, "fhBaxisVsNofHits", kH2RmsHist, "RMS.Baxis [cm]", 0.3, 0.75, 1);
    //    DrawMeanRms(kGeoTestFile, "fhDRVsNofHits", kH2RmsHist, "RMS.dR [cm]", 0.3, 0.45, 1);
    //
    // DrawMeanEff(kGeoTestFile, "fhAccMomEl", "fhMcMomEl", "Geometrical acceptance [%]", 75, 95, 1, nofFilesPartAll);
    //    DrawMeanEff(kGeoTestFile, "fhNofHitsCircleFit", "fhNofHitsAll", "Circle fit efficiency [%]", 93, 99);
    //    DrawMeanEff(kGeoTestFile, "fhNofHitsEllipseFit", "fhNofHitsAll", "Ellipse fit efficiency [%]", 84, 99);
    //
    //
    //    DrawMeanRms(kUrqmdTestFile, "fh_nof_hits_per_event", kH1MeanHist, "Nof hits per event", 500, 850, 1);
    //    DrawMeanRms(kUrqmdTestFile, "fh_nof_rings_1hit", kH1MeanHist, "Nof rings (>=1 hit) per event", 22, 30, 1);
    //    DrawMeanRms(kUrqmdTestFile, "fh_nof_rings_7hits", kH1MeanHist, "Nof rings (>=7hits) per event", 16, 28, 1);
    //    DrawMeanRms(kUrqmdTestFile, "fh_gamma_nontarget_mom", kHEntriesHist, "e^{#pm}_{not target} from #gamma per event", 15, 25, 1);
    //
    //    DrawMeanEff(kRecoQaFile, "hte_Rich_Rich_Electron_Rec_p", "hte_Rich_Rich_Electron_Acc_p", "Ring reco efficiency [%]", 70, 99);
    //    DrawMeanEff(kRecoQaFile, "hte_Rich_Rich_ElectronReference_Rec_p", "hte_Rich_Rich_ElectronReference_Acc_p", "Ring reco efficiency (ref) [%]", 70, 99);
    //    DrawMeanEff(kRecoQaFile, "hte_StsRich_StsRich_Electron_Rec_p", "hte_StsRich_StsRich_Electron_Acc_p", "STS-RICH matching efficiency [%]", 70, 99);


        DrawMeanRms(iS, kGeoTestFile, "fhNofHits_hits", kH1MeanHist, "#hits in ring", 23, 32, 5, nofFilesPartCamRadius);
        DrawMeanRms(iS, kGeoTestFile, "fhBoverAVsMom_points", kH2MeanHist, "B/A (hit fit)", 0.87, 0.96, 5, nofFilesPartCamRadius);
        DrawMeanEff(iS, kGeoTestFile, "fhAccMomEl", "fhMcMomEl", "Geometrical acceptance [%]", 75, 95, 5, nofFilesPartCamRadius);
        DrawMeanRms(iS, kUrqmdTestFile, "fh_nof_hits_per_event", kH1MeanHist, "Nof hits per event", 500, 850, 5, nofFilesPartCamRadius);
        DrawMeanRms(iS, kUrqmdTestFile, "fh_gamma_nontarget_mom", kHEntriesHist, "e^{#pm}_{not target} from #gamma per event", 18, 26, 5, nofFilesPartCamRadius);
        DrawMeanEff(iS, kRecoQaFile, "hte_Rich_Rich_Electron_Rec_p", "hte_Rich_Rich_Electron_Acc_p", "Ring reco efficiency [%]", 70, 99, 5, nofFilesPartCamRadius);
        DrawMeanEff(iS, kRecoQaFile, "hte_StsRich_StsRich_Electron_Rec_p", "hte_StsRich_StsRich_Electron_Acc_p", "STS-RICH matching efficiency [%]", 70, 99, 5, nofFilesPartCamRadius);

        DrawMeanRms(iS, kGeoTestFile, "fhNofHits_hits", kH1MeanHist, "#hits in ring", 23, 32, 6, nofFilesPartCamTilt);
        DrawMeanRms(iS, kGeoTestFile, "fhBoverAVsMom_points", kH2MeanHist, "B/A (hit fit)", 0.87, 0.96, 6, nofFilesPartCamTilt);
        DrawMeanEff(iS, kGeoTestFile, "fhAccMomEl", "fhMcMomEl", "Geometrical acceptance [%]", 75, 95, 6, nofFilesPartCamTilt);
        DrawMeanRms(iS, kUrqmdTestFile, "fh_nof_hits_per_event", kH1MeanHist, "Nof hits per event", 500, 850, 6, nofFilesPartCamTilt);
        DrawMeanRms(iS, kUrqmdTestFile, "fh_gamma_nontarget_mom", kHEntriesHist, "e^{#pm}_{not target} from #gamma per event", 18, 26, 6, nofFilesPartCamTilt);
        DrawMeanEff(iS, kRecoQaFile, "hte_Rich_Rich_Electron_Rec_p", "hte_Rich_Rich_Electron_Acc_p", "Ring reco efficiency [%]", 70, 99, 6, nofFilesPartCamTilt);
        DrawMeanEff(iS, kRecoQaFile, "hte_StsRich_StsRich_Electron_Rec_p", "hte_StsRich_StsRich_Electron_Acc_p", "STS-RICH matching efficiency [%]", 70, 99, 6, nofFilesPartCamTilt);

        DrawMeanRms2D(iS, kGeoTestFile, "fhNofHits_hits", kH1MeanHist, "#hits in ring", 23, 32);
        DrawMeanRms2D(iS, kGeoTestFile, "fhBoverAVsMom_points", kH2MeanHist, "B/A (point fit)", 0.87, 0.96);
        DrawMeanEff2D(iS, kGeoTestFile, "fhAccMomEl", "fhMcMomEl", "Geometrical acceptance [%]", 82, 92);
        DrawMeanRms2D(iS, kUrqmdTestFile, "fh_nof_hits_per_event", kH1MeanHist, "Nof hits per event", 500, 850);
        DrawMeanRms2D(iS, kUrqmdTestFile, "fh_gamma_nontarget_mom", kHEntriesHist, "e^{#pm}_{not target} from #gamma per event", 18, 26);
        DrawMeanEff2D(iS, kRecoQaFile, "hte_Rich_Rich_Electron_Rec_p", "hte_Rich_Rich_Electron_Acc_p", "Ring reco efficiency [%]", 70, 96);
        DrawMeanEff2D(iS, kRecoQaFile, "hte_StsRich_StsRich_Electron_Rec_p", "hte_StsRich_StsRich_Electron_Acc_p", "STS-RICH matching efficiency [%]", 70, 90);

    }
    fHM->SaveCanvasToImage(fOutputDir);
}


ClassImp(CbmRichGeoTestOpt)