const int MIN_BCAL_HITS = 1;
const int MIN_FCAL_HITS = 1;
const int MIN_TOF_HITS = 1;
const int MIN_CDC_HITS = 8;
const int MIN_FDC_HITS = 8;

const double mom_min = 0;
const double mom_max = 9;
const int mom_nbin = 90;

const double ang_min = 0;
const double ang_max = 90;
const double ang_nbin = 90;

void makeplot(TTree* T, char *name, char *title, TCanvas* c1, int i, int j,
		int k, bool charge = 1) {
	char hname[100];
	char htitle[200];
	char plotstr[500];
	char cutstr[500];

	/*
	 * Phase space plot
	 */
	c1->cd(i);
	gPad->SetMargin(0.13,0.13,0.13,0.07);
	sprintf(hname, "%s_2d", name);
	sprintf(htitle, "%s distribution;Angle [deg];Momentum [GeV]", name);
	TH2I* phase_histo = new TH2I(hname, htitle, ang_nbin, ang_min, ang_max,
			mom_nbin, mom_min, mom_max);

	sprintf(plotstr, "%s.P:%s.Th*180./TMath::Pi()>>%s", name, name, hname);
	T->Draw(plotstr, "", "colz");

	/*
	 * Draw a Kaon PID box
	 */
	TBox *b1 = new TBox(1,0,10,2);
	b1->SetLineWidth(2);
	b1->SetLineColor(kOrange);
	b1->SetFillStyle(0);
	b1->Draw();

	/*
	 * Draw label
	 */
	TLatex ll;
	ll.SetTextAlign(22);
	ll.SetTextSize(0.06);
	ll.SetNDC(true);
	ll.DrawLatex(0.5, 0.96, title);

	/*
	 * Momentum acceptance
	 */
	c1->cd(j);
	gPad->SetMargin(0.13,0.07,0.13,0.07);
	sprintf(hname, "%s_ph", name);
	sprintf(htitle, "%s momentum distribution;Momentum [GeV]", name);
	TH1I* mom_dist = new TH1I(hname, htitle, mom_nbin, mom_min, mom_max);
	mom_dist->SetLineColor(kBlack);

	sprintf(plotstr, "%s.P>>%s", name, hname);
	T->Draw(plotstr, "", "");

	sprintf(hname, "%s_pch", name);
	sprintf(htitle, "%s momentum acceptance in BCAL;Momentum [GeV]", name);
	TH1I* mom_accp_bcal = new TH1I(hname, htitle, mom_nbin, mom_min, mom_max);
	mom_accp_bcal->SetLineColor(kRed);
	mom_accp_bcal->SetFillColor(kRed);
	mom_accp_bcal->SetFillStyle(3004);

	sprintf(plotstr, "%s.P>>%s", name, hname);
	if (charge)
		sprintf(cutstr,
				"(%s.hits_cdc>=%d||%s.hits_fdc>=%d)&&(%s.hits_bcal>=%d)", name,
				MIN_CDC_HITS, name, MIN_FDC_HITS, name, MIN_BCAL_HITS);
	else
		sprintf(cutstr, "%s.hits_bcal>=%d", name, MIN_BCAL_HITS);
	T->Draw(plotstr, cutstr, "same");

	sprintf(hname, "%s_pfh", name);
	sprintf(htitle, "%s momentum acceptance in TOF;Momentum [GeV]", name);
	TH1I* mom_accp_tof = new TH1I(hname, htitle, mom_nbin, mom_min, mom_max);
	mom_accp_tof->SetLineColor(kBlue);
	mom_accp_tof->SetFillColor(kBlue);
	mom_accp_tof->SetFillStyle(3005);

	sprintf(plotstr, "%s.P>>%s", name, hname);
	if (charge)
		sprintf(cutstr,
				"(%s.hits_cdc>=%d||%s.hits_fdc>=%d)&&(%s.hits_tof>=%d)&&(%s.hits_fcal>=%d)",
				name, MIN_CDC_HITS, name, MIN_FDC_HITS, name, MIN_TOF_HITS,
				name, MIN_FCAL_HITS);
	else
		sprintf(cutstr, "%s.hits_fcal>=%d", name, MIN_FCAL_HITS);
	T->Draw(plotstr, cutstr, "same");

	/*
	 * Draw Label
	 */
	ll.SetTextColor(kBlue);
	ll.DrawLatex(0.43, 0.96, "BCAL");
	ll.SetTextColor(kRed);
	ll.DrawLatex(0.63, 0.96, "TOF");

	/*
	 * Angular acceptance
	 */
	c1->cd(k);
	gPad->SetMargin(0.13,0.07,0.13,0.07);
	sprintf(hname, "%s_ah", name);
	sprintf(htitle, "%s angular distribution;Angle [deg]", name);
	TH1I* angle_dist = new TH1I(hname, htitle, ang_nbin, ang_min, ang_max);
	angle_dist->SetLineColor(kBlack);

	sprintf(plotstr, "%s.Th*180./TMath::Pi()>>%s", name, hname);
	T->Draw(plotstr, "", "");

	sprintf(hname, "%s_ach", name);
	sprintf(htitle, "%s angular acceptance in BCAL;Angle [deg]", name);
	TH1I* angle_accp_bcal = new TH1I(hname, htitle, ang_nbin, ang_min, ang_max);
	angle_accp_bcal->SetLineColor(kRed);
	angle_accp_bcal->SetFillColor(kRed);
	angle_accp_bcal->SetFillStyle(3004);

	sprintf(plotstr, "%s.Th*180./TMath::Pi()>>%s", name, hname);
	if (charge)
		sprintf(cutstr,
				"(%s.hits_cdc>=%d||%s.hits_fdc>=%d)&&(%s.hits_bcal>=%d)", name,
				MIN_CDC_HITS, name, MIN_FDC_HITS, name, MIN_BCAL_HITS);
	else
		sprintf(cutstr, "%s.hits_bcal>=%d", name, MIN_BCAL_HITS);
	T->Draw(plotstr, cutstr, "same");

	sprintf(hname, "%s_afh", name);
	sprintf(htitle, "%s angular acceptance in TOF;Angle [deg]", name);
	TH1I* angle_accp_tof = new TH1I(hname, htitle, ang_nbin, ang_min, ang_max);
	angle_accp_tof->SetLineColor(kBlue);
	angle_accp_tof->SetFillColor(kBlue);
	angle_accp_tof->SetFillStyle(3005);

	sprintf(plotstr, "%s.Th*180./TMath::Pi()>>%s", name, hname);
	if (charge)
		sprintf(cutstr,
				"(%s.hits_cdc>=%d||%s.hits_fdc>=%d)&&(%s.hits_tof>=%d)&&(%s.hits_fcal>=%d)",
				name, MIN_CDC_HITS, name, MIN_FDC_HITS, name, MIN_TOF_HITS,
				name, MIN_FCAL_HITS);
	else
		sprintf(cutstr, "%s.hits_fcal>=%d", name, MIN_FCAL_HITS);
	T->Draw(plotstr, cutstr, "same");
	/*
	 * Draw Label
	 */
	ll.SetTextColor(kBlue);
	ll.DrawLatex(0.43, 0.96, "BCAL");
	ll.SetTextColor(kRed);
	ll.DrawLatex(0.63, 0.96, "TOF");
}

void corrplot(TTree* T, char *name1, char *name2, TCanvas* c1, int i, int j) {
	char hname[100];
	char htitle[200];
	char plotstr[500];
	char cutstr[500];

	/*
	 * Angular plot
	 */
	c1->cd(i);
	gPad->SetMargin(0.13,0.13,0.13,0.07);
	sprintf(hname, "ang_%s_vs_%s", name1, name2);
	sprintf(htitle, "Angular correlation;Angle of %s [deg];Angle of %s [deg]",
			name2, name1);
	TH2I* phase_histo = new TH2I(hname, htitle, ang_nbin, ang_min, ang_max,
			ang_nbin, ang_min, ang_max);

	sprintf(plotstr, "%s.Th*180./TMath::Pi():%s.Th*180./TMath::Pi()>>%s", name1,
			name2, hname);
	T->Draw(plotstr, "", "colz");

	/*
	 * Momentum plot
	 */
	c1->cd(j);
	gPad->SetMargin(0.13,0.13,0.13,0.07);
	sprintf(hname, "mom_%s_vs_%s", name1, name2);
	sprintf(htitle,
			"Momentum correlation;Momentum of %s [GeV];Momentum of %s [GeV]",
			name2, name1);
	TH2I* phase_histo = new TH2I(hname, htitle, mom_nbin, mom_min, mom_max,
			mom_nbin, mom_min, mom_max);

	sprintf(plotstr, "%s.P:%s.P>>%s", name1, name2, hname);
	T->Draw(plotstr, "", "colz");
}

TCanvas *drawplot(TTree *T, char *name, char *title, char *cname, char*ctitle,
		bool charge = 1, int width = 1500, int height = 400) {
	TCanvas *tmpc = new TCanvas(cname, ctitle, width, height);
	tmpc->Divide(3, 1);
	makeplot(T, name, title, tmpc, 1, 2, 3, charge);
	tmpc->Update();
	return tmpc;
}

TCanvas *drawplot2(TTree *T, char *name1, char *title1, char *name2,
		char *title2, char *cname, char*ctitle, bool charge = 1, int width =
				1500, int height = 800) {
	TCanvas *tmpc = new TCanvas(cname, ctitle, width, height);
	tmpc->Divide(3, 2);
	makeplot(T, name1, title1, tmpc, 1, 2, 3, charge);
	makeplot(T, name2, title2, tmpc, 4, 5, 6, charge);
	tmpc->Update();
	return tmpc;
}

TCanvas *drawcorr(TTree *T, char *name, char *title, char *cname, char*ctitle,
		int width = 1000, int height = 400) {
	TCanvas *tmpc = new TCanvas(cname, ctitle, width, height);
	tmpc->Divide(2, 1);
	corrplot(T, name, title, tmpc, 1, 2);
	tmpc->Update();
	return tmpc;
}

TCanvas *drawcorr2(TTree *T, char *name1, char *title1, char *name2,
		char *title2, char *cname, char*ctitle, int width = 1000, int height =
				800) {
	TCanvas *tmpc = new TCanvas(cname, ctitle, width, height);
	tmpc->Divide(2, 2);
	corrplot(T, name1, title1, tmpc, 1, 2);
	corrplot(T, name2, title2, tmpc, 3, 4);
	tmpc->Update();
	return tmpc;
}

void mcplot() {
	gStyle->SetOptStat(0);
	double r2d = 180. / TMath::Pi();
	TFile *f1 = new TFile("hd_root.root", "r");
	TDirectory *dir = gROOT->FindObject("DMC");
	dir->cd();
	TTree *T = (TTree*) gROOT->FindObject("thrown");

	/*
	 * Set Alias
	 */
//	T->SetAlias("pi0.M",
//			"sqrt(pow(photon[0].p.fE+photon[1].p.fE,2)-pow(photon[0].p.fP.fX+photon[1].p.fP.fX,2)-pow(photon[0].p.fP.fY+photon[1].p.fP.fY,2)-pow(photon[0].p.fP.fZ+photon[1].p.fP.fZ,2))");
	/*
	 * Make plots
	 */
	TCanvas *c1 = drawplot2(T, "Kp", "K plus", "Km", "K minus", "c1",
			"Kaon distribution");
	c1->Print("./plot/kaon.png");

	TCanvas *c2 = drawplot2(T, "pip", "pi plus", "pim", "pi minus", "c2",
			"Pion distribution");
	c2->Print("./plot/pion.png");

	TCanvas *c3 = drawplot(T, "proton", "proton", "c3", "Proton distribution");
	c3->Print("./plot/proton.png");

	TCanvas *c4 = drawplot2(T, "electron", "electron", "positron", "positron",
			"c4", "Lepton distribution");
	c4->Print("./plot/lepton.png");

	TCanvas *c5 = drawplot(T, "photon", "photon", "c5", "Photon distribution",
			0);
	c5->Print("./plot/photon.png");

	TCanvas *c6 = drawcorr(T, "Kp", "Km", "c6", "Kaon Correlation");
	c6->Print("./plot/kaon_corr.png");

	TCanvas *c7 = drawcorr(T, "pip", "pim", "c7", "Pion Correlation");
	c7->Print("./plot/pion_corr.png");

	TCanvas *c8 = drawcorr(T, "positron", "electron", "c8",
			"Lepton Correlation");
	c8->Print("./plot/lepton_corr.png");
}