//
//  This was copied originally from revision 18550 of pion_polarization_experiments.C
// For this, the projected point from the JLab measurement is added and some
// experiments with larger error bars are dropped so that the y-axis can be zoomed in
// to show better detail around the 5.7 range we're interested in.
//

// This produces a plot summarizing the experimental and theoretical 
// landscsape of alpha_pi - beta_pi.
//
// For some reason, cint crashes on this so it needs to be compiled
// using the external C++ compiler like this:
//
//   .x pion_polarization_experiments.C+
//
// Note that a similar summary plot can be found in a talk by 
// Jan Friedrich titled "Hadron Spectroscopy & Primakoff Reactions at COMPASS"
// shown at the International Workshop on Hadron Structure and Spectroscopy
// in Paris on April 5, 20122 on slide 40.
//=============================================================================
//
// References:
//
// PLUTO @ DESY       (from Morgan review article)
//                    PLUTO Collab., C. Berger et al., Z. Phys C 26 (1984) 199
//
// DM1 @ DCI          (from Morgan review article)
//                    DM1 Collab., A. Courau et al., Nucl. Phys. B 271 (1986) 1
//
// DM2 @ DCI          (from Morgan review article)
//                    DM2 Collb., X. Ajaltouni et al., contrib. VIIth Intern.
//                      Workshop on Photon-photon collisions (Pairs, April 1986)
//                      cited by J.H. Field, DESY report 86-117 (September 1986)
//
// MARK II @ SLAC     (from Babusci review article)
//                    MARK-II Collab., J. Boyer et al., Phys Rev D 42 (1990) 1350
//
// MAMI @ MAINZ
//                    J. Ahrens et al, Eur. Phys. J. A 23, 113–127 (2005)
//
// PACHRA @ LEBEDEV   (from Babusci review article)
//                    T A Aibergenov et al, Czech J Phys 36 (1986) 948
//
// SIGMA @ Serpukov   (from Babusci review article)
//                    Yu M. Antipov et al, Phys Lett B V121 (1983) 445
//                    Yu M. Antipov et al, Z Phys C 26 (1985) 495
//
// COMPASS @ CERN     (from 2015 Phys. Rev. Lett.)
//                    C. Adolph et al, Phys Rev Lett 114 (2015) 062002
//
//
// Review articles:
// ----------------
// D. Morgan, M R Pennington Phys Lett. B V192 N1,2 (1987) 207-211
//
// D. Babusci et al., Phys. Lett. B V277 (1992) 158-162
//


#include <vector>
#include <iostream>
using namespace std;

#include <TH2.h>
#include <TGraphErrors.h>
#include <TLatex.h>
#include <TLine.h>
#include <TCanvas.h>
#include <TColor.h>
#include <TStyle.h>
#include <TArrow.h>

class PiPol{
	public:
		PiPol(string my_name, string my_subtitle, double my_alpha_minus_beta, double my_err_stat, double my_err_sys):name(my_name),subtitle(my_subtitle),alpha_minus_beta(my_alpha_minus_beta),err_stat(my_err_stat),err_sys(my_err_sys){}
	
		string name;
		string subtitle;
		double alpha_minus_beta;
		double err_stat;
		double err_sys;
	
};


TGraphErrors* MakeTGraph(vector<PiPol> &vals, int offset, int color, int style);
void DrawLabels(vector<PiPol> &vals, int offset);
void DrawBracket(TH2D *axes, int first_bin, int Nbins, int color, string slabel, string sublabel="");


