#include "TH1.h"
#include "TH2.h"
#include "THnSparse.h"
#include "TAxis.h"
#include "TString.h"

// Usage example: TH1D *h1 = Project1D(his3D, 451.5, 670.5, 2.0);
// Note: in principle, you should set eg_low <= eg_hi
TH1D *Project1D(THnSparse *h, // the ("totally unsymmetric") "cubic" histogram
                Double_t eg_low = 451.5, // E_gamma1
                Double_t eg_hi = 670.5, // E_gamma2
		Double_t w = 2.0,
                Option_t *option = "") { // THnSparse::Projection option

  Double_t low1 = eg_low - w;
  Double_t up1  = eg_low + w;

  Double_t low2 = eg_hi - w;
  Double_t up2  = eg_hi + w;

  if (!h) return 0; // just a precaution
  TAxis *a0 = h->GetAxis(0);
  TAxis *a1 = h->GetAxis(1);
  TAxis *a2 = h->GetAxis(2);
  a0->SetRange(a0->FindFixBin( low1 ), a0->FindFixBin( up1 ));
  a1->SetRange(a1->FindFixBin( low2 ), a1->FindFixBin( up2 ));
  a2->SetRange(-1, -1);
  TH1D *h1 = h->Projection(2, option);
  h1->SetName(TString::Format("%s_Project1D", h->GetName()));
  a0->SetRange(a0->FindFixBin( low2 ), a0->FindFixBin( up2 ));
  a1->SetRange(a1->FindFixBin( low1 ), a1->FindFixBin( up1 ));
  TH1D *hi = h->Projection(2, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(a0->FindFixBin( low1 ), a0->FindFixBin( up1 ));
  a1->SetRange(-1, -1);
  a2->SetRange(a2->FindFixBin( low2 ), a2->FindFixBin( up2 ));
  hi = h->Projection(1, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(a0->FindFixBin( low2 ), a0->FindFixBin( up2 ));
  a2->SetRange(a2->FindFixBin( low1 ), a2->FindFixBin( up1 ));
  hi = h->Projection(1, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(a1->FindFixBin( low1 ), a1->FindFixBin( up1 ));
  a2->SetRange(a2->FindFixBin( low2 ), a2->FindFixBin( up2 ));
  hi = h->Projection(0, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a1->SetRange(a1->FindFixBin( low2 ), a1->FindFixBin( up2 ));
  a2->SetRange(a2->FindFixBin( low1 ), a2->FindFixBin( up1 ));
  hi = h->Projection(0, option);
  h1->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(-1, -1);
  a2->SetRange(-1, -1);

  h1->SetName(TString::Format("his_%d_%d", Int_t(eg_low), Int_t(eg_hi)));

  return h1;
}

// Usage example: TH2D *h2 = Project2D(his3D, 670.5, 2.0);
// Note: in principle, you should set low <= up
TH2D *Project2D(THnSparse *h, // the ("totally unsymmetric") "cubic" histogram
                Double_t eg = 670.5, // E_gamma
		Double_t w = 2.0,
                Option_t *option = "") { // THnSparse::Projection option
  if (!h) return 0; // just a precaution

  Double_t low = eg - w;
  Double_t up  = eg + w;

  TAxis *a0 = h->GetAxis(0);
  TAxis *a1 = h->GetAxis(1);
  TAxis *a2 = h->GetAxis(2);
  a0->SetRange(a0->FindFixBin( low ), a0->FindFixBin( up ));
  a1->SetRange(-1, -1);
  a2->SetRange(-1, -1);
  TH2D *h2 = h->Projection(1, 2, option);
  h2->SetName(TString::Format("%s_Project2D", h->GetName()));
  TH2D *hi = h->Projection(2, 1, option);
  h2->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(a1->FindFixBin( low ), a1->FindFixBin( up ));
  a2->SetRange(-1, -1);
  hi = h->Projection(0, 2, option);
  h2->Add(hi);
  delete hi; // no longer needed
  hi = h->Projection(2, 0, option);
  h2->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(-1, -1);
  a2->SetRange(a2->FindFixBin( low ), a2->FindFixBin( up ));
  hi = h->Projection(0, 1, option);
  h2->Add(hi);
  delete hi; // no longer needed
  hi = h->Projection(1, 0, option);
  h2->Add(hi);
  delete hi; // no longer needed
  a0->SetRange(-1, -1);
  a1->SetRange(-1, -1);
  a2->SetRange(-1, -1);

  h2->SetName(TString::Format("his2D_%d", Int_t(eg)));

  return h2;
}