# read histograms and save as yaml file
#
#
import ROOT
from hepdata_lib import RootFileReader

# See https://hepdata-lib.readthedocs.io for how to use the "hepdata_lib" tool.
from hepdata_lib import Submission, Table, Variable, Uncertainty
import numpy as np

#######################################################################################
# Define the main Submission object and add a comment and abstract
submission = Submission()

#submission.read_abstract("abstract.txt")
submission.add_link("arXiv", "https://arxiv.org/abs/2208.02740")
submission.add_record_id(2132332, "inspire")
submission.add_additional_resource("Original abstract file", "abstract.txt", copy_file=True)
submission.add_additional_resource("Original submission python-file", "makeSubmission.py", copy_file=True)
submission.add_related_recid(102468)

#######################################################################################
# root file:
# HADES_Flow_AuAu123_vn_pdt_stat_sys_EPJA59-80.root
# root-file containes for each harmonic v1, v2, v3, v4, # each particle p, d, t
# and for the four centrality classes two 2D-histogramms (one with stat unc. 
# and the other with syst. unc. in the error bars)
# additional for v4 two pt-bin-widths (200 and 100 MeV/c) are provided

descriptionRoot = r"Original root-file (full data-set) contains for each harmonic v1, v2, v3, v4, each particle p, d, t and for the four 10% centrality classes two 2D-histograms (one with statistical and the other with systematic uncertainties as error values) additional for v4 two pt-bin-widths (200 and 100 MeV/c) are provided."
filenameRoot  = "data_root/HADES_Flow_AuAu123_vn_pdt_stat_sys_EPJA59-80.root"
readerRoot    = RootFileReader(filenameRoot)
submission.add_additional_resource(descriptionRoot, filenameRoot, copy_file=True)

#######################################################################################

PARTICLE  = ["Protons", "Deuterons", "Tritons"]
PARTICLES = ["p", "d", "t"]
PARTREAC  = ["AU AU --> P X" , "AU AU --> DEUT X" , "AU AU --> TRIT X"]
PARTLOC   = ["left", "middle", "right"]

HARMONI  = ["v1", "v2", "v3"]
HARMONIC = ["v1", "v2", "v3", "v4"]
HARMLOC  = ["top", "second", "third", "bottom"]

PTRANGE   = [(150, 2000), (200, 2000), (300, 2000)]

RAPYRANGE = [(-0.65, 0.85), (-0.65, 0.65), (-0.65, 0.35)]

MOMENTBINS = [r"600-900", "1500-1800"]
MOMENTTEXT = [r"0.6 < $p_{t}$ < 0.9", "1.5 < $p_{t}$ < 1.8"]

HARMVAR = [r"d$v_{1}$/d$y^{\prime}|_{y^{\prime} = 0}$", r"$v_{2}$", r"d$v_{3}$/d$y^{\prime}|_{y^{\prime} = 0}$", r"$v_{4}$"]
VNPANEL =["upper left","upper right","lower left","lower right"]

YBINUP    = [0.05, -0.15, 0.25, -0.35, 0.45, -0.55, 0.65]
ybinsize  = 0.1

PTBINUP   = [0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]
ptbinsize = 0.049

PTBINUPV4   = [0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]
ptbinsizeV4 = 0.199


#######################################################################################

table_fig5_10 = Table("Figure 5 Event Plane Resolution 10pct")
table_fig5_10.description = r"The resolution $\Re_{n}$ of the first-order spectator event plane for flow coefficients of different orders $n$ as a function of the event centrality (Adamczewski-Musch:2020iio). The circles correspond to centrality intervals of $5 \%$ width and the squares to $10 \%$ width (curves are meant to guide the eye)."
table_fig5_10.keywords["observables"] = ["EP"]
table_fig5_10.keywords["phrases"]     = ["Flow", "Event Plane Resolution"]
table_fig5_10.location = "Data from Figure 5 (square marker), located on page 6."
table_fig5_10.add_image("figure/EPres_vsCentAuAu1230_0420_ar2.png")

### load data
data = np.loadtxt("data/EPres_vsCent10_AuAu1230_0420.dat", skiprows=2)

### EP1
# x-axis: Centrality
EP10 = Variable("Centrality", is_independent=True, is_binned=False, units="%")
EP10.values = data[:,0]
table_fig5_10.add_variable(EP10)