void pion_polarization_experiments_JLab_projected(void)
{
	TColor::CreateColorWheel();
	gStyle->SetLineWidth(2);
	gStyle->SetPadTopMargin(0.05);
	gStyle->SetPadBottomMargin(0.15);
	
	vector<PiPol> experiments1;
	vector<PiPol> experiments2;
	vector<PiPol> experiments3;
	vector<PiPol> experiments4;
	vector<PiPol> theory;
	
	// gamma* gamma* -> pi+ pi-
//	experiments1.push_back(PiPol("PLUTO", "DESY", 2.0*19.1, 2.0*4.8, 2.0*5.7));
//	experiments1.push_back(PiPol("DM1", "DCI", 2.0*17.2, 2.0*4.6, 0.0));
//	experiments1.push_back(PiPol("DM2", "DCI", 2.0*26.3, 2.0*7.4, 0.0));
	experiments1.push_back(PiPol("MARK II", "SLAC", 2.0*2.2, 2.0*1.6, 0.0));
	
	// gamma pi+* -> n pi+ gamma
	experiments2.push_back(PiPol("MAMI", "MAINZ", 11.6, 1.5, sqrt(3.0*3.0+0.5+0.5)));
//	experiments2.push_back(PiPol("PACHRA", "LEBEDEV", 2.0*20.0, 2.0*12.0, 0.0));
	
	// pi- A -> A pi- gamma
	experiments3.push_back(PiPol("SIGMA", "Serpukov", 13.6, 2.8, 2.4));
	experiments3.push_back(PiPol("COMPASS", "CERN", 2.0*2.0, 2.0*0.6, 2.0*0.7));

	// gamma gamma* -> pi+ pi-
	// JLab CPP has 0.6% stat, 1.5% sys on cross-section but that projects
	// to an estimated 10% on alpha-beta. Break the 10% into stat and sys
	// parts that have the same ratios
	double stat = 0.10 * 0.006 / sqrt(pow(0.006,2.0) + pow(0.015, 2.0));
	double sys  = 0.10 * 0.015 / sqrt(pow(0.006,2.0) + pow(0.015, 2.0));
	cout << "total JLab CPP error: " << sqrt(stat*stat + sys*sys) << endl;
	experiments4.push_back(PiPol("GlueX(projected)", "JLab", 5.7, stat*5.7, sys*5.7));

	// theory
	theory.push_back(PiPol("Gasser", "ChPT", 5.7, 1.0, 0.0));
	theory.push_back(PiPol("Fil'kov", "DR", 13.6, 2.6, 1.9));
	theory.push_back(PiPol("Pasquini", "DR", 5.7, 0.0, 0.0));

	unsigned int Nvals = 0;

	TGraphErrors *gexperiments1 = MakeTGraph(experiments1, Nvals, kBlack, 23);
	Nvals += experiments1.size()+1;

	TGraphErrors *gexperiments2 = MakeTGraph(experiments2, Nvals, kBlue, 20);
	Nvals += experiments2.size()+1;
	
	TGraphErrors *gexperiments3 = MakeTGraph(experiments3, Nvals, kMagenta, 22);
	Nvals += experiments3.size()+1;
	
	TGraphErrors *gexperiments4 = MakeTGraph(experiments4, Nvals, kRed, 8);
	Nvals += experiments4.size()+1;
	
	Nvals+=3; // for space between theory and experiment
	TGraphErrors *gtheory = MakeTGraph(theory, Nvals, kBlack, 21);
	Nvals += theory.size()+1;
	
	double ymax = 18.0;
	TH2D *axes = new TH2D("axes","", Nvals, 0.0, (double)Nvals+0.5, 100, 0.0, ymax);
	axes->SetStats(0);
	axes->SetYTitle("#alpha_{#pi} - #beta_{#pi} (#times10^{-4} fm^{3})");

	TCanvas *c1 = new TCanvas("c1","", 800, 600);
	c1->SetTicks();
	
	// Setting bin label to empty string prevents numeric values
	// from being drawn on X-axis
	axes->GetXaxis()->SetBinLabel(1,"");
	
	axes->Draw();

	// Draw dashed horizontal line at Gasser's value
	double x_ChPT = (double)Nvals + 1.0;
	TLine *lin = new TLine(0.0, 5.7, x_ChPT - 0.5, 5.7);
	lin->SetLineStyle(1);
	lin->SetLineWidth(2);
	lin->SetLineColor(kGray+2);
	lin->Draw();

	lin = new TLine(0.0, 5.7+1.0, x_ChPT + 0.5, 5.7+1.0);
	lin->SetLineStyle(2);
	lin->SetLineWidth(1);
	lin->SetLineColor(kGray+2);
	lin->Draw();
	
	lin = new TLine(0.0, 5.7-1.0, x_ChPT + 0.5, 5.7-1.0);
	lin->SetLineStyle(2);
	lin->SetLineWidth(1);
	lin->SetLineColor(kGray+2);
	lin->Draw();
	
	// Label ChPT prediction on left hand side
	TLatex *lab_ChPT = new TLatex(x_ChPT-0.4, 5.7, "ChPT");
	lab_ChPT->SetTextSize(0.030);
	lab_ChPT->SetTextAlign(12);
	lab_ChPT->Draw();
	
	double xarrows = x_ChPT + 0.1;
	TArrow *arrow_ChPT_top = new TArrow(xarrows, 5.7+1.0+1.0, xarrows, 5.7+1.0+0.0, 0.012, "-|>");
	arrow_ChPT_top->Draw();
	TArrow *arrow_ChPT_bot = new TArrow(xarrows, 5.7-1.0-1.0, xarrows, 5.7-1.0-0.0, 0.012, "-|>");
	arrow_ChPT_bot->Draw();

	
	// Draw vertical line to separate experiment and theory
	double x_separator = (double)(Nvals-theory.size()) - 1.25;
	lin = new TLine(x_separator, 0.0, x_separator, axes->GetYaxis()->GetXmax());
	lin->SetLineStyle(7);
	lin->SetLineWidth(3);
	lin->SetLineColor(kGray+1);
	lin->Draw();
	
	TLatex *lab_theory = new TLatex(x_separator+0.1, 0.85*ymax, "theory");
	lab_theory->SetTextAngle(270.0);
	lab_theory->SetTextSize(0.03);
	lab_theory->SetTextAlign(21);
	lab_theory->Draw();

	TLatex *lab_experiment = new TLatex(x_separator-0.1, 0.85*ymax, "experiment");
	lab_experiment->SetTextAngle(90.0);
	lab_experiment->SetTextSize(0.03);
	lab_experiment->SetTextAlign(21);
	lab_experiment->Draw();
	
	TArrow *arrow_theory = new TArrow(x_separator+0.40, 0.85*ymax, x_separator+0.75, 0.85*ymax, 0.012, "-|>");
	arrow_theory->SetLineWidth(5);
	arrow_theory->Draw();
	TArrow *arrow_experiment = new TArrow(x_separator-0.40, 0.85*ymax, x_separator-0.75, 0.85*ymax, 0.012, "-|>");
	arrow_experiment->SetLineWidth(5);
	arrow_experiment->Draw();

	// Draw data points w/ error bars
	gexperiments1->Draw("Psame");
	gexperiments2->Draw("Psame");
	gexperiments3->Draw("Psame");
	gexperiments4->Draw("Psame");
	gtheory->Draw("Psame");
	
	// Draw labels next to data points and reactions on bottom
	Nvals = 0;

	DrawLabels(experiments1, Nvals);
	DrawBracket(axes, Nvals+1, experiments1.size(), kBlack, "#gamma#gamma#rightarrow#pi^{+}#pi^{-}", "Collider");
	Nvals += experiments1.size() + 1;

	DrawLabels(experiments2, Nvals);
	DrawBracket(axes, Nvals+1, experiments2.size(), kBlue, "#gammap#rightarrown#pi^{+}#gamma");
	Nvals += experiments2.size() + 1;

	DrawLabels(experiments3, Nvals);
	DrawBracket(axes, Nvals+1, experiments3.size(), kMagenta, "#piA#rightarrow#pi'#gammaA", "Primakoff");
	Nvals += experiments3.size() + 1;

	DrawLabels(experiments4, Nvals);
	DrawBracket(axes, Nvals+1, experiments4.size()+2, kRed, "#gammaA#rightarrow#pi^{+}#pi^{-}A", "Primakoff");
	Nvals += experiments4.size() + 1;

	Nvals+=3;
	DrawLabels(theory, Nvals);
	DrawBracket(axes, Nvals+1, theory.size(), kBlack, "Theory Predictions");
	Nvals += theory.size() + 1;

	// Special notes at bottom
//	TLatex *note1 = new TLatex(15.0, -6.0, "#it{*COMPASS result from DNP2013 abstract NJ00005}");
//	note1->SetTextSize(0.020);
//	note1->SetTextAlign(32);
//	note1->Draw();
	
	// Save output to file
	c1->SaveAs("pion_polarization_experiments_JLab_projected.pdf");
	c1->SaveAs("pion_polarization_experiments_JLab_projected.png");
}

