#include "Functions.C"

void Draw_graphs(TString inFileName="", TString outFileName="./test_graphs.root")
{
    TFile *fiCorr = new TFile(inFileName.Data(),"read");
    if (inFileName == "" || !fiCorr)
    {
        std::cerr << "No input file was provided!" << std::endl;
        return;
    }

    const std::pair<float, float> bcent_range = {4.06, 8.27};

    const int npid = 6; // h+-, proton, K+, K-, pi+, pi-
    const std::string pidname [] = {"hadrons", "proton", "kaon", "akaon", "pion", "apion"};
    const int ncorrsteps = 4;
    const std::string corrnames [ncorrsteps] = {"PLAIN", "RECENTERED", "TWIST", "RESCALED"};
    const std::string RPcorrname = "PLAIN";
    const int ndatatypes = 3;
    const std::string datatypes [ndatatypes] = {"mc", "reco", "mcreco"};
    // const std::string datatype = "mcreco";
    const int nProj = 2;
    const std::string projv2EP [nProj] = {"cos2", "sin2"};
    const std::string projv1EP [nProj] = {"cos1", "sin1"};
    const std::string projv2SP [nProj] = {"x2",   "y2"};
    const std::string projv1SP [nProj] = {"x1",   "y1"};

    std::string projPtName = "";
    std::string projYName = "";

    TFile *foGraphs = new TFile(outFileName.Data(),"recreate");

    Qn::DataContainerStatCalculate q2_tpcEP[nProj][nProj], q2_tpcSP[nProj][nProj];
    Qn::DataContainerStatCalculate q1_fhcalEP[nProj][nProj], q1_fhcalSP[nProj][nProj];
    Qn::DataContainerStatCalculate u2_tpcEP1[nProj][nProj][npid], u2_tpcEP2[nProj][nProj][npid];
    Qn::DataContainerStatCalculate u2_tpcSP1[nProj][nProj][npid], u2_tpcSP2[nProj][nProj][npid];
    // Qn::DataContainerStatCalculate u2_rpSP1[nProj][nProj][npid], u2_rpSP2[nProj][nProj][npid];
    Qn::DataContainerStatCalculate u2_fhcalEP1[nProj][nProj][nProj][npid], u2_fhcalEP2[nProj][nProj][nProj][npid];
    Qn::DataContainerStatCalculate u2_fhcalSP1[nProj][nProj][nProj][npid], u2_fhcalSP2[nProj][nProj][nProj][npid];
    Qn::DataContainerStatCalculate u2_rpSP1[nProj][nProj][nProj][npid], u2_rpSP2[nProj][nProj][nProj][npid];
    Qn::DataContainerStatCalculate u1_fhcalEP1[nProj][nProj][npid], u1_fhcalEP2[nProj][nProj][npid];
    Qn::DataContainerStatCalculate u1_fhcalSP1[nProj][nProj][npid], u1_fhcalSP2[nProj][nProj][npid];
    Qn::DataContainerStatCalculate u1_rpSP1[nProj][nProj][npid], u1_rpSP2[nProj][nProj][npid];
    TGraphErrors *graph;

    // Setting up the name of the Qn*Qn correlations
    //      Qn*Qn = "dir/nameQn1_correction.nameQn2_correction.component"
    //      For example: "QQ/EP/mc_TPC_EP_L_PLAIN.mc_TPC_EP_R_PLAIN.cos2cos2"
    std::string name;

    // Define the correction step that is in the input file
    std::string corrname = "";
    std::string datatype = "";
    for (int idtype=0; idtype<ndatatypes; idtype++)
    {
        for (int icorr=0; icorr<ncorrsteps; icorr++)
        {
            name = datatypes[idtype]+"_TPC_EP_L_"+corrnames[icorr]+"."+datatypes[idtype]+"_TPC_EP_R_"+corrnames[icorr]+".x2x2";
            if (!fiCorr->FindKeyAny(name.c_str())) continue;
            else
            {
                corrname = corrnames[icorr];
                datatype = datatypes[idtype];
                std::cout << "Data type found in the input file      : " << datatype << std::endl;
                std::cout << "Correction step found in the input file: " << corrname << std::endl;
            }
        }
    }
    if (corrname == "" || datatype == "")
    {
        std::cerr << "Couldn't find correction or datatype from the file. Exiting..." << std::endl;
        return;
    }


    if (datatype == "mc") projPtName = "_pT";
    if (datatype == "reco") projPtName = "_pT";
    if (datatype == "mcreco") projPtName = "_mc_pT";
    if (projPtName == "")
    {
        std::cerr << "Data type is not correctly set!" << std::endl;
        return;
    }

    if (datatype == "mc") projYName = "_rapidity";
    if (datatype == "reco") projYName = "_rapidity_pdg";
    if (datatype == "mcreco") projYName = "_mc_rapidity";
    if (projYName == "")
    {
        std::cerr << "Data type is not correctly set!" << std::endl;
        return;
    }

    std::string trackBranchName = "";
    if (datatype == "mc") trackBranchName = "McTracks";
    if (datatype == "reco") trackBranchName = "TpcTracks";
    if (datatype == "mcreco") trackBranchName = "TpcTracks";
    if (trackBranchName == "")
    {
        std::cerr << "Data type is not correctly set!" << std::endl;
        return;
    }


    // Getting correlations from the TFile
    // Qn*Qn
    for (int iproj1=0; iproj1<nProj; iproj1++)
    {
        for (int iproj2=0; iproj2<nProj; iproj2++)
        {
            name = "v2/QQ/EP/"+datatype+"_TPC_EP_L_"+corrname+"."+datatype+"_TPC_EP_R_"+corrname+"."+projv2EP[iproj1]+projv2EP[iproj2];
            if (!GetDCStatCalculate(fiCorr, name, q2_tpcEP[iproj1][iproj2]))
            {
                std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
            }
            name = "v2/QQ/SP/"+datatype+"_TPC_EP_L_"+corrname+"."+datatype+"_TPC_EP_R_"+corrname+"."+projv2SP[iproj1]+projv2SP[iproj2];
            if (!GetDCStatCalculate(fiCorr, name, q2_tpcSP[iproj1][iproj2]))
            {
                std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
            }
            name = "v1/QQ/EP/"+datatype+"_FHCal_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv1EP[iproj1]+projv1EP[iproj2];
            if (!GetDCStatCalculate(fiCorr, name, q1_fhcalEP[iproj1][iproj2]))
            {
                std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
            }
            name = "v1/QQ/SP/"+datatype+"_FHCal_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv1SP[iproj1]+projv1SP[iproj2];
            if (!GetDCStatCalculate(fiCorr, name, q1_fhcalSP[iproj1][iproj2]))
            {
                std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
            }
        }
    }

    // un*Qn
    for (int iproj1=0; iproj1<nProj; iproj1++)
    {
        for (int iproj2=0; iproj2<nProj; iproj2++)
        {
            for (int ipid=0; ipid<npid; ipid++)
            {
                name = "v2/uQ/EP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_TPC_EP_R_"+corrname+"."+projv2EP[iproj1]+projv2EP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u2_tpcEP1[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v2/uQ/EP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_TPC_EP_L_"+corrname+"."+projv2EP[iproj1]+projv2EP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u2_tpcEP2[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_TPC_EP_R_"+corrname+"."+projv2SP[iproj1]+projv2SP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u2_tpcSP1[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_TPC_EP_L_"+corrname+"."+projv2SP[iproj1]+projv2SP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u2_tpcSP2[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_RP_"+RPcorrname+"."+projv2SP[iproj1]+projv2SP[iproj2];
                // if (!GetDCStatCalculate(fiCorr, name, u2_rpSP1[iproj1][iproj2][ipid]))
                // {
                //     std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                // }
                // name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_RP_"+RPcorrname+"."+projv2SP[iproj1]+projv2SP[iproj2];
                // if (!GetDCStatCalculate(fiCorr, name, u2_rpSP2[iproj1][iproj2][ipid]))
                // {
                //     std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                // }

                name = "v1/uQ/EP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv1EP[iproj1]+projv1EP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u1_fhcalEP1[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v1/uQ/EP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_FHCal_L_"+corrname+"."+projv1EP[iproj1]+projv1EP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u1_fhcalEP2[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v1/uQ/SP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv1SP[iproj1]+projv1SP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u1_fhcalSP1[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v1/uQ/SP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_FHCal_L_"+corrname+"."+projv1SP[iproj1]+projv1SP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u1_fhcalSP2[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v1/uQ/SP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_RP_"+RPcorrname+"."+projv1SP[iproj1]+projv1SP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u1_rpSP1[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
                name = "v1/uQ/SP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_RP_"+RPcorrname+"."+projv1SP[iproj1]+projv1SP[iproj2];
                if (!GetDCStatCalculate(fiCorr, name, u1_rpSP2[iproj1][iproj2][ipid]))
                {
                    std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                }
            }
        }
    }

    // u2n*Qn*Qn
    for (int iproj1=0; iproj1<nProj; iproj1++)
    {
        for (int iproj2=0; iproj2<nProj; iproj2++)
        {
            for (int iproj3=0; iproj3<nProj; iproj3++)
            {
                for (int ipid=0; ipid<npid; ipid++)
                {
                    name = "v2/uQ/EP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_FHCal_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv2EP[iproj1]+projv1EP[iproj2]+projv1EP[iproj3];
                    if (!GetDCStatCalculate(fiCorr, name, u2_fhcalEP1[iproj1][iproj2][iproj3][ipid]))
                    {
                        std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                    }
                    name = "v2/uQ/EP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_FHCal_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv2EP[iproj1]+projv1EP[iproj2]+projv1EP[iproj3];
                    if (!GetDCStatCalculate(fiCorr, name, u2_fhcalEP2[iproj1][iproj2][iproj3][ipid]))
                    {
                        std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                    }
                    name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_FHCal_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv2SP[iproj1]+projv1SP[iproj2]+projv1SP[iproj3];
                    if (!GetDCStatCalculate(fiCorr, name, u2_fhcalSP1[iproj1][iproj2][iproj3][ipid]))
                    {
                        std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                    }
                    name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_FHCal_L_"+corrname+"."+datatype+"_FHCal_R_"+corrname+"."+projv2SP[iproj1]+projv1SP[iproj2]+projv1SP[iproj3];
                    if (!GetDCStatCalculate(fiCorr, name, u2_fhcalSP2[iproj1][iproj2][iproj3][ipid]))
                    {
                        std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                    }

                    name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_L_"+corrname+"."+datatype+"_RP_"+RPcorrname+"."+datatype+"_RP_"+RPcorrname+"."+projv2SP[iproj1]+projv1SP[iproj2]+projv1SP[iproj3];
                    if (!GetDCStatCalculate(fiCorr, name, u2_rpSP1[iproj1][iproj2][iproj3][ipid]))
                    {
                        std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                    }
                    name = "v2/uQ/SP/"+datatype+"_" + pidname[ipid] + "_R_"+corrname+"."+datatype+"_RP_"+RPcorrname+"."+datatype+"_RP_"+RPcorrname+"."+projv2SP[iproj1]+projv1SP[iproj2]+projv1SP[iproj3];
                    if (!GetDCStatCalculate(fiCorr, name, u2_rpSP2[iproj1][iproj2][iproj3][ipid]))
                    {
                        std::cerr << "Error: Cannot get " << name << " from the file!" << std::endl;
                    }
                }
            }
        }
    }

    // Taking into account that EP2/SP2 is rotated by pi for v1/Q1
    q1_fhcalEP[0][0] = -1. * q1_fhcalEP[0][0];
    q1_fhcalEP[1][1] = -1. * q1_fhcalEP[1][1];
    q1_fhcalSP[0][0] = -1. * q1_fhcalSP[0][0];
    q1_fhcalSP[1][1] = -1. * q1_fhcalSP[1][1];
    
    // Symmetrization for the pT-dependence of v1
    // for (int ipid=0; ipid<npid; ipid++)
    // {
    //     u1_fhcalEP2[0][0][ipid] = -1. * u1_fhcalEP2[0][0][ipid];
    //     u1_fhcalEP2[1][1][ipid] = -1. * u1_fhcalEP2[1][1][ipid];
    //     u1_fhcalSP2[0][0][ipid] = -1. * u1_fhcalSP2[0][0][ipid];
    //     u1_fhcalSP2[1][1][ipid] = -1. * u1_fhcalSP2[1][1][ipid];
    // }

    // Calculate resolution: res = Sqrt(<cos(Psi1 - Psi2)>) = Sqrt(<xx> + <yy>)
    Qn::DataContainerStatCalculate res2_tpcEP = Merge(q2_tpcEP[0][0], q2_tpcEP[1][1]);
    Qn::DataContainerStatCalculate res_tpcEP = Sqrt( res2_tpcEP );
    Qn::DataContainerStatCalculate res2_fhcalEP = Merge(q1_fhcalEP[0][0], q1_fhcalEP[1][1]);
    Qn::DataContainerStatCalculate res_fhcalEP = Sqrt( res2_fhcalEP );

    Qn::DataContainerStatCalculate res2_tpcSP = Merge(q2_tpcSP[0][0], q2_tpcSP[1][1]);
    Qn::DataContainerStatCalculate res_tpcSP = Sqrt( res2_tpcSP );
    Qn::DataContainerStatCalculate res2_fhcalSP = Merge(q1_fhcalSP[0][0], q1_fhcalSP[1][1]);
    Qn::DataContainerStatCalculate res_fhcalSP = Sqrt( res2_fhcalSP );

    foGraphs->mkdir("res");
    foGraphs->cd("res");
    name = "res_v2_tpcEP_" + corrname;
    graph = Qn::DataContainerHelper::ToTGraph(res_tpcEP);
    graph->SetName(Form("%s", name.c_str()));
    name += ";b, fm;Resolution";
    graph->SetTitle(Form("%s", name.c_str()));
    graph->Write();
    name = "res_v1_fhcalEP_" + corrname;
    graph = Qn::DataContainerHelper::ToTGraph(res_fhcalEP);
    graph->SetName(Form("%s", name.c_str()));
    name += ";b, fm;Resolution";
    graph->SetTitle(Form("%s", name.c_str()));
    graph->Write();

    foGraphs->mkdir("resSP");
    foGraphs->cd("resSP");
    name = "res_v2_tpcSP_" + corrname;
    graph = Qn::DataContainerHelper::ToTGraph(res_tpcSP);
    graph->SetName(Form("%s", name.c_str()));
    name += ";b, fm;Resolution";
    graph->SetTitle(Form("%s", name.c_str()));
    graph->Write();
    name = "res_v1_fhcalSP_" + corrname;
    graph = Qn::DataContainerHelper::ToTGraph(res_fhcalSP);
    graph->SetName(Form("%s", name.c_str()));
    name += ";b, fm;Resolution";
    graph->SetTitle(Form("%s", name.c_str()));
    graph->Write();

    Qn::DataContainerStatCalculate vn_pT_obs, vn_pT_div, vn_pT_cent, vn_pT;
    Qn::DataContainerStatCalculate vn_y_obs, vn_y_div, vn_y_cent, vn_y;
    Qn::DataContainerStatCalculate uq_cent, uq_pT, uq_y;
    Qn::DataContainerStatCalculate vn_comp_pT_obs, vn_comp_pT_div, vn_comp_pT_cent, vn_comp_pT;
    Qn::DataContainerStatCalculate vn_comp_y_obs, vn_comp_y_div, vn_comp_y_cent, vn_comp_y;

    // v2 TPC EP
    foGraphs->mkdir("v2tpcEP");
    foGraphs->cd("v2tpcEP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        vn_pT_obs = Merge(Merge(u2_tpcEP1[0][0][ipid], u2_tpcEP1[1][1][ipid]), Merge(u2_tpcEP2[0][0][ipid], u2_tpcEP2[1][1][ipid]));
        vn_pT_div = vn_pT_obs / res_tpcEP;
        vn_pT_cent = vn_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projPtName;
        vn_pT = vn_pT_cent.Projection({name});
        vn_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v2_pT_tpcEP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_pT);
        graph->SetName(name.c_str());
        name += ";p_{T}, GeV/c;v_{2}";
        graph->SetTitle(name.c_str());
        graph->Write();
        vn_y_obs = Merge(Merge(u2_tpcEP1[0][0][ipid], u2_tpcEP1[1][1][ipid]), Merge(u2_tpcEP2[0][0][ipid], u2_tpcEP2[1][1][ipid]));
        vn_y_div = vn_y_obs / res_tpcEP;
        vn_y_cent = vn_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projYName;
        vn_y = vn_y_cent.Projection({name});
        vn_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v2_y_tpcEP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_y);
        graph->SetName(name.c_str());
        name += ";y;v_{2}";
        graph->SetTitle(name.c_str());
        graph->Write();
    }

    // v2 TPC SP
    foGraphs->mkdir("v2tpcSP");
    foGraphs->cd("v2tpcSP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        vn_pT_obs = Merge(Merge(u2_tpcSP1[0][0][ipid], u2_tpcSP1[1][1][ipid]), Merge(u2_tpcSP2[0][0][ipid], u2_tpcSP2[1][1][ipid]));
        vn_pT_div = vn_pT_obs / res_tpcSP;
        vn_pT_cent = vn_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projPtName;
        vn_pT = vn_pT_cent.Projection({name});
        vn_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v2_pT_tpcSP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_pT);
        graph->SetName(name.c_str());
        name += ";p_{T}, GeV/c;v_{2}";
        graph->SetTitle(name.c_str());
        graph->Write();
        vn_y_obs = Merge(Merge(u2_tpcSP1[0][0][ipid], u2_tpcSP1[1][1][ipid]), Merge(u2_tpcSP2[0][0][ipid], u2_tpcSP2[1][1][ipid]));
        vn_y_div = vn_y_obs / res_tpcSP;
        vn_y_cent = vn_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projYName;
        vn_y = vn_y_cent.Projection({name});
        vn_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v2_y_tpcSP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_y);
        graph->SetName(name.c_str());
        name += ";y;v_{2}";
        graph->SetTitle(name.c_str());
        graph->Write();
    }

    // v1 FHCal EP
    foGraphs->mkdir("v1fhcalEP");
    foGraphs->cd("v1fhcalEP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        vn_pT_obs = Merge(Merge(u1_fhcalEP1[0][0][ipid], u1_fhcalEP1[1][1][ipid]), Merge(u1_fhcalEP2[0][0][ipid], u1_fhcalEP2[1][1][ipid]));
        if (corrname == "mc")
        {
            vn_pT_obs = -1.*vn_pT_obs;
        }
        vn_pT_div = vn_pT_obs / res_fhcalEP;
        vn_pT_cent = vn_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projPtName;
        vn_pT = vn_pT_cent.Projection({name});
        vn_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v1_pT_fhcalEP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_pT);
        graph->SetName(name.c_str());
        name += ";p_{T}, GeV/c;v_{1}";
        graph->SetTitle(name.c_str());
        graph->Write();
        vn_y_obs = Merge(Merge(-1.*u1_fhcalEP1[0][0][ipid], -1.*u1_fhcalEP1[1][1][ipid]), Merge(u1_fhcalEP2[0][0][ipid], u1_fhcalEP2[1][1][ipid]));
        if (corrname == "mc")
        {
            vn_y_obs = -1.*vn_y_obs;
        }
        vn_y_div = vn_y_obs / res_fhcalEP;
        vn_y_cent = vn_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projYName;
        vn_y = vn_y_cent.Projection({name});
        vn_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v1_y_fhcalEP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_y);
        graph->SetName(name.c_str());
        name += ";y;v_{1}";
        graph->SetTitle(name.c_str());
        graph->Write();
    }

    // v1 FHCal SP
    foGraphs->mkdir("v1fhcalSP");
    foGraphs->cd("v1fhcalSP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        vn_pT_obs = Merge(Merge(u1_fhcalSP1[0][0][ipid], u1_fhcalSP1[1][1][ipid]), Merge(u1_fhcalSP2[0][0][ipid], u1_fhcalSP2[1][1][ipid]));
        if (datatype == "mc")
        {
            vn_pT_obs = -1. * vn_pT_obs;
        }
        vn_pT_div = vn_pT_obs / res_fhcalSP;
        vn_pT_cent = vn_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projPtName;
        vn_pT = vn_pT_cent.Projection({name});
        vn_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v1_pT_fhcalSP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_pT);
        graph->SetName(name.c_str());
        name += ";p_{T}, GeV/c;v_{1}";
        graph->SetTitle(name.c_str());
        graph->Write();
        vn_y_obs = Merge(Merge(-1.*u1_fhcalSP1[0][0][ipid], -1.*u1_fhcalSP1[1][1][ipid]), Merge(u1_fhcalSP2[0][0][ipid], u1_fhcalSP2[1][1][ipid]));
        if (datatype == "mc")
        {
            vn_y_obs = -1. * vn_y_obs;
        }
        vn_y_div = vn_y_obs / res_fhcalSP;
        vn_y_cent = vn_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projYName;
        vn_y = vn_y_cent.Projection({name});
        vn_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v1_y_fhcalSP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_y);
        graph->SetName(name.c_str());
        name += ";y;v_{1}";
        graph->SetTitle(name.c_str());
        graph->Write();
    }

    // v1 RP SP
    foGraphs->mkdir("v1rpSP");
    foGraphs->cd("v1rpSP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        vn_pT_obs = Merge(Merge(-1.*u1_rpSP1[0][0][ipid], -1.*u1_rpSP1[1][1][ipid]), Merge(u1_rpSP2[0][0][ipid], u1_rpSP2[1][1][ipid]));
        vn_pT_div = vn_pT_obs;
        vn_pT_cent = vn_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projPtName;
        vn_pT = vn_pT_cent.Projection({name});
        vn_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v1_pT_rpSP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_pT);
        graph->SetName(name.c_str());
        name += ";p_{T}, GeV/c;v_{1}";
        graph->SetTitle(name.c_str());
        graph->Write();
        vn_y_obs = Merge(Merge(u1_rpSP1[0][0][ipid], u1_rpSP1[1][1][ipid]), Merge(u1_rpSP2[0][0][ipid], u1_rpSP2[1][1][ipid]));
        vn_y_div = vn_y_obs;
        vn_y_cent = vn_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
        name = trackBranchName+projYName;
        vn_y = vn_y_cent.Projection({name});
        vn_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
        name = "v1_y_rpSP_" + pidname[ipid] + "_" + corrname;
        graph = Qn::DataContainerHelper::ToTGraph(vn_y);
        graph->SetName(name.c_str());
        name += ";y;v_{1}";
        graph->SetTitle(name.c_str());
        graph->Write();
    }

    // QQ correlations
    foGraphs->mkdir("QQ");
    foGraphs->cd("QQ");
    for (int iproj1=0; iproj1<nProj; iproj1++)
    {
        for (int iproj2=0; iproj2<nProj; iproj2++)
        {
            name = "QQ_tpc_"+ corrname + "_"+projv2EP[iproj1]+projv2EP[iproj2];
            graph = Qn::DataContainerHelper::ToTGraph(q2_tpcEP[iproj1][iproj2]);
            graph->SetName(name.c_str());
            name += ";b, fm;<"+projv2EP[iproj1]+projv2EP[iproj2]+">";
            graph->SetTitle(name.c_str());
            graph->Write();

            name = "QQ_tpc_"+ corrname + "_"+projv2SP[iproj1]+projv2SP[iproj2];
            graph = Qn::DataContainerHelper::ToTGraph(q2_tpcSP[iproj1][iproj2]);
            graph->SetName(name.c_str());
            name += ";b, fm;<"+projv2SP[iproj1]+projv2SP[iproj2]+">";
            graph->SetTitle(name.c_str());
            graph->Write();

            name = "QQ_fhcal_"+ corrname + "_"+projv1EP[iproj1]+projv1EP[iproj2];
            graph = Qn::DataContainerHelper::ToTGraph(q1_fhcalEP[iproj1][iproj2]);
            graph->SetName(name.c_str());
            name += ";b, fm;<"+projv1EP[iproj1]+projv1EP[iproj2]+">";
            graph->SetTitle(name.c_str());
            graph->Write();

            name = "QQ_fhcal_"+ corrname + "_"+projv1SP[iproj1]+projv1SP[iproj2];
            graph = Qn::DataContainerHelper::ToTGraph(q1_fhcalSP[iproj1][iproj2]);
            graph->SetName(name.c_str());
            name += ";b, fm;<"+projv1SP[iproj1]+projv1SP[iproj2]+">";
            graph->SetTitle(name.c_str());
            graph->Write();
        }
    }

    // uQ/QQ correlations
    foGraphs->mkdir("v1Components_EP");
    foGraphs->cd("v1Components_EP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        for (int iproj1=0; iproj1<nProj; iproj1++)
        {
            for (int iproj2=0; iproj2<nProj; iproj2++)
            {
                vn_comp_pT_obs = Merge(u1_fhcalEP1[iproj1][iproj2][ipid], u1_fhcalEP2[iproj1][iproj2][ipid]);
                vn_comp_pT_div = vn_comp_pT_obs / Sqrt(q1_fhcalEP[iproj1][iproj2]);
                vn_comp_pT_cent = vn_comp_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                vn_comp_pT = vn_comp_pT_cent.Projection({name});
                vn_comp_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v1_pT_fhcalEP_" + pidname[ipid] + "_" + corrname +"_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;v_{1}";
                graph->SetTitle(name.c_str());
                graph->Write();
                vn_comp_y_obs = Merge(-1.*u1_fhcalEP1[iproj1][iproj2][ipid], u1_fhcalEP2[iproj1][iproj2][ipid]);
                vn_comp_y_div = vn_comp_y_obs / Sqrt(q1_fhcalEP[iproj1][iproj2]);
                vn_comp_y_cent = vn_comp_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projYName;
                vn_comp_y = vn_comp_y_cent.Projection({name});
                vn_comp_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v1_y_fhcalEP_" + pidname[ipid] + "_" + corrname +"_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_y);
                graph->SetName(name.c_str());
                name += ";y;v_{1}";
                graph->SetTitle(name.c_str());
                graph->Write();
            }
        }
    }

    foGraphs->mkdir("v2Components_EP");
    foGraphs->cd("v2Components_EP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        for (int iproj1=0; iproj1<nProj; iproj1++)
        {
            for (int iproj2=0; iproj2<nProj; iproj2++)
            {
                vn_comp_pT_obs = Merge(u1_fhcalEP1[iproj1][iproj2][ipid], u1_fhcalEP2[iproj1][iproj2][ipid]);
                vn_comp_pT_div = vn_comp_pT_obs / Sqrt(q1_fhcalEP[iproj1][iproj2]);
                vn_comp_pT_cent = vn_comp_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                vn_comp_pT = vn_comp_pT_cent.Projection({name});
                vn_comp_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v1_pT_fhcalEP_" + pidname[ipid] + "_" + corrname +"_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;v_{1}";
                graph->SetTitle(name.c_str());
                graph->Write();
                vn_comp_y_obs = Merge(-1.*u1_fhcalEP1[iproj1][iproj2][ipid], u1_fhcalEP2[iproj1][iproj2][ipid]);
                vn_comp_y_div = vn_comp_y_obs / Sqrt(q1_fhcalEP[iproj1][iproj2]);
                vn_comp_y_cent = vn_comp_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projYName;
                vn_comp_y = vn_comp_y_cent.Projection({name});
                vn_comp_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v1_y_fhcalEP_" + pidname[ipid] + "_" + corrname +"_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_y);
                graph->SetName(name.c_str());
                name += ";y;v_{1}";
                graph->SetTitle(name.c_str());
                graph->Write();
            }
        }
    }

    foGraphs->mkdir("v1Components_SP");
    foGraphs->cd("v1Components_SP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        for (int iproj1=0; iproj1<nProj; iproj1++)
        {
            for (int iproj2=0; iproj2<nProj; iproj2++)
            {
                vn_comp_pT_obs = Merge(u1_fhcalSP1[iproj1][iproj2][ipid], u1_fhcalSP2[iproj1][iproj2][ipid]);
                vn_comp_pT_div = vn_comp_pT_obs / Sqrt(q1_fhcalSP[iproj1][iproj2]);
                vn_comp_pT_cent = vn_comp_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                vn_comp_pT = vn_comp_pT_cent.Projection({name});
                vn_comp_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v1_pT_fhcalSP_" + pidname[ipid] + "_" + corrname +"_"+projv1SP[iproj1]+projv1SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;v_{1}";
                graph->SetTitle(name.c_str());
                graph->Write();
                vn_comp_y_obs = Merge(-1.*u1_fhcalSP1[iproj1][iproj2][ipid], u1_fhcalSP2[iproj1][iproj2][ipid]);
                vn_comp_y_div = vn_comp_y_obs / Sqrt(q1_fhcalSP[iproj1][iproj2]);
                vn_comp_y_cent = vn_comp_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projYName;
                vn_comp_y = vn_comp_y_cent.Projection({name});
                vn_comp_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v1_y_fhcalSP_" + pidname[ipid] + "_" + corrname +"_"+projv1SP[iproj1]+projv1SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_y);
                graph->SetName(name.c_str());
                name += ";y;v_{1}";
                graph->SetTitle(name.c_str());
                graph->Write();
            }
        }
    }

    foGraphs->mkdir("v2Components_SP");
    foGraphs->cd("v2Components_SP");
    for (int ipid=0; ipid<npid; ipid++)
    {
        for (int iproj1=0; iproj1<nProj; iproj1++)
        {
            for (int iproj2=0; iproj2<nProj; iproj2++)
            {
                vn_comp_pT_obs = Merge(u2_tpcSP1[iproj1][iproj2][ipid], u2_tpcSP2[iproj1][iproj2][ipid]);
                vn_comp_pT_div = vn_comp_pT_obs / Sqrt(q2_tpcSP[iproj1][iproj2]);
                vn_comp_pT_cent = vn_comp_pT_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                vn_comp_pT = vn_comp_pT_cent.Projection({name});
                vn_comp_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v2_pT_tpcSP_" + pidname[ipid] + "_" + corrname +"_"+projv2SP[iproj1]+projv2SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;v_{2}";
                graph->SetTitle(name.c_str());
                graph->Write();
                vn_comp_y_obs = Merge(u2_tpcSP1[iproj1][iproj2][ipid], u2_tpcSP2[iproj1][iproj2][ipid]);
                vn_comp_y_div = vn_comp_y_obs / Sqrt(q2_tpcSP[iproj1][iproj2]);
                vn_comp_y_cent = vn_comp_y_div.Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projYName;
                vn_comp_y = vn_comp_y_cent.Projection({name});
                vn_comp_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "v2_y_tpcSP_" + pidname[ipid] + "_" + corrname +"_"+projv2SP[iproj1]+projv2SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(vn_comp_y);
                graph->SetName(name.c_str());
                name += ";y;v_{2}";
                graph->SetTitle(name.c_str());
                graph->Write();
            }
        }
    }

    // uQ correlations
    foGraphs->mkdir("uQ");
    foGraphs->cd("uQ");
    for (int ipid=0; ipid<npid; ipid++)
    {
        for (int iproj1=0; iproj1<nProj; iproj1++)
        {
            for (int iproj2=0; iproj2<nProj; iproj2++)
            {
                uq_cent = u2_tpcEP1[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_tpcEP1_" + pidname[ipid] + "_" + corrname + "_"+projv2EP[iproj1]+projv2EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv2EP[iproj1]+projv2EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_tpcEP1_" + pidname[ipid] + "_" + corrname + "_"+projv2EP[iproj1]+projv2EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv2EP[iproj1]+projv2EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u2_tpcEP2[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_tpcEP2_" + pidname[ipid] + "_" + corrname + "_"+projv2EP[iproj1]+projv2EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv2EP[iproj1]+projv2EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_tpcEP2_" + pidname[ipid] + "_" + corrname + "_"+projv2EP[iproj1]+projv2EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv2EP[iproj1]+projv2EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u2_tpcSP1[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_tpcSP1_" + pidname[ipid] + "_" + corrname + "_"+projv2SP[iproj1]+projv2SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv2SP[iproj1]+projv2SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_tpcSP1_" + pidname[ipid] + "_" + corrname + "_"+projv2SP[iproj1]+projv2SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv2SP[iproj1]+projv2SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u2_tpcSP2[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_tpcSP2_" + pidname[ipid] + "_" + corrname + "_"+projv2SP[iproj1]+projv2SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv2SP[iproj1]+projv2SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_tpcSP2_" + pidname[ipid] + "_" + corrname + "_"+projv2SP[iproj1]+projv2SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv2SP[iproj1]+projv2SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u1_fhcalEP1[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_fhcalEP1_" + pidname[ipid] + "_" + corrname + "_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv1EP[iproj1]+projv1EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_fhcalEP1_" + pidname[ipid] + "_" + corrname + "_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv1EP[iproj1]+projv1EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u1_fhcalEP2[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_fhcalEP2_" + pidname[ipid] + "_" + corrname + "_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv1EP[iproj1]+projv1EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_fhcalEP2_" + pidname[ipid] + "_" + corrname + "_"+projv1EP[iproj1]+projv1EP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv1EP[iproj1]+projv1EP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u1_fhcalSP1[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_fhcalSP1_" + pidname[ipid] + "_" + corrname + "_"+projv1SP[iproj1]+projv1SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv1SP[iproj1]+projv1SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_fhcalSP1_" + pidname[ipid] + "_" + corrname + "_"+projv1SP[iproj1]+projv1SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv1SP[iproj1]+projv1SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();

                uq_cent = u1_fhcalSP2[iproj1][iproj2][ipid].Rebin({"McEvent_B", 1, bcent_range.first, bcent_range.second});
                name = trackBranchName+projPtName;
                uq_pT = uq_cent.Projection({name});
                uq_pT.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_pT_fhcalSP2_" + pidname[ipid] + "_" + corrname + "_"+projv1SP[iproj1]+projv1SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_pT);
                graph->SetName(name.c_str());
                name += ";p_{T}, GeV/c;<"+projv1SP[iproj1]+projv1SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
                name = trackBranchName+projYName;
                uq_y = uq_cent.Projection({name});
                uq_y.SetErrors(Qn::Stat::ErrorType::BOOTSTRAP);
                name = "uQ_y_fhcalSP2_" + pidname[ipid] + "_" + corrname + "_"+projv1SP[iproj1]+projv1SP[iproj2];
                graph = Qn::DataContainerHelper::ToTGraph(uq_y);
                graph->SetName(name.c_str());
                name += ";y;<"+projv1SP[iproj1]+projv1SP[iproj2]+">";
                graph->SetTitle(name.c_str());
                graph->Write();
            }
        }
    }

    foGraphs->Close();
}