#include <stdio.h>
#include <iostream>

#include "TLorentzVector.h"
#include "TGenPhaseSpace.h"
#include "TF1.h"
#include "TString.h"
#include "TH2D.h"
#include "TROOT.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TProfile.h"
#include "TGraph.h"
#include "TGaxis.h"
#include "TLegend.h"
#include "TMultiGraph.h"

#include "TStyle.h"
#include "TLatex.h"
#include "TGraphPolargram.h"
#include "TPaletteAxis.h"
#include <getopt.h>
#include <glob.h>
#include <fstream>
#include "TGraphPolar.h"

int Test()
{

  double p0 = 25;
  const char* fileOut = "bla.pdf";
 
  //gROOT->ProcessLine(".L ~/.rootlogon.C");
  gStyle->SetLabelSize(0.06,"xyz");
  gStyle->SetLabelOffset(0.012,"xy");
  gStyle->SetTitleOffset(1.05,"x");
  gStyle->SetTitleOffset(1.15,"y");
  gStyle->SetTitleSize(0.065,"xy");
  gStyle->SetTitleBorderSize(0);
  gStyle->SetOptStat("");

  TLorentzVector* p4Pi = new TLorentzVector();
  double mPi = 0.13957;
  p4Pi->SetXYZM(0,0,p0,mPi);

  double binW = 200/2595.;
  double thMax =1.5;
  TH2D hPol("hPol", "",40,-TMath::Pi(),TMath::Pi(),90,0,thMax);
  TH1D hTh("hTh"  , ";#theta_{#nu#pi} [#circ]; evts per sr",90,0,thMax);

  hPol.GetZaxis()->SetMaxDigits(3);
  hTh.GetYaxis()->SetMaxDigits(3);

  double bW = hTh.GetBinWidth(2)*TMath::DegToRad();
  double bWPhi = 2*TMath::Pi()/40;

  TCanvas cv;
  double mMu = 0.10565837;
  double massDecay[2] = {mMu,0.};
  TGenPhaseSpace generator;

  int nP=100;

  for (int ip(0); ip < nP; ip++) {

    int nPh = int(gRandom->Poisson(1000.));
    p4Pi->SetXYZM(0,0,p0,mPi);
    generator.SetDecay(*p4Pi, 2, massDecay);
  
      for (int i(0); i<nPh ; i++) {
	// generate
	Double_t weight = generator.Generate();

	// true
	TLorentzVector* p4Mu = generator.GetDecay(0);
	TLorentzVector* p4Nu = generator.GetDecay(1);
	double th  = p4Nu->Vect().Theta();
	double phi = p4Nu->Vect().Phi();

	hPol.Fill(phi,th*TMath::RadToDeg(),weight/nP / bW / bWPhi);
	//hPol.Fill(phi,th*TMath::RadToDeg(),weight/nP);
	hTh.Fill(th*TMath::RadToDeg()     ,weight/nP / bW / (2*TMath::Pi()));
      }
      
    }
  


  // plot_histogram_in_polar(&hPol)->SaveAs("test.png");
  hPol.Draw("lego2 polz");
  //  cv.SetLeftMargin(0.15);
  // cv.SetRightMargin(0.15);

  hPol.SetLineColorAlpha(kBlack, 0.0);
  cv.SetGrid(0,0);
  //  cv.SetRightMargin(0.15);
  cv.SetTheta(90.);
  cv.SetPhi(0.);
  gPad->Update();

  TGraphPolargram * gp = new TGraphPolargram("g",0,thMax,                                             -180., 180.);
  gp->SetToDegree ();
  gp->SetNdivPolar(8);
  gp->SetNdivRadial(4);
  gp->SetLineColor(1);
  gp->Draw();

  gPad->Update();
  TPaletteAxis *palette = (TPaletteAxis*)hPol.GetListOfFunctions()->FindObject("palette");

// the following lines moe the paletter. Choose the values you need for the position.
   palette->SetX1NDC(0.90);
   palette->SetX2NDC(0.92);
   palette->SetY1NDC(0.2);
   palette->SetY2NDC(0.8);

   gPad->Modified();
   gPad->Update();

   TString bOut(fileOut);
   TString fOut1(bOut);
   fOut1.Remove(fOut1.Last('.'));

   cv.SaveAs(fOut1+"Polar.pdf");
   cv.SaveAs(fOut1+"Polar.png");

  cv.SetLeftMargin(0.15);
  cv.SetBottomMargin(0.15);
  hTh.Draw("hist");
  TString fOut2(bOut);
  fOut2.Remove(fOut2.Last('.'));
  cv.SaveAs(fOut2+"SR.pdf");

  std::ofstream ofs(fOut2+"SR.txt");
  for (int i(0); i < hTh.GetXaxis()->GetNbins(); i++) {
    ofs<<hTh.GetXaxis()->GetBinCenter(i+1)<<" "<<hTh.GetBinContent(i+1)<<std::endl;
  }
  ofs.close();
  double reso = 4;
  hTh.Scale(TMath::Pi()*reso*1e-3*reso*1e-3);
  hTh.GetYaxis()->SetTitle(Form("evt in a %1.0f mard cone",reso));
  hTh.Draw("hist");
  TString fOut3(bOut);
  fOut3.Remove(fOut3.Last('.'));
  cv.SaveAs(fOut3+"Mrad.pdf");

  std::ofstream ofs2(fOut3+"Mrad.txt");
  for (int i(0); i < hTh.GetXaxis()->GetNbins(); i++) {
    ofs2<<hTh.GetXaxis()->GetBinCenter(i+1)<<" "<<hTh.GetBinContent(i+1)<<std::endl;
  }
  ofs2.close();



  return 0;
}