# y-axis
for x in range(1,7):
    text = r"$R_{" + str(x) + r"}$"
    EP10_R = Variable(text, is_independent=False, is_binned=False, units="")
    EP10_R.values = data[:,x]
    table_fig5_10.add_variable(EP10_R)

table_fig5_10.add_additional_resource("Original data file", "data/EPres_vsCent10_AuAu1230_0420.dat", copy_file=True)
submission.add_table(table_fig5_10)

#######################################################################################

table_fig5_5 = Table("Figure 5 Event Plane Resolution 5pct")
table_fig5_5.description = r"The resolution $\Re_{n}$ of the first-order spectator event plane for flow coefficients of different orders $n$ as a function of the event centrality (Adamczewski-Musch:2020iio). The circles correspond to centrality intervals of $5 \%$ width and the squares to $10 \%$ width (curves are meant to guide the eye)."
table_fig5_5.keywords["observables"] = ["EP"]
table_fig5_5.keywords["phrases"]     = ["Flow", "Event Plane Resolution"]
table_fig5_5.location = "Data from Figure 5 (circle marker), located on page 6."
table_fig5_5.add_image("figure/EPres_vsCentAuAu1230_0420_ar2.png")

### load data
data = np.loadtxt("data/EPres_vsCent5_AuAu1230_0420.dat", skiprows=2)

### EP1
# x-axis: Centrality
EP5 = Variable("Centrality", is_independent=True, is_binned=False, units="%")
EP5.values = data[:,0]
table_fig5_5.add_variable(EP5)

# y-axis
for x in range(1,7):
    text = r"$R_{" + str(x) + r"}$"
    EP5_R = Variable(text, is_independent=False, is_binned=False, units="")
    EP5_R.values = data[:,x]
    table_fig5_5.add_variable(EP5_R)

table_fig5_5.add_additional_resource("Original data file", "data/EPres_vsCent5_AuAu1230_0420.dat", copy_file=True)
submission.add_table(table_fig5_5)

#######################################################################################
## Fig. 8  (4x3x8)  y_cm (pt-slice)

description = r"The flow coefficients $v_{1}$, $v_{2}$, $v_{3}$, and $v_{4}$ (from top to bottom panels) of protons, deuterons and tritons (from left to right panels) in semi-central ($20 - 30 \%$) Au+Au collisions at $\sqrt{s_{NN}} = 2.4$ GeV as a function of the centre-of-mass rapidity $y_{cm}$ in transverse momentum intervals of $50$ MeV$/c$ width. Systematic uncertainties are displayed as boxes. Lines are to guide the eye."

for iharm, harmonic in enumerate(HARMONI):
    
    for ipart, particle in enumerate(PARTICLE):
        
        # Define a table by the independent and dependent variables.
        table = Table(f"Figure 8 {HARMONIC[iharm]} ycm {particle}")

        # Add some metadata to the table.
        table.description = description
        table.location = f"Data from Figure 8 ({PARTLOC[ipart]} column {HARMLOC[iharm]} row), located on page 10."

        # Add a thumbnail image of the original figure from the paper.
        table.add_image(f"figure/vn_frag_2030Cent_Y_3x1_n{iharm+1}.pdf")
        
        table.keywords["observables"] = [HARMONIC[iharm].upper()]
        table.keywords["reactions"]   = ["AU AU" , PARTREAC[ipart] ]
        table.keywords["cmenergies"]  = ["2.4"]
        table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy","Inclusive", "Rapidity Dependence"]
        
        ### x-axis  (y_cm)
        histnameData     = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_data"
        histnameSys      = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_sys"
        
        histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=RAPYRANGE[ipart] , ylim=(500,550))

        # Create variable objects
        x = Variable(r"$y_{cm}$", is_independent=True, is_binned=True)
        x.values = histPtY_data["x_edges"]
        #print(histPtY_data["x_edges"])
        #print(x.values)
        table.add_variable(x)

        for ipt, ptbins in enumerate(PTBINUP):

            # reload pt-range
            ptmin = 1000.*(ptbins-ptbinsize)
            ptmax = 1000.*(ptbins)
            histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=RAPYRANGE[ipart] , ylim=(ptmin,ptmax))
            histPtY_sys  = readerRoot.read_hist_2d(histnameSys  , xlim=RAPYRANGE[ipart] , ylim=(ptmin,ptmax))
            
            # y-axis
            vn = Variable(f"$v_{iharm+1}$", is_independent=False, is_binned=False, zero_uncertainties_warning=False)
            vn.values = histPtY_data["z"]
            #print(vn.values)
        
            vnStat = Uncertainty("stat", is_symmetric=True)
            vnStat.values = histPtY_data["dz"]
            vn.add_uncertainty(vnStat)
            #print(vnStat.values)
        
            vnSys = Uncertainty("syst", is_symmetric=True)
            vnSys.values = histPtY_sys["dz"]
            vn.add_uncertainty(vnSys)
            
            vn.add_qualifier("$p_{t}$ [GeV/c]", f"{round(ptmin/1000.,2)} - {round(ptmax/1000.,2)}")
            vn.add_qualifier("Centrality", "20 - 30 %")
                            
            #table.add_variable(y)
            table.add_variable(vn)
        
        # Add the table to the Submission object.
        submission.add_table(table)