//-------------------
// MakeTGraph
//-------------------
TGraphErrors* MakeTGraph(vector<PiPol> &vals, int offset, int color, int style)
{
	double x[100];
	double y[100];
	double ey[100];
	for(unsigned int i=0; i<vals.size(); i++){
		x[i] = (double)(i+offset)+0.5;
		y[i] = vals[i].alpha_minus_beta;
		double stat = vals[i].err_stat;
		double sys = vals[i].err_sys;
		ey[i] = sqrt(stat*stat + sys*sys);
		
	}

	TGraphErrors *graph = new TGraphErrors(vals.size(), x, y, 0, ey);
	graph->SetMarkerStyle(style);
	graph->SetMarkerColor(color);
	graph->SetLineColor(color);

	return graph;
}

//-------------------
// DrawLabels
//-------------------
void DrawLabels(vector<PiPol> &vals, int offset)
{
	for(unsigned int i=0; i<vals.size(); i++){
		double x = (double)(i+offset)+0.5+0.085;
		double y = vals[i].alpha_minus_beta+0.1;
		
		bool has_subtitle = vals[i].subtitle.size()>0;
		
		TLatex *lab = new TLatex(x+0.02, y, vals[i].name.c_str());
		lab->SetTextSize(0.03);
		lab->SetTextAlign(has_subtitle ? 11:12);
		lab->Draw();

		if(has_subtitle){
			string subtitle = string("(")+vals[i].subtitle+")";
			lab = new TLatex(x+0.02, y-0.2, subtitle.c_str());
			lab->SetTextSize(0.025);
			lab->SetTextAlign(13);
			lab->Draw();
		}
	}
}

//-------------------
// DrawBracket
//-------------------
void DrawBracket(TH2D *axes, int first_bin, int Nbins, int color, string slabel, string sublabel)
{
	double binwidth = axes->GetXaxis()->GetBinWidth(1);
	double x1 = axes->GetXaxis()->GetBinLowEdge(first_bin);
	double x2 = axes->GetXaxis()->GetBinLowEdge(first_bin+Nbins) + 0.25*binwidth;
	
	static int Nbrackets = 0;
	double yshift = 0.65*((double)(Nbrackets%2)-0.5);
	double y = -0.75 + yshift;
	TLine *lin = new TLine(x1, y, x2, y);
	lin->SetLineWidth(3);
	lin->SetLineColor(color);
	lin->Draw();

	lin = (TLine*)lin->Clone();
	lin->SetX1(x1);  lin->SetX2(x1);
	lin->SetY1(y);  lin->SetY2(0.75*y);
	lin->Draw();

	lin = (TLine*)lin->Clone();
	lin->SetX1(x2);  lin->SetX2(x2);
	lin->SetY1(y);  lin->SetY2(0.75*y);
	lin->Draw();
	
	double x = (x1+x2)/2.0;
	TLatex *lab = new TLatex(x, 1.25*y, slabel.c_str());
	lab->SetTextAlign(23);
	lab->SetTextSize(0.035);
	lab->Draw();
	
	if(sublabel!=""){
		TLatex *lab2 = new TLatex();
		lab2->SetTextAlign(23);
		lab2->SetTextSize(0.030);
		lab2->DrawLatex(x, y-0.95, sublabel.c_str());
	}

	Nbrackets++;
}