#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_pkk()
{
  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,"b")) 
    {
      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* hk1 = new TH2F("hk1","Kaon+ Phase Space; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);
  TH2F* hk2 = new TH2F("hk2","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* hk1d = new TH2F("hk1d","Kaon+ Detection; Momentum [GeV]; Angle [degree]",50,0,1.5,60,0,180);
  TH2F* hk2d = new TH2F("hk2d","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* hk1a = new TH2F("hk1a","Kaon+ Angular Detection; Phi Angle [degree]; Theta Angle [degree]",100,-200,200,60,0,180);
  TH2F* hk2a = new TH2F("hk2a","Kaon- Angular Detection; Phi Angle [degree]; Theta Angle [degree]",100,-200,200,60,0,180);

  TString cut, cut_det;
  if(!strcmp(t_name,"b"))
    {
      cut = "w.pf*w.nf*w.pkk";
    }
  else if(!strcmp(t_name,"d")) 
    {
      cut = "w.pkk";
    }
  else if(!strcmp(t_name,"nb"))
    {
      cut = "w.pf*w.kk";
    }
  else if(!strcmp(t_name,"af"))
    {
      cut = "w.paf*w.kk";
   } 
  else if(!strcmp(t_name,"l"))
    {
      cut = "w.lk*w.pk";
    }

  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;
  
  // kaon+
  cut_det = cut+"*k1.det";
  c1->cd(1);
  T1->Draw("k1.af.Th*180/TMath::Pi():k1.af.P>>hk1",cut,"cont");
  c1->cd(2);
  T1->Draw("k1.af.Th*180/TMath::Pi():k1.af.P>>hk1d",cut_det,"cont");
  gPad->Update();
  pxmax=gPad->GetUxmax();
  pxmin=gPad->GetUxmin();
  pymax=gPad->GetUymax();
  pymin=gPad->GetUymin();
  eff = hk1d->Integral()/hk1->Integral();
  sprintf(ss, "Efficiency: %.2f", eff);
  ll.DrawLatex(pxmin+(pxmax-pxmin)*0.5,pymin+(pymax-pymin)*0.7,ss);
  c1->cd(3);
  T1->Draw("k1.af.Th*180/TMath::Pi():k1.af.Ph*180/TMath::Pi()>>hk1a",cut_det,"cont");
  c1->Update();
  c1->cd(2);

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

  // proton 1
  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);
  //  c2->SetLogy(1);
  cut_det = cut+"*p1.det*k1.det*k2.det";
  c2->cd(1);
  if(flag_p2)
    {
      c2_1->SetLogy(1);
    }
  else c2->SetLogy(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);
  ll.DrawLatex(pxmin+(pxmax-pxmin)*0.05,pymin+(pymax-pymin)*50,ss);
  if(flag_p2)
    {
      cut_det = cut+"*p2.det*k1.det*k2.det";
      c2->cd(2);
      c2_2->SetLogy(1);
      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);
      ll.DrawLatex(pxmin+(pxmax-pxmin)*0.05,pymin+(pymax-pymin)*50,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);
}