#######################################################################################
## Fig. 8  (4x3x8)  y_cm (pt-slice)  for v4
description = r"The flow coefficients $v_{1}$, $v_{2}$, $v_{3}$, and $v_{4}$ (from top to bottom panels) of protons, deuterons and tritons (from left to right panels) in semi-central ($20 - 30 \%$) Au+Au collisions at $\sqrt{s_{NN}} = 2.4$ GeV as a function of the centre-of-mass rapidity $y_{cm}$ in transverse momentum intervals of $200$ MeV$/c$ width. Systematic uncertainties are displayed as boxes. Lines are to guide the eye."

iharm = 3
for ipart, particle in enumerate(PARTICLE):
    
    # Define a table by the independent and dependent variables.
    table = Table(f"Figure 8 {HARMONIC[iharm]} ycm {particle}")

    # Add some metadata to the table.
    table.description = description
    table.location = f"Data from Figure 8 ({PARTLOC[ipart]} column {HARMLOC[iharm]} row), located on page 10."

    # Add a thumbnail image of the original figure from the paper.
    table.add_image(f"figure/vn_frag_2030Cent_Y_3x1_n{iharm+1}.pdf")
    
    table.keywords["observables"] = [HARMONIC[iharm].upper()]
    table.keywords["reactions"]   = ["AU AU" , PARTREAC[ipart] ]
    table.keywords["cmenergies"]  = ["2.4"]
    table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy","Inclusive", "Rapidity Dependence"]
    
    ### x-axis  (y_cm)
    histnameData     = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_data"
    histnameSys      = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_sys"
    histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=RAPYRANGE[ipart] , ylim=(600,800))

    # Create variable objects
    x = Variable(r"$y_{cm}$", is_independent=True, is_binned=True)
    x.values = histPtY_data["x_edges"]
    #print(histPtY_data["x_edges"])
    #print(x.values)
    table.add_variable(x)

    for ipt, ptbins in enumerate(PTBINUPV4):

        # reload pt-range
        ptmin = 1000.*(ptbins-ptbinsizeV4)
        ptmax = 1000.*(ptbins)
        
        histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=RAPYRANGE[ipart] , ylim=(ptmin,ptmax))
        histPtY_sys  = readerRoot.read_hist_2d(histnameSys , xlim=RAPYRANGE[ipart] ,  ylim=(ptmin,ptmax))
        
        # y-axis
        vn = Variable(f"$v_{iharm+1}$", is_independent=False, is_binned=False, zero_uncertainties_warning=False)
        vn.values = histPtY_data["z"]
    
        vnStat = Uncertainty("stat", is_symmetric=True)
        vnStat.values = histPtY_data["dz"]
        vn.add_uncertainty(vnStat)
    
        vnSys = Uncertainty("syst", is_symmetric=True)
        vnSys.values = histPtY_sys["dz"]
        vn.add_uncertainty(vnSys)
        
        vn.add_qualifier("$p_{t}$ [GeV/c]", f"{round(ptmin/1000.,2)} - {round(ptmax/1000.,2)}")
        vn.add_qualifier("Centrality", "20 - 30 %")
    
        #table.add_variable(y)
        table.add_variable(vn)
    
    # Add the table to the Submission object.
    submission.add_table(table)

