#include<cstdio>
#include<TTree.h>
#include<TSystem.h>
#include<TGenPhaseSpace.h>
#include<TApplication.h>
#include<TFile.h>
#include<TMath.h>
#include<TROOT.h>
#include<TH1.h>
#include<TH2.h>
#include<TCanvas.h>
#include<TStyle.h>
#include<TString.h>
#include<TLatex.h>

void show()
{
  char* t_name = gApplication->Argv(3);

  gStyle->SetPadGridX(1);
  gStyle->SetPadGridY(1);
  gStyle->SetOptStat(0);
  
  TFile* ff1 = new TFile("tree_cu63.root","r");
  ff1->cd();
  char tree_name[20], tree_name_addon[20];
  sprintf(tree_name, "T_%s", t_name);
  sprintf(tree_name_addon, "T_%s_addon", t_name);
  TTree* T1 = (TTree*)gROOT->FindObject(tree_name);
  if(!T1) return;
  T1->AddFriend(tree_name_addon,"tree_cu63_addon.root");

  // flag for p2 existance
  int flag_p2 = 0;
  printf("tree name is %s\n", t_name);
  if(!strcmp(t_name,"bl")||!strcmp(t_name,"bs")) 
    {
      flag_p2 = 1;
      printf("Plotting proton 2 as well.\n");
    }

  TH2F* hp1 = new TH2F("hp1","Recoil Proton Phase Space; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);
  TH2F* hp2 = new TH2F("hp2","Decay Proton Phase Space; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);
  TH2F* hpi = new TH2F("hpi","Pion Phase Space; Momentum [GeV]; Angle [degree]",50,0,0.5,60,0,180);
  TH2F* hk = new TH2F("hk","Kaon Phase Space; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);

  TH2F* hp1d = new TH2F("hp1d","Recoil Proton Detection; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);
  TH2F* hp2d = new TH2F("hp2d","Decay Proton Detection; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);
  TH2F* hpid = new TH2F("hpid","Pion Detection; Momentum [GeV]; Angle [degree]",50,0,0.5,60,0,180);
  TH2F* hkd= new TH2F("hkd","Kaon Detection; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);

  TH2F* hp1a = new TH2F("hp1a","Recoil Proton Angular Detection; Phi Angle [degree]; Theta Angle [degree]",100,-200,200,60,0,180);
  TH2F* hp2a = new TH2F("hp2a","Decay Proton Angular Detection; Phi Angle [degree]; Theta Angle [degree]",100,-200,200,60,0,180);
  TH2F* hpia = new TH2F("hpia","Pion Angular Detection; Phi Angle [degree]; Theta Angle [degree]",100,-200,200,60,0,180);
  TH2F* hka = new TH2F("hka","Kaon Angular Detection; Phi Angle [degree]; Theta Angle [degree]",100,-200,200,60,0,180);

  TString cut, cut_det;
  if(!strcmp(t_name,"bl"))
    {
      cut = "w.pf*w.nf*w.lk*w.ppi";
    }
  else if(!strcmp(t_name,"dl")) 
    {
      cut = "w.lk*w.ppi";
    }
  else if(!strcmp(t_name,"bs"))
    {
      cut = "w.pf*w.nf*w.sk*w.lg*w.ppi";
    }
  else if(!strcmp(t_name,"ds"))
    {
      cut = "w.sk*w.lg*w.ppi";
    }

  TCanvas* c1;
  if(flag_p2)
    {
      c1 = new TCanvas("c1","Phase Space",900,1000);
      c1->Divide(3,4);
    }
  else
    {
      c1 = new TCanvas("c1","Phase Space",900,750);
      c1->Divide(3,3);
    }

  double pxmax, pxmin, pymax, pymin;
  char ss[100];
  TLatex ll;
  double eff;
  
  cut_det = cut+"*pi.det";
  c1->cd(1);
  T1->Draw("pi.af.Th*180/TMath::Pi():pi.af.P>>hpi",cut,"cont");
  c1->cd(2);
  T1->Draw("pi.af.Th*180/TMath::Pi():pi.af.P>>hpid",cut_det,"cont");
  gPad->Update();
  pxmax=gPad->GetUxmax();
  pxmin=gPad->GetUxmin();
  pymax=gPad->GetUymax();
  pymin=gPad->GetUymin();
  eff = hpid->Integral()/hpi->Integral();
  sprintf(ss, "Efficiency: %.2f", eff);
  ll.DrawLatex(pxmin+(pxmax-pxmin)*0.5,pymin+(pymax-pymin)*0.7,ss);
  c1->cd(3);
  T1->Draw("pi.af.Th*180/TMath::Pi():pi.af.Ph*180/TMath::Pi()>>hpia",cut_det,"cont");
  c1->Update();
  c1->cd(2);

  cut_det = cut+"*k.det";
  c1->cd(4);
  T1->Draw("k.af.Th*180/TMath::Pi():k.af.P>>hk",cut,"cont");
  c1->cd(5);
  T1->Draw("k.af.Th*180/TMath::Pi():k.af.P>>hkd",cut_det,"cont");
  gPad->Update();
  pxmax=gPad->GetUxmax();
  pxmin=gPad->GetUxmin();
  pymax=gPad->GetUymax();
  pymin=gPad->GetUymin();
  eff = hkd->Integral()/hk->Integral();
  sprintf(ss, "Efficiency: %.2f", eff);
  ll.DrawLatex(pxmin+(pxmax-pxmin)*0.5,pymin+(pymax-pymin)*0.7,ss);
  c1->cd(6);
  T1->Draw("k.af.Th*180/TMath::Pi():k.af.Ph*180/TMath::Pi()>>hka",cut_det,"cont");
  c1->Update();

  cut_det = cut+"*p1.det";
  c1->cd(7);
  T1->Draw("p1.af.Th*180/TMath::Pi():p1.af.P>>hp1",cut,"cont");
  c1->cd(8);
  T1->Draw("p1.af.Th*180/TMath::Pi():p1.af.P>>hp1d",cut_det,"cont");
  gPad->Update();
  pxmax=gPad->GetUxmax();
  pxmin=gPad->GetUxmin();
  pymax=gPad->GetUymax();
  pymin=gPad->GetUymin();
  eff = hp1d->Integral()/hp1->Integral();
  sprintf(ss, "Efficiency: %.2f", eff);
  ll.DrawLatex(pxmin+(pxmax-pxmin)*0.5,pymin+(pymax-pymin)*0.7,ss);
  c1->cd(9);
  T1->Draw("p1.af.Th*180/TMath::Pi():p1.af.Ph*180/TMath::Pi()>>hp1a",cut_det,"cont");
  c1->Update();

  if(flag_p2)
    {
      cut_det = cut+"*p2.det";
      c1->cd(10);
      T1->Draw("p2.af.Th*180/TMath::Pi():p2.af.P>>hp2",cut,"cont");
      c1->cd(11);
      T1->Draw("p2.af.Th*180/TMath::Pi():p2.af.P>>hp2d",cut_det,"cont");
      gPad->Update();
      pxmax=gPad->GetUxmax();
      pxmin=gPad->GetUxmin();
      pymax=gPad->GetUymax();
      pymin=gPad->GetUymin();
      eff = hp2d->Integral()/hp2->Integral();
      sprintf(ss, "Efficiency: %.2f", eff);
      ll.DrawLatex(pxmin+(pxmax-pxmin)*0.5,pymin+(pymax-pymin)*0.7,ss);
      c1->cd(12);
      T1->Draw("p2.af.Th*180/TMath::Pi():p2.af.Ph*180/TMath::Pi()>>hp2a",cut_det,"cont");
      c1->Update();
    }
  char plot_name[200];
  sprintf(plot_name,"./plot/phase_%s.png",t_name);
  c1->Print(plot_name);
  sprintf(plot_name,"./plot/phase_%s.eps",t_name);
  c1->Print(plot_name);

  TH1F* hm1 = new TH1F("hm1","Invariant mass with recoil proton; W (GeV/c2)",500,1.4,2.4);
  TH1F* hm1d = new TH1F("hm1d","Invariant mass with recoil proton; W (GeV/c2)",500,1.4,2.4);
  TH1F* hm2 = new TH1F("hm2","Invariant mass with decay proton; W (GeV/c2)",500,1.4,2.4);
  TH1F* hm2d = new TH1F("hm2d","Invariant mass with decay proton; W (GeV/c2)",500,1.4,2.4);

  hm1->SetLineColor(4);
  hm1->SetFillColor(4);
  hm1->SetFillStyle(3002);
  hm1d->SetLineColor(2);
  hm1d->SetFillColor(2);
  hm1d->SetFillStyle(3001);

  hm2->SetLineColor(4);
  hm2->SetFillColor(4);
  hm2->SetFillStyle(3002);
  hm2d->SetLineColor(2);
  hm2d->SetFillColor(2);
  hm2d->SetFillStyle(3001);

  TCanvas* c2;
  if(flag_p2)
    {
      c2 = new TCanvas("c2","Invariant Mass",600,800);
      c2->Divide(1,2);
    }
  else c2 = new TCanvas("c2","Invariant Mass",600,400);

  cut_det = cut+"*p1.det*k.det*pi.det";
  c2->cd(1);
  T1->Draw("M1>>hm1",cut);
  T1->Draw("M1>>hm1d",cut_det,"same");
  gPad->Update();
  pxmax=gPad->GetUxmax();
  pxmin=gPad->GetUxmin();
  pymax=gPad->GetUymax();
  pymin=gPad->GetUymin();
  eff = hm1d->Integral()/hm1->Integral();
  sprintf(ss, "Efficiency: %.4f", eff);
  ll.DrawLatex(pxmin+(pxmax-pxmin)*0.05,pymin+(pymax-pymin)*0.85,ss);
  if(flag_p2)
    {
      cut_det = cut+"*p2.det*k.det*pi.det";
      c2->cd(2);
      T1->Draw("M2>>hm2",cut);
      T1->Draw("M2>>hm2d",cut_det,"same");
      gPad->Update();
      pxmax=gPad->GetUxmax();
      pxmin=gPad->GetUxmin();
      pymax=gPad->GetUymax();
      pymin=gPad->GetUymin();
      eff = hm2d->Integral()/hm2->Integral();
      sprintf(ss, "Efficiency: %.4f", eff);
      ll.DrawLatex(pxmin+(pxmax-pxmin)*0.05,pymin+(pymax-pymin)*0.85,ss);
    }
  sprintf(plot_name,"./plot/inv_%s.png",t_name);
  c2->Print(plot_name);
  sprintf(plot_name,"./plot/inv_%s.eps",t_name);
  c2->Print(plot_name);
}