#include <iostream>

#include <TString.h>
#include <TStyle.h>
#include <TMath.h>
#include <TFile.h>
#include <TList.h>
#include <TH3.h>
#include <TF1.h>
#include <TLatex.h>
#include <TGAxis.h>
#include <TPad.h>


namespace
{

//_______________________________________________________________________
TH3 * clone_sizes(const TH3 * src)
{
    TString name(src->GetName() + TString("_clone"));
    //
    const double xMin = src->GetXaxis()->GetXmin();
    const double xMax = src->GetXaxis()->GetXmax();
    const int xBins = src->GetNbinsX();
    //
    const double yMin = src->GetYaxis()->GetXmin();
    const double yMax = src->GetYaxis()->GetXmax();
    const int yBins = src->GetNbinsY();
    //
    const double zMin = src->GetZaxis()->GetXmin();
    const double zMax = src->GetZaxis()->GetXmax();
    const int zBins = src->GetNbinsZ();

    return new TH3D(name.Data(), name.Data(), xBins, xMin, xMax, yBins, yMin, yMax, zBins, zMin, zMax);
}

//_______________________________________________________________________
TH3 * transform_hist(const TH3 * src, double (*map)(double))
{
    TH3 * dst = clone_sizes(src);

    for(int i = 1; i <= src->GetNbinsX(); ++i)
    {
        for(int j = 1; j <= src->GetNbinsY(); ++j)
        {
            for(int k = 1; k <= src->GetNbinsZ(); ++k)
            {
                const double val = src->GetBinContent(i, j, k);
                dst->SetBinContent(i, j, k, map(val));
            }
        }
    }

    return dst;
}

//_______________________________________________________________________
double convert_energy(double val)
{
  double val_Joules = val*1.6e-16;
    return val_Joules;
}

}

//_______________________________________________________________________
void png_issue()
{
    TFile input("rel_energy.root", "READ");
    TH3 * energyHist = dynamic_cast<TH3*>(input.Get("energy_hist"));
    if(!energyHist)
    {
        std::cout<<"Could not find energy_hist or type is not TH3 derived\n";
        return;
    }

    gStyle->SetCanvasPreferGL(1);
    TH3 * energyHistModified = transform_hist(energyHist, convert_energy); 


 
    double density = 8.91e-6;  // kg/mm3
    double specific_heat = 440.;  // J/kgK
    double bin_volume = 1.; // 1 mm3

    for(int i = 1; i <= energyHistModified->GetNbinsX(); ++i)
    {
        for(int j = 1; j <= energyHistModified->GetNbinsY(); ++j)
        {
            for(int k = 1; k <= energyHistModified->GetNbinsZ(); ++k)
            {
                const double val = energyHistModified->GetBinContent(i, j, k);
                const double new_val =  val / (density*bin_volume*specific_heat);
                energyHistModified->SetBinContent(i, j, k, new_val);
            }
        }
    }

    ///// finished conversion ////////

 energyHistModified->Scale(1./(150.*5.e7));

    gStyle->SetPalette(1);
    energyHistModified->SetDirectory(0);
    energyHistModified->GetXaxis()->SetTitle("X (mm)");
    energyHistModified->GetYaxis()->SetTitle("Y (mm)");
    energyHistModified->GetZaxis()->SetTitle("Z (mm)");
    energyHistModified->GetXaxis()->SetTitleOffset(2.);
    energyHistModified->GetYaxis()->SetTitleOffset(2.);
    energyHistModified->GetZaxis()->SetTitleOffset(1.5);
    energyHistModified->GetXaxis()->CenterTitle();
    energyHistModified->GetYaxis()->CenterTitle();
    energyHistModified->GetZaxis()->CenterTitle();
    energyHistModified->Draw("glcolz");

  TLatex * text = new TLatex(0.915,0.85,"#DeltaT (K)");
  text->SetNDC(kTRUE);
  text->Draw();

}