#######################################################################################
## Fig. 9     (4x3x7)  pt   (y_cm-slice)

description = r"The flow coefficients $v_{1}$, $v_{2}$, $v_{3}$, and $v_{4}$ (from top to bottom panels) of protons, deuterons and tritons (from left to right panels) in semi-central ($20 - 30 \%$) Au+Au  collisions at $\sqrt{s_{NN}} = 2.4$ GeV as a function of $p_{t}$ in several rapidity intervals chosen symmetrically around mid-rapidity. Systematic uncertainties are displayed as boxes."
  
for iharm, harmonic in enumerate(HARMONI):
    
    for ipart, particle in enumerate(PARTICLE):
        
        # Define a table by the independent and dependent variables.
        table = Table(f"Figure 9 {HARMONIC[iharm]} pt {particle}")

        # Add some metadata to the table.
        table.description = description
        table.location = f"Data from Figure 9 ({PARTLOC[ipart]} column {HARMLOC[iharm]} row), located on page 11."

        # Add a thumbnail image of the original figure from the paper.
        table.add_image(f"figure/vn_frag_2030Cent_Pt_3x1_n{iharm+1}.pdf")
        
        table.keywords["observables"] = [HARMONIC[iharm].upper()]
        table.keywords["reactions"]   = ["AU AU" , PARTREAC[ipart] ]
        table.keywords["cmenergies"]  = ["2.4"]
        table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy","Inclusive", "Transverse Momentum Dependence"]
        
        ### x-axis  (p_t)
        histnameData     = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_data"
        histnameSys      = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_sys"

        histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=(-0.05,0.05) , ylim=PTRANGE[ipart])

        # Create variable objects
        x = Variable(r"$p_{t}$", is_independent=True, is_binned=True)
        myArray = np.array(histPtY_data["y_edges"])/1000.
        #print(myArray)
        x.values = myArray
        table.add_variable(x)

        for iy, ybins in enumerate(YBINUP):

            # reload pt-range
            ymin = (ybins-ybinsize)
            ymax = (ybins)
            histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=(ymin,ymax) , ylim=PTRANGE[ipart])
            histPtY_sys  = readerRoot.read_hist_2d(histnameSys , xlim=(ymin,ymax) ,  ylim=PTRANGE[ipart])
            
            # y-axis
            vn = Variable(f"$v_{iharm+1}$", is_independent=False, is_binned=False, zero_uncertainties_warning=False)
            vn.values = histPtY_data["z"]
        
            vnStat = Uncertainty("stat", is_symmetric=True)
            vnStat.values = histPtY_data["dz"]
            vn.add_uncertainty(vnStat)
        
            vnSys = Uncertainty("syst", is_symmetric=True)
            vnSys.values = histPtY_sys["dz"]
            vn.add_uncertainty(vnSys)
            
            vn.add_qualifier("$y_{cm}$", f"{round(ymin,2)} - {round(ymax,2)}")
            vn.add_qualifier("Centrality", "20 - 30 %")

            table.add_variable(vn)
        
        # Add the table to the Submission object.
        submission.add_table(table)

#######################################################################################
## Fig. 9     (4x3x7)  pt   (y_cm-slice) for v4
iharm = 3
for ipart, particle in enumerate(PARTICLE):
    
    # Define a table by the independent and dependent variables.
    table = Table(f"Figure 9 {HARMONIC[iharm]} pt {particle}")

    # Add some metadata to the table.
    table.description = description
    table.location = f"Data from Figure 9 ({PARTLOC[ipart]} column {HARMLOC[iharm]} row), located on page 11."

    # Add a thumbnail image of the original figure from the paper.
    table.add_image(f"figure/vn_frag_2030Cent_Pt_3x1_n{iharm+1}.pdf")
    
    table.keywords["observables"] = [HARMONIC[iharm].upper()]
    table.keywords["reactions"]   = ["AU AU" , PARTREAC[ipart] ]
    table.keywords["cmenergies"]  = ["2.4"]
    table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy","Inclusive", "Transverse Momentum Dependence"]
    
    ### x-axis  (y_cm)
    histnameData     = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_data_pt100MeVc"
    histnameSys     = f"PtY_v{iharm+1}EPcent2030_COSv{iharm+1}_{PARTICLES[ipart]}_sys_pt100MeVc"
    histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=(-0.05,0.05) , ylim=PTRANGE[ipart])

    # Create variable objects
    x = Variable(r"$p_{cm}$", is_independent=True, is_binned=True)
    myArray = np.array(histPtY_data["y_edges"])/1000.
    #print(myArray)
    x.values = myArray
    
    table.add_variable(x)

    for iy, ybins in enumerate(YBINUP):

        # reload pt-range
        ymin = (ybins-ybinsize)
        ymax = (ybins)
        histPtY_data = readerRoot.read_hist_2d(histnameData , xlim=(ymin,ymax) , ylim=PTRANGE[ipart])
        histPtY_sys  = readerRoot.read_hist_2d(histnameSys  , xlim=(ymin,ymax) , ylim=PTRANGE[ipart])
        
        # y-axis
        vn = Variable(f"$v_{iharm+1}$", is_independent=False, is_binned=False, zero_uncertainties_warning=False)
        vn.values = histPtY_data["z"]
    
        vnStat = Uncertainty("stat", is_symmetric=True)
        vnStat.values = histPtY_data["dz"]
        vn.add_uncertainty(vnStat)
    
        vnSys = Uncertainty("syst", is_symmetric=True)
        vnSys.values = histPtY_sys["dz"]
        vn.add_uncertainty(vnSys)
        
        vn.add_qualifier("$y_{cm}$", f"{round(ymin,2)} - {round(ymax,2)}")
        vn.add_qualifier("Centrality", "20 - 30 %")

        table.add_variable(vn)
    
    # Add the table to the Submission object.
    submission.add_table(table)


#######################################################################################
## Fig. 10    (4x6)    int. in two pt-bin v1-v4

description = r"Directed (d$v_{1}$/d$y^{\prime}|_{y^{\prime} = 0}$, upper left panel), triangular (d$v_{3}$/d$y^{\prime}|_{y^{\prime} = 0}$, upper right panel), elliptic ($v_{2}$, lower left panel) and quadrangular ($v_{4}$, lower right panel) flow of protons, deuterons and tritons in two transverse momentum intervals (open symbols: $0.6 < p_{t} < 0.9$ GeV$/c$ and filled symbols: $1.5 <  p_{t} < 1.8$ GeV$/c$) at mid-rapidity in Au+Au collisions at $\sqrt{s_{NN}} = 2.4$ GeV for four centrality classes. Systematic uncertainties are displayed as boxes."

centrality = np.array([5 ,15 ,25 ,35])

for iharm, harmonic in enumerate(HARMONIC):
    
    # Define a table by the independent and dependent variables.
    table = Table(f"Figure 10 {harmonic} Centrality")

    # Add some metadata to the table.
    table.description = description
    table.location = f"Data from Figure 10 ({VNPANEL[iharm]} panel), located on page 12."

    # Add a thumbnail image of the original figure from the paper.
    table.add_image(f"figure/v{iharm+1}_midrap_pdt_allCent.pdf")
    
    table.keywords["observables"] = [HARMONIC[iharm].upper()]
    table.keywords["reactions"]   = ["AU AU" , "AU AU --> P X" , "AU AU --> DEUT X" , "AU AU --> TRIT X" ]
    table.keywords["cmenergies"]  = ["2.4"]
    table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy","Inclusive", "Transverse Momentum Dependence" , "Centrality Dependence"]
    
    # x-axis: Centrality
    varVnX = Variable("Centrality", is_independent=True, is_binned=False, units="%")
    varVnX.values = centrality
    table.add_variable(varVnX)
    
    # Data for Particle and pt intervals
    for ipart, particle in enumerate(PARTICLE):
        
        for imom, momentum in enumerate(MOMENTBINS):
            
            ### load data
            desc     = f"Data file: v{iharm+1} pt{momentum} {PARTICLES[ipart]}"
            filename = f"data/HADES_v{iharm+1}_fit_au1230_pt{momentum}_{PARTICLES[ipart]}.dat"
            data = np.loadtxt(filename)

            # y-axis
            varVnY = Variable(HARMVAR[iharm], is_independent=False, is_binned=False, units="")
            if((iharm+1)%2) :
                varVnY.values = data[:,1]*0.74
            else:
                varVnY.values = data[:,1]

            varVnY.add_qualifier("Particle", particle)
            varVnY.add_qualifier("Pt [GeV/c]", MOMENTTEXT[imom])
        
            varVnYSys = Uncertainty("syst", is_symmetric=True)
            if((iharm+1)%2) :
                varVnYSys.values = data[:,2]*0.74
            else:
                varVnYSys.values = data[:,2]

            varVnY.add_uncertainty(varVnYSys)
            table.add_variable(varVnY)
            table.add_additional_resource(desc, filename , copy_file=True)
       
    # Add the table to the Submission object.
    submission.add_table(table)

#######################################################################################
## Fig. 11 v1 v2

description = r"Shown as red points are the slope of $v_{1}$ at mid-rapidity (left panel), d$v_{1}$/d$y^{\prime}|_{y^{\prime} = 0}$, and the $p_{t}$ integrated $v_{2}$ at mid-rapidity (right panel) for protons in Au+Au collisions at $\sqrt{s_{NN}} = 2.4$ GeV ($10 - 30 \%$ centrality)."

energy = np.array([2.4])
varVn  = np.array([2.4, 2.4])
varUnc = np.array([2.4, 2.4])

# Define a table by the independent and dependent variables.
table = Table(f"Figure 11 v1v2 Energy")

# Add some metadata to the table.
table.description = description
table.location = f"Data from Figure 11, located on page 14."

# Add a thumbnail image of the original figure from the paper.
# FIXME table.add_image(f"figure/v{iharm+1}_midrap_pdt_allCent.pdf")

table.keywords["observables"] = ["V1", "V2"]
table.keywords["reactions"]   = ["AU AU" , "AU AU --> P X" ]
table.keywords["cmenergies"]  = ["2.4"]
table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy","Inclusive", "Transverse Momentum Dependence" , "Centrality Dependence"]

# x-axis: Centrality
varVnX = Variable("Energy", is_independent=True, is_binned=False, units="GeV")
varVnX.values = np.array([2.4])
table.add_variable(varVnX)

# y-axis
varV1 = Variable(HARMVAR[0], is_independent=False, is_binned=False, units="")
varV1.values = np.array([0.459])
varV1.add_qualifier("Particle", "Protons")
varV1.add_qualifier("Centrality", "10 - 30 %")

varV1Sys = Uncertainty("syst", is_symmetric=True)
varV1Sys.values = np.array([0.022])
varV1.add_uncertainty(varV1Sys)
table.add_variable(varV1)

# y-axis
varV2 = Variable(HARMVAR[1], is_independent=False, is_binned=False, units="")
varV2.values = np.array([-0.06])
varV2.add_qualifier("Particle", "Protons")
varV2.add_qualifier("Centrality", "10 - 30 %")

varV2Sys = Uncertainty("syst", is_symmetric=True)
varV2Sys.values = np.array([0.001])
varV2.add_uncertainty(varV2Sys)
table.add_variable(varV2)

# Add the table to the Submission object.
submission.add_table(table)


#######################################################################################

table_tab4 = Table("Table 4 Eccentricities")
table_tab4.description = r"Parameters describing the initial nucleon distribution for the different centrality classes as calculated within the Glauber-MC approach (Adamczewski-Musch:2017sdk).  Listed are the corresponding average impact parameter $\langle b \rangle$ and the average participant eccentricities $\langle \epsilon_{2} \rangle$ and $\langle \epsilon_{4} \rangle$."
table_tab4.keywords["observables"] = ["Eccentricities"]
table_tab4.keywords["phrases"]     = ["Flow", "Eccentricities", "Glauber MC"]
table_tab4.location = "Data from Table 4, located on page 16."

### load data
data = np.loadtxt("data/Eccentricities_GlauberMC_vsCent10_AuAu1230.dat", skiprows=2)
#print(data[:,0:2])

# x-axis: Centrality
cent = Variable("Centrality", is_independent=True, is_binned=True, units="%")
cent.values = data[:,0:2]
table_tab4.add_variable(cent)

ecc = Variable(r"$\langle b \rangle$", is_independent=False, is_binned=False, units="fm")
ecc.values = data[:,2]
table_tab4.add_variable(ecc)


ecc = Variable(r"$\langle \varepsilon_{2} \rangle$", is_independent=False, is_binned=False, units="")
ecc.values = data[:,3]
eccUnc = Uncertainty("model uncert.", is_symmetric=True)
eccUnc.values = data[:,4]
ecc.add_uncertainty(eccUnc)
table_tab4.add_variable(ecc)

ecc = Variable(r"$\langle \varepsilon_{4} \rangle$", is_independent=False, is_binned=False, units="")
ecc.values = data[:,5]
eccUnc = Uncertainty("model uncert.", is_symmetric=True)
eccUnc.values = data[:,6]
ecc.add_uncertainty(eccUnc)
table_tab4.add_variable(ecc)

table_tab4.add_additional_resource("Data File: Eccentricities_GlauberMC_vsCent10_AuAu1230", "data/Eccentricities_GlauberMC_vsCent10_AuAu1230.dat", copy_file=True)
submission.add_table(table_tab4)

#######################################################################################
## Fig. 15    

description = r"Scaled elliptic flow of protons, deuterons and tritons in two transverse momentum intervals at mid-rapidity in Au+Au collisions at $\sqrt{s_{NN}} = 2.4$ GeV for four centrality classes. The values are divided by the second order eccentricity, $v_{2} / \langle \epsilon_{2} \rangle$, as calculated within the Glauber-MC approach for the corresponding centrality range (see Table 4 Eccentricities). Systematic uncertainties are displayed as boxes."

centrality = np.array([5 ,15 ,25 ,35])

# Define a table by the independent and dependent variables.
table = Table(f"Figure 15 v2eps2 Centrality")

# Add some metadata to the table.
table.description = description
table.location = f"Data from Figure 15, located on page 16."

# Add a thumbnail image of the original figure from the paper.
table.add_image(f"figure/v2_midrap_pdt_allCent_scaledEps2.pdf")

table.keywords["observables"] = ["V2/EPS2"]
table.keywords["reactions"]   = ["AU AU" , "AU AU --> P X" , "AU AU --> DEUT X" , "AU AU --> TRIT X" ]
table.keywords["cmenergies"]  = ["2.4"]
table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy", "Eccentricities", "Inclusive", "Transverse Momentum Dependence" , "Centrality Dependence"]

# x-axis: Centrality
varVnX = Variable("Centrality", is_independent=True, is_binned=False, units="%")
varVnX.values = centrality
table.add_variable(varVnX)

# Data for Particle and pt intervals
for ipart, particle in enumerate(PARTICLE):
 
 for imom, momentum in enumerate(MOMENTBINS):
     
     ### load data
     desc     = f"Data file: v2 pt{momentum} {PARTICLES[ipart]} scaledEps2"
     filename = f"data/HADES_v2_fit_au1230_pt{momentum}_{PARTICLES[ipart]}_scaledEps2.dat"
     data = np.loadtxt(filename)

     # y-axis
     varVnY = Variable(r"$v_{2} / \langle \epsilon_{2} \rangle$", is_independent=False, is_binned=False, units="")
     varVnY.values = data[:,1]

     varVnY.add_qualifier("Particle", particle)
     varVnY.add_qualifier("Pt [GeV/c]", MOMENTTEXT[imom])
 
     varVnYSys = Uncertainty("syst", is_symmetric=True)
     varVnYSys.values = data[:,2]

     varVnY.add_uncertainty(varVnYSys)
     table.add_variable(varVnY)
     table.add_additional_resource(desc, filename , copy_file=True)

# Add the table to the Submission object.
submission.add_table(table)

#######################################################################################
## Fig. 16 left

description = r"Same as in Figure 15, but for the scaled quadrangular flow. The values are divided by the square of second order eccentricity, $v_{4} / \langle \epsilon_{2} \rangle^{2}$ (left panel), and by the fourth order eccentricity, $v_{4} / \langle \epsilon_{4} \rangle$ (right panel)."

centrality = np.array([5 ,15 ,25 ,35])

# Define a table by the independent and dependent variables.
table = Table(f"Figure 16 v4eps22 Centrality")

# Add some metadata to the table.
table.description = description
table.location = f"Data from Figure 16 (left panel), located on page 16."

# Add a thumbnail image of the original figure from the paper.
table.add_image(f"figure/v4_midrap_pdt_allCent_scaledEps22.pdf")

table.keywords["observables"] = ["V4/EPS22"]
table.keywords["reactions"]   = ["AU AU" , "AU AU --> P X" , "AU AU --> DEUT X" , "AU AU --> TRIT X" ]
table.keywords["cmenergies"]  = ["2.4"]
table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy", "Eccentricities", "Inclusive", "Transverse Momentum Dependence" , "Centrality Dependence"]

# x-axis: Centrality
varVnX = Variable("Centrality", is_independent=True, is_binned=False, units="%")
varVnX.values = centrality
table.add_variable(varVnX)

# Data for Particle and pt intervals
for ipart, particle in enumerate(PARTICLE):
 
 for imom, momentum in enumerate(MOMENTBINS):
     
     ### load data
     desc     = f"Data file: v4 pt{momentum} {PARTICLES[ipart]} scaledEps22"
     filename = f"data/HADES_v4_fit_au1230_pt{momentum}_{PARTICLES[ipart]}_scaledEps22.dat"
     data = np.loadtxt(filename)

     # y-axis
     varVnY = Variable(r"$v_{4} / \langle \epsilon_{2} \rangle^{2}$", is_independent=False, is_binned=False, units="")
     varVnY.values = data[:,1]

     varVnY.add_qualifier("Particle", particle)
     varVnY.add_qualifier("Pt [GeV/c]", MOMENTTEXT[imom])
 
     varVnYSys = Uncertainty("syst", is_symmetric=True)
     varVnYSys.values = data[:,2]

     varVnY.add_uncertainty(varVnYSys)
     table.add_variable(varVnY)
     table.add_additional_resource(desc, filename , copy_file=True)

# Add the table to the Submission object.
submission.add_table(table)

#######################################################################################
## Fig. 16 right

# description = same as in upper table
centrality = np.array([5 ,15 ,25 ,35])

# Define a table by the independent and dependent variables.
table = Table(f"Figure 16 v4eps4 Centrality")

# Add some metadata to the table.
table.description = description
table.location = f"Data from Figure 16 (right panel), located on page 16."

# Add a thumbnail image of the original figure from the paper.
table.add_image(f"figure/v4_midrap_pdt_allCent_scaledEps4.pdf")

table.keywords["observables"] = ["V4/EPS4"]
table.keywords["reactions"]   = ["AU AU" , "AU AU --> P X" , "AU AU --> DEUT X" , "AU AU --> TRIT X" ]
table.keywords["cmenergies"]  = ["2.4"]
table.keywords["phrases"]     = ["Flow", "Azimuthal Anisotropy", "Eccentricities", "Inclusive", "Transverse Momentum Dependence" , "Centrality Dependence"]

# x-axis: Centrality
varVnX = Variable("Centrality", is_independent=True, is_binned=False, units="%")
varVnX.values = centrality
table.add_variable(varVnX)

# Data for Particle and pt intervals
for ipart, particle in enumerate(PARTICLE):
 
 for imom, momentum in enumerate(MOMENTBINS):
     
     ### load data
     desc     = f"Data file: v4 pt{momentum} {PARTICLES[ipart]} scaledEps4"
     filename = f"data/HADES_v4_fit_au1230_pt{momentum}_{PARTICLES[ipart]}_scaledEps4.dat"
     data = np.loadtxt(filename)

     # y-axis
     varVnY = Variable(r"$v_{4} / \langle \epsilon_{4} \rangle$", is_independent=False, is_binned=False, units="")
     varVnY.values = data[:,1]

     varVnY.add_qualifier("Particle", particle)
     varVnY.add_qualifier("Pt [GeV/c]", MOMENTTEXT[imom])
 
     varVnYSys = Uncertainty("syst", is_symmetric=True)
     varVnYSys.values = data[:,2]

     varVnY.add_uncertainty(varVnYSys)
     table.add_variable(varVnY)
     table.add_additional_resource(desc, filename , copy_file=True)

# Add the table to the Submission object.
submission.add_table(table)

#######################################################################################

outdir = "HEPdata_HADESflow"
submission.create_files(outdir, remove_old=True)