#include <math.h>

void weights(char* name="weights_mu800_900mev")
{
  gROOT->Reset();
  gStyle->SetOptStat(0000);
  gStyle->SetOptFit(111);
  gStyle->SetCanvasBorderMode(-1);
  gStyle->SetCanvasBorderSize(1);
  gStyle->SetCanvasColor(10);
  gStyle->SetTitleOffset(1.20,"Y"); // move the title

  TString inputname(name);
  inputname+=".root";
  TFile f(inputname);
  TTree *tree = (TTree*)f.Get("W");

  double ch,pl;
  tree->SetBranchAddress("ch",&ch);
  tree->SetBranchAddress("pl",&pl);

  double xlow=0.025;
  double xup=0.11;
  TH2F *h1 = new TH2F("h1","",100,xlow,xup,100,0.,0.8);
  TH2F *h2 = new TH2F("h2","",100,xlow,xup,100,0.,3.);

  int nbins=6;
  TProfile *prof1 = new TProfile("p1","",nbins,xlow,xup);
  TProfile *prof2 = new TProfile("p2","",nbins,xlow,xup);

  for (unsigned int i = 0; i < tree->GetEntries(); i++) {
    tree->GetEntry(i);
    h1->Fill(pl,ch);
    prof1->Fill(pl,ch);
    if (ch>0.) h2->Fill(pl,-log(ch));
    if (ch>0.) prof2->Fill(pl,-log(ch));
  }
  prof1->SetErrorOption("s"); // to have spread, and not spread/sqrt(N), as error.
  prof2->SetErrorOption("s"); // to have spread, and not spread/sqrt(N), as error.

  // graph of spreads:
  double binwidth=(xup-xlow)/nbins;
  double* x = NULL; 
  double* y1 = NULL; 
  double* y2 = NULL; 
  x = new double[nbins];
  y1 = new double[nbins];
  y2 = new double[nbins];
  for (unsigned int i=0; i<nbins; i++) {
    x[i]=0.;
    y1[i]=0.;
    y2[i]=0.;
  }
  TGraph* gr1 = new TGraph(nbins,x,y1);
  TGraph* gr2 = new TGraph(nbins,x,y2);
  for (unsigned int i=0; i<nbins; i++) {
    x[i]=xlow+binwidth*i;
    if (prof1->GetBinContent(i+1)>0.) y1[i]=prof1->GetBinContent(i+1)/prof1->GetBinError(i+1);
    if (prof2->GetBinContent(i+1)>0.) y2[i]=prof2->GetBinContent(i+1)/prof2->GetBinError(i+1);
    gr1->SetPoint(i,x[i],y1[i]);
    gr2->SetPoint(i,x[i],y2[i]);
  }

  
  // plot:
  TCanvas* myCanvas = new TCanvas("myCanvas","",200,10,900,800);
  myCanvas->Divide(3,2);
  
  h1->GetXaxis()->SetTitle("path lenght [cm]");
  h2->GetXaxis()->SetTitle("path lenght [cm]");
  prof1->GetXaxis()->SetTitle("path lenght [cm]");
  prof2->GetXaxis()->SetTitle("path lenght [cm]");
  gr1->GetXaxis()->SetTitle("path lenght [cm]");
  gr2->GetXaxis()->SetTitle("path lenght [cm]");

  h1->GetYaxis()->SetTitle("charge [MeV]");
  prof1->GetYaxis()->SetTitle("charge [MeV]");
  h2->GetYaxis()->SetTitle("-log(charge) [MeV]");
  prof2->GetYaxis()->SetTitle("-log(charge) [MeV]");
  gr1->GetYaxis()->SetTitle("mean/RMS");
  gr2->GetYaxis()->SetTitle("mean/RMS");

  h1->SetMarkerStyle(23);
  h2->SetMarkerStyle(23);
  h1->SetMarkerSize(0.5);
  h2->SetMarkerSize(0.5);

  gr1->SetMarkerStyle(21);
  gr2->SetMarkerStyle(21);
  gr1->SetTitle("");
  gr2->SetTitle("");

  myCanvas->cd(1);
  h1->DrawCopy();

  myCanvas->cd(2);
  prof1->DrawCopy();

  myCanvas->cd(3);
  gr1->Draw("AP");
  // fit:
  gr1->Fit("pol1");
  TF1 *fit = gr1->GetFunction("pol1");
  unsigned ndecimals=2;
  Double_t p0 = fround(fit->GetParameter(0),ndecimals);
  Double_t p1 = fround(fit->GetParameter(1),ndecimals);
  cout << "--------------------" << endl;
  cout << "Result of the fit:" << endl;
  cout << "--------------------" << endl;
  cout << "p0 = " << p0 << endl;
  cout << "p1 = " << p1 << endl;
  cout << "--------------------" << endl;
  // save results in ascii:
  TString ascii("weights/");
  ascii+=name;
  ascii+=".txt";
  fstream f1;
  f1.open(ascii,ios::out);
  f1 << p0 << endl;
  f1 << p1 << endl;
  f1.close();
  cout << "Fit results written into: " << ascii << endl;


  myCanvas->cd(4);
  h2->DrawCopy();

  myCanvas->cd(5);
  prof2->DrawCopy();

  myCanvas->cd(6);
  gr2->Draw("AP");

  TString gif("weights/");
  gif+=name;
  gif+=".png";
  myCanvas->Print(gif);


  delete [] x;
  delete [] y1;
  delete [] y2;

  delete h1;
  delete h2;
  delete prof1;
  delete prof2;
  delete gr1;
  delete gr2;
  delete myCanvas;
}

/* round number n to d decimal points */
double fround(double n, unsigned d)
{
      return floor(n * pow(10., d) + .5) / pow(10., d);
}