Browse Source

First commit

PeterParfenov 2 months ago
commit
08ff5bf904
6 changed files with 3173 additions and 0 deletions
  1. 16 0
      CutsMC.C
  2. 839 0
      Draw_graphs.C
  3. 42 0
      Functions.C
  4. 547 0
      MpdDst2AT.C
  5. 1075 0
      mpd-analysis-config.yml
  6. 654 0
      mpd-correlation-plain.yml

+ 16 - 0
CutsMC.C

@@ -0,0 +1,16 @@
+#include "CutsRegistry.hpp"
+
+int CutsMC() {
+    using AnalysisTree::RegisterCuts;
+    using namespace AnalysisTree::Version1;
+
+    SimpleCut vtx_z(Variable("McEvent","vtx_z"), -500, 500);
+    SimpleCut b_mc(Variable("McEvent","B"), 0., 14.);
+
+    RegisterCuts("mpd/auau/7gev/mc", Cuts("McEvent", {
+                vtx_z,
+                b_mc
+    }));
+
+    return 0;
+}

+ 839 - 0
Draw_graphs.C

@@ -0,0 +1,839 @@
+#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();
+}

+ 42 - 0
Functions.C

@@ -0,0 +1,42 @@
+// Function to read either Qn::DataContainerStatCollect or Qn::DataContainerStatCalculate from the file into dc_calc
+Bool_t GetDCStatCalculate(TFile *const& file, std::string name, Qn::DataContainerStatCalculate& dc_calc)
+{
+    Qn::DataContainerStatCollect *tmp_dc_collect{nullptr};
+    Qn::DataContainerStatCalculate *tmp_dc_calculate{nullptr};
+
+    file->GetObject(name.c_str(), tmp_dc_calculate);
+    if (!tmp_dc_calculate)
+    {
+        file->GetObject(name.c_str(), tmp_dc_collect);
+        if (!tmp_dc_collect)
+        {
+            std::cerr << "Cannot get the object " << name << "!" << std::endl;
+            return false;
+        }
+        dc_calc = (Qn::DataContainerStatCalculate) *tmp_dc_collect;
+    }
+    else
+    {
+        dc_calc = (Qn::DataContainerStatCalculate) *tmp_dc_calculate;
+    }
+
+    return true;
+}
+
+// Returns mean value from dc1 and dc2
+Qn::DataContainerStatCalculate Merge(Qn::DataContainerStatCalculate dc1, Qn::DataContainerStatCalculate dc2)
+{
+    Qn::DataContainerStatCalculate result = dc1;
+    if (dc1.size() != dc2.size())
+    {
+        std::cerr << "Merge: dc1 and dc2 has different sizes!" << std::endl;
+        return result;
+    }
+
+    for (int i=0; i<result.size(); i++)
+    {
+        result.At(i) = Qn::Merge(dc1.At(i), dc2.At(i));
+    }
+
+    return result;
+}

+ 547 - 0
MpdDst2AT.C

@@ -0,0 +1,547 @@
+/*// Standard c++ headers
+#include <iostream>
+#include <string>
+#include <utility>
+#include <set>
+#include <map>
+#include <functional>
+
+// Standard ROOT headers
+#include <Rtypes.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TString.h>
+#include <TStopwatch.h>
+#include <TClonesArray.h>
+#include <TObject.h>
+#include <TMath.h>
+#include <TDatabasePDG.h>
+#include <TParticlePDG.h>
+
+// FairRoot/MpdRoot headers
+#include "FairMCEventHeader.h"
+#ifdef _MCSTACK_
+#include "FairMCTrack.h"
+#endif
+#ifdef _MPDMCSTACK_
+#include "MpdMCTrack.h"
+#endif
+#include "MpdEvent.h"
+#include "MpdZdcDigi.h"
+#include "MpdPid.h"
+#include "MpdTrack.h"
+#include "MpdKalmanTrack.h"
+#include "MpdVertex.h"
+
+// AnalysisTree headers
+#include "AnalysisTree/Configuration.hpp"
+#include "AnalysisTree/DataHeader.hpp"
+#include "AnalysisTree/EventHeader.hpp"
+#include "AnalysisTree/Detector.hpp"
+#include "AnalysisTree/Matching.hpp"
+*/
+
+#define _MPDMCSTACK_
+
+float get_beamP(float sqrtSnn, float m_target = 0.938, float m_beam = 0.938);
+float GetFHCalPhi(int iModule);
+TVector3 GetFHCalPos(int iModule);
+Bool_t IsCharged(Int_t pdg);
+Float_t GetRapidity(Float_t p, Float_t pz, Int_t pid_type);
+Float_t GetRapidityPDG(Float_t p, Float_t pz, Int_t pdg);
+
+void MpdDst2AT(TString iFileName, TString oFileName, std::string system = "", float sqrtSnn = -1.)
+{
+  bool isDataHeaderDefined = false;
+  if (system != "" && sqrtSnn != -1) isDataHeaderDefined = true;
+
+  TStopwatch timer;
+  timer.Start();
+
+  // Open input dst tree
+  TFile *fi = new TFile(iFileName.Data(),"read");
+  if (!fi || fi->IsZombie())
+  {
+    std::cerr << "ERROR: Input file probably is empty. Exit the program." << std::endl;
+    return 100;
+  }
+
+  TTree  *dstTree = (TTree*) fi->Get("mpdsim");
+
+  // PID related parameters
+  const Double_t PIDsigM = 4.0;
+  const Double_t PIDsigE = 4.0;
+  const Double_t PIDenergy = 11.;
+  const Double_t PIDkoeff = 1.;
+  const TString PIDgenerator = "URQMD";
+  const TString PIDtracking = "CF";
+  const TString PIDparticles = "pikapr";
+  const Int_t   Num_Of_Modules = 90;
+
+  MpdPid *pid = new MpdPid(PIDsigM, PIDsigE, PIDenergy, PIDkoeff, PIDgenerator, PIDtracking, PIDparticles);
+
+  // Prepare to read input dst
+  FairMCEventHeader *MCHeader;
+  TClonesArray *MCTracks;
+  MpdEvent *MPDEvent;
+  TClonesArray *FHCalHits;
+  TClonesArray *MpdGlobalTracks;
+  MpdZdcDigi *FHCalHit;
+  //TClonesArray *mpdKalmanTracks;
+  TClonesArray *vertexes;
+
+  MCHeader = nullptr;
+  MCTracks = nullptr;
+  MPDEvent = nullptr;
+  FHCalHits = nullptr;
+  //mpdKalmanTracks = nullptr;
+  vertexes = nullptr;
+
+  dstTree->SetBranchAddress("MCEventHeader.", &MCHeader);
+  dstTree->SetBranchAddress("MCTrack", &MCTracks);
+  dstTree->SetBranchAddress("MPDEvent.", &MPDEvent);
+  dstTree->SetBranchAddress("ZdcDigi", &FHCalHits);
+  //dstTree->SetBranchAddress("TpcKalmanTrack", &mpdKalmanTracks);
+  dstTree->SetBranchAddress("Vertex", &vertexes);
+
+  // Set up output dst
+  TFile *outFile = new TFile(oFileName.Data(), "RECREATE");
+  TTree *outTree = new TTree("aTree","AnalysisTree Dst at MPD");
+
+  // Set up AnalysisTree data header
+  AnalysisTree::DataHeader *dataHeader = new AnalysisTree::DataHeader;
+  if (system != "") dataHeader->SetSystem(system);
+  if (sqrtSnn > 0) dataHeader->SetBeamMomentum(get_beamP(sqrtSnn));
+  auto &fhcal_mod_pos = dataHeader->AddDetector();
+  TVector3 modulePos;
+  for (int imodule=0; imodule<Num_Of_Modules; imodule++)
+  {
+    auto *module = fhcal_mod_pos.AddChannel();
+    modulePos = GetFHCalPos(imodule);
+    module->SetPosition(modulePos);
+  }
+
+  // Set up AnalysisTree configureation
+  AnalysisTree::Configuration *out_config = new AnalysisTree::Configuration;
+
+  // Set up AnalysisTree configuration
+  std::string str_reco_event_branch = "RecoEvent";
+  std::string str_mc_event_branch = "McEvent";
+  std::string str_tpc_tracks_branch = "TpcTracks";
+  std::string str_fhcal_branch = "FHCalModules";
+  std::string str_mc_tracks_branch = "McTracks";
+  std::string str_tpc2mc_tracks_branch = "TpcTracks2McTracks";
+
+  AnalysisTree::BranchConfig reco_event_branch(str_reco_event_branch.c_str(), AnalysisTree::DetType::kEventHeader); 
+  AnalysisTree::BranchConfig mc_event_branch(str_mc_event_branch.c_str(), AnalysisTree::DetType::kEventHeader); 
+  AnalysisTree::BranchConfig fhcal_branch(str_fhcal_branch.c_str(), AnalysisTree::DetType::kModule); 
+  AnalysisTree::BranchConfig tpc_tracks_branch(str_tpc_tracks_branch.c_str(), AnalysisTree::DetType::kTrack); 
+  AnalysisTree::BranchConfig mc_tracks_branch(str_mc_tracks_branch.c_str(), AnalysisTree::DetType::kParticle);
+
+  mc_event_branch.AddField<float>("B");
+  mc_event_branch.AddField<float>("PhiRp");
+  // mc_event_branch.AddField<float>("Centrality_B");
+
+  // reco_event_branch.AddField<float>("refMult");
+  // reco_event_branch.AddField<float>("Centrality_refMult");
+
+  tpc_tracks_branch.AddField<int>("nhits");
+  tpc_tracks_branch.AddField<int>("nhits_poss");
+  tpc_tracks_branch.AddField<int>("charge");
+  tpc_tracks_branch.AddFields<float>({"dca_x","dca_y","dca_z"});
+  tpc_tracks_branch.AddField<float>("chi2");
+  tpc_tracks_branch.AddField<float>("tof_mass2");
+  tpc_tracks_branch.AddField<float>("tof_flag");
+  tpc_tracks_branch.AddField<float>("dedx");
+  tpc_tracks_branch.AddField<float>("pid_prob_pion");
+  tpc_tracks_branch.AddField<float>("pid_prob_kaon");
+  tpc_tracks_branch.AddField<float>("pid_prob_proton");
+  tpc_tracks_branch.AddField<float>("mc_E");
+  tpc_tracks_branch.AddField<float>("mc_pT");
+  tpc_tracks_branch.AddField<float>("mc_eta");
+  tpc_tracks_branch.AddField<float>("mc_phi");
+  tpc_tracks_branch.AddField<float>("mc_rapidity");
+  tpc_tracks_branch.AddField<int>("mc_mother_id");
+  tpc_tracks_branch.AddField<int>("mc_pdg");
+  tpc_tracks_branch.AddField<int>("eta_sign");
+  tpc_tracks_branch.AddField<float>("rapidity_pdg");
+  tpc_tracks_branch.AddField<float>("rapidity_pion");
+  tpc_tracks_branch.AddField<float>("rapidity_kaon");
+  tpc_tracks_branch.AddField<float>("rapidity_proton");
+
+  fhcal_branch.AddField<float>("phi");
+  fhcal_branch.AddField<float>("signal_eta_signed");
+
+  mc_tracks_branch.AddField<int>("mother_id");
+  mc_tracks_branch.AddField<bool>("is_charged");
+  mc_tracks_branch.AddField<int>("eta_sign");
+
+  auto hasher = std::hash<std::string>();
+
+  // Initialize AnalysisTree Dst components
+  out_config->AddBranchConfig(reco_event_branch);
+  AnalysisTree::EventHeader *reco_event = new AnalysisTree::EventHeader( Short_t(hasher(reco_event_branch.GetName())) );
+  out_config->AddBranchConfig(mc_event_branch);
+  AnalysisTree::EventHeader *mc_event = new AnalysisTree::EventHeader( Short_t(hasher(mc_event_branch.GetName())) );
+  out_config->AddBranchConfig(tpc_tracks_branch);
+  AnalysisTree::TrackDetector *tpc_tracks = new AnalysisTree::TrackDetector( Short_t(hasher(tpc_tracks_branch.GetName())) ); 
+  out_config->AddBranchConfig(fhcal_branch);
+  AnalysisTree::ModuleDetector *fhcal_modules = new AnalysisTree::ModuleDetector( Short_t(hasher(fhcal_branch.GetName())) );
+  out_config->AddBranchConfig(mc_tracks_branch);
+  AnalysisTree::Particles *mc_tracks = new AnalysisTree::Particles( Short_t(hasher(mc_tracks_branch.GetName())) ); 
+  AnalysisTree::Matching *tpc2mc_tracks = new AnalysisTree::Matching(out_config->GetBranchConfig(str_tpc_tracks_branch).GetId(), out_config->GetBranchConfig(str_mc_tracks_branch).GetId());
+  out_config->AddMatch(tpc2mc_tracks);
+
+  reco_event->Init(reco_event_branch);
+  mc_event->Init(mc_event_branch);
+
+  // mc_event's additional field ids
+  //const int mceventid = mc_event->GetId();
+  const int iB = out_config->GetBranchConfig(str_mc_event_branch).GetFieldId("B");
+  const int iPhiRp = out_config->GetBranchConfig(str_mc_event_branch).GetFieldId("PhiRp");
+
+  // fhcal_modules' additional field ids
+  //const int fhcalid = fhcal_modules->GetId();
+  const int ifhcalphi = out_config->GetBranchConfig(str_fhcal_branch).GetFieldId("phi");
+  const int ifhcalsign = out_config->GetBranchConfig(str_fhcal_branch).GetFieldId("signal_eta_signed");
+
+  // tpc_tracks' additional field ids
+  const int inhits = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nhits");
+  const int inhits_poss = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("nhits_poss");
+  const int icharge = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("charge");
+  const int idcax = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_x");
+  const int idcay = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_y");
+  const int idcaz = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dca_z");
+  const int ichi2 = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("chi2");
+  const int itof_mass2 = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("tof_mass2");
+  const int itof_flag = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("tof_flag");
+  const int idedx = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("dedx");
+  const int ipid_prob_pion = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("pid_prob_pion");
+  const int ipid_prob_kaon = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("pid_prob_kaon");
+  const int ipid_prob_proton = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("pid_prob_proton");
+  const int imc_E = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_E");
+  const int imc_pT = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_pT");
+  const int imc_eta = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_eta");
+  const int imc_phi = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_phi");
+  const int imc_rapidity = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_rapidity");
+  const int imc_mother_id = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_mother_id");
+  const int imc_pdg = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("mc_pdg");
+  const int ietasign = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("eta_sign");
+  const int irapidity_pion = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("rapidity_pion");
+  const int irapidity_pdg = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("rapidity_pdg");
+  const int irapidity_kaon = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("rapidity_kaon");
+  const int irapidity_proton = out_config->GetBranchConfig(str_tpc_tracks_branch).GetFieldId("rapidity_proton");
+
+  // mc_tracks' additional field ids
+  const int imother_id = out_config->GetBranchConfig(str_mc_tracks_branch).GetFieldId("mother_id");
+  const int icharge_mc = out_config->GetBranchConfig(str_mc_tracks_branch).GetFieldId("is_charged");
+  const int ietasign_mc = out_config->GetBranchConfig(str_mc_tracks_branch).GetFieldId("eta_sign");
+
+  // Create branches in the output tree
+  outTree->Branch(str_reco_event_branch.c_str(), "AnalysisTree::EventHeader", &reco_event, 32000, 99);
+  outTree->Branch(str_mc_event_branch.c_str(), "AnalysisTree::EventHeader", &mc_event, 32000, 99);
+  outTree->Branch(str_tpc_tracks_branch.c_str(), "AnalysisTree::TrackDetector", &tpc_tracks, 256000, 99);
+  outTree->Branch(str_fhcal_branch.c_str(), "AnalysisTree::ModuleDetector",  &fhcal_modules, 128000, 99);
+  outTree->Branch(str_mc_tracks_branch.c_str(), "AnalysisTree::Particles",  &mc_tracks, 256000, 99);
+  outTree->Branch(str_tpc2mc_tracks_branch.c_str(), "AnalysisTree::Matching", &tpc2mc_tracks, 32000, 99);
+
+  // Printout basic configuration info
+  out_config->Print();
+
+  // Starting event loop
+  TVector3 primaryVertex, mc_assoc_mom;
+  std::set <Int_t> UsedMCTracks;        // using to remap mc-reco track matching
+  std::map <Int_t,Int_t> InitMcNewMcId; // map[old-mc-id] = new-mc-id
+  Float_t FHCalSumEnergy[Num_Of_Modules];
+  Int_t   FHCalNumOfHits[Num_Of_Modules];
+  Int_t n_entries = dstTree->GetEntriesFast();
+
+  bool isGoodPID;
+  int charge;
+
+  for (int iEv = 0; iEv < n_entries; iEv++)
+  {
+    std::cout << "Event [" << iEv << "/" << n_entries << "]" << std::endl;
+    dstTree->GetEntry(iEv);
+
+    UsedMCTracks.clear();
+    InitMcNewMcId.clear();
+    tpc2mc_tracks->Clear();
+    for (int i=0; i<Num_Of_Modules; i++)
+    {
+      FHCalSumEnergy[i] = 0.;
+      FHCalNumOfHits[i] = 0;
+    }
+
+    tpc_tracks->ClearChannels();
+    fhcal_modules->ClearChannels();
+    mc_tracks->ClearChannels();
+
+    // Reading Reco Event
+    MpdVertex *vertex = (MpdVertex *)vertexes->First();
+    vertex->Position(primaryVertex);
+    reco_event->SetVertexPosition3(primaryVertex);
+
+    // Read MC Event
+    TVector3 vtx(MCHeader->GetX(),MCHeader->GetY(),MCHeader->GetZ());
+    mc_event->SetVertexPosition3(vtx);
+    mc_event->SetField(float(MCHeader->GetB()), iB);
+    mc_event->SetField(float(MCHeader->GetRotZ()), iPhiRp);
+
+    // Read energy in FHCal modules 
+   fhcal_modules->Reserve(Num_Of_Modules);
+    for (int imodule=0; imodule<Num_Of_Modules; imodule++)
+    {
+      auto *module = fhcal_modules->AddChannel();
+      module->Init(out_config->GetBranchConfig(str_fhcal_branch));
+      module->SetSignal(0.f);
+    }
+    Int_t number_of_FHCal_hits = FHCalHits->GetEntriesFast();
+    for(int ihit=0; ihit<number_of_FHCal_hits; ihit++)
+    {
+      FHCalHit = (MpdZdcDigi*) FHCalHits->UncheckedAt(ihit);
+      Int_t DetId = FHCalHit->GetDetectorID();
+      Int_t ModId = FHCalHit->GetModuleID()-1;
+      Int_t ModNumber = ModId + (Num_Of_Modules/2) * (DetId-1);
+
+      FHCalSumEnergy[ModNumber] += FHCalHit->GetELoss();
+      FHCalNumOfHits[ModNumber]++;
+    }
+    for (int imodule=0; imodule<Num_Of_Modules; imodule++)
+    {
+      auto& module = fhcal_modules->Channel(imodule);
+      int fhcal_sign = (imodule < 45) ? -1 : 1;
+      module.SetNumber(FHCalNumOfHits[imodule]); // Number of hits that got in the module
+      module.SetSignal(FHCalSumEnergy[imodule]); // Total energy from hits in the module
+      module.SetField(float(GetFHCalPhi(imodule)), ifhcalphi);
+      module.SetField(float(FHCalSumEnergy[imodule]*(float)fhcal_sign), ifhcalsign);
+    }
+
+    // Reading Reco Tracks
+    MpdGlobalTracks = (TClonesArray*)MPDEvent->GetGlobalTracks();
+    Int_t Num_of_tpc_tracks = MpdGlobalTracks->GetEntriesFast();
+    tpc_tracks->Reserve(Num_of_tpc_tracks);
+
+    for (int itrack=0; itrack<Num_of_tpc_tracks; itrack++)
+    {
+      MpdTrack* mpdtrack = (MpdTrack*) MpdGlobalTracks->UncheckedAt(itrack);
+      UsedMCTracks.insert(mpdtrack->GetID());
+      auto *track = tpc_tracks->AddChannel();
+      track->Init(out_config->GetBranchConfig(str_tpc_tracks_branch));
+
+      int tpc_eta_sign = (mpdtrack->GetEta() < 0) ? -1 : 1;
+
+      track->SetMomentum(mpdtrack->GetPx(),mpdtrack->GetPy(),mpdtrack->GetPz());
+      track->SetField(int(mpdtrack->GetNofHits()), inhits);
+      track->SetField(int(mpdtrack->GetNofHitsPossTpc()), inhits_poss);
+      track->SetField(float(mpdtrack->GetTofMass2()), itof_mass2);
+      track->SetField(float(mpdtrack->GetTofFlag()), itof_flag);
+      track->SetField(float(mpdtrack->GetdEdXTPC()), idedx);
+      track->SetField(float(mpdtrack->GetChi2()), ichi2);
+      track->SetField(float(mpdtrack->GetDCAX()), idcax);
+      track->SetField(float(mpdtrack->GetDCAY()), idcay);
+      track->SetField(float(mpdtrack->GetDCAZ()), idcaz);
+      charge = (mpdtrack->GetPt() < 0) ? 1 : -1;
+      track->SetField(int(charge), icharge);
+      track->SetField(int(tpc_eta_sign), ietasign);
+      Float_t mpdtrack_p = TMath::Sqrt(TMath::Power(mpdtrack->GetPx(),2) + TMath::Power(mpdtrack->GetPy(),2) + TMath::Power(mpdtrack->GetPz(),2));
+      track->SetField(float(GetRapidity(mpdtrack_p, mpdtrack->GetPz(), 0)), irapidity_pion);
+      track->SetField(float(GetRapidity(mpdtrack_p, mpdtrack->GetPz(), 1)), irapidity_kaon);
+      track->SetField(float(GetRapidity(mpdtrack_p, mpdtrack->GetPz(), 2)), irapidity_proton);
+
+      // PID
+      if (mpdtrack->GetTofFlag() == 2 || mpdtrack->GetTofFlag() == 6)
+      {
+        // TOF + TPC
+        isGoodPID = pid->FillProbs(TMath::Abs(mpdtrack->GetPt()) * TMath::CosH(mpdtrack->GetEta()),
+                                   mpdtrack->GetdEdXTPC(), mpdtrack->GetTofMass2(), charge);
+      }
+      else
+      {
+        // TPC only
+        isGoodPID = pid->FillProbs(TMath::Abs(mpdtrack->GetPt()) * TMath::CosH(mpdtrack->GetEta()),
+                                   mpdtrack->GetdEdXTPC(), charge);
+      }
+
+      if (isGoodPID)
+      {
+        track->SetField(float(pid->GetProbPi()), ipid_prob_pion);
+        track->SetField(float(pid->GetProbKa()), ipid_prob_kaon);
+        track->SetField(float(pid->GetProbPr()), ipid_prob_proton);
+      }
+      else
+      {
+        track->SetField(float(-999.), ipid_prob_pion);
+        track->SetField(float(-999.), ipid_prob_kaon);
+        track->SetField(float(-999.), ipid_prob_proton);
+      }
+
+#ifdef _MCSTACK_
+      FairMCTrack *mctrack = (FairMCTrack*) MCTracks->UncheckedAt(mpdtrack->GetID());
+#endif
+#ifdef _MPDMCSTACK_
+      MpdMCTrack *mctrack = (MpdMCTrack*) MCTracks->UncheckedAt(mpdtrack->GetID());
+#endif
+      mc_assoc_mom.SetXYZ(mctrack->GetPx(), mctrack->GetPy(), mctrack->GetPz());
+      track->SetField(float(mctrack->GetEnergy()), imc_E);
+      track->SetField(float(mc_assoc_mom.Pt()), imc_pT);
+      track->SetField(float(mc_assoc_mom.Eta()), imc_eta);
+      track->SetField(float(mc_assoc_mom.Phi()), imc_phi);
+      track->SetField(float(mctrack->GetRapidity()), imc_rapidity);
+      track->SetField(int(mctrack->GetMotherId()), imc_mother_id);
+      track->SetField(int(mctrack->GetPdgCode()), imc_pdg);
+      track->SetField(float(GetRapidityPDG(mpdtrack_p, mpdtrack->GetPz(), mctrack->GetPdgCode())), irapidity_pdg);
+
+    } // End of the tpc track loop
+
+    // Read Mc tracks
+    Int_t Num_of_mc_tracks = MCTracks->GetEntriesFast();
+    mc_tracks->Reserve(Num_of_mc_tracks);
+
+    for (int imctrack=0; imctrack<Num_of_mc_tracks; imctrack++)
+    {
+#ifdef _MCSTACK_
+      FairMCTrack *mctrack = (FairMCTrack*) MCTracks->UncheckedAt(imctrack);
+#endif
+#ifdef _MPDMCSTACK_
+      MpdMCTrack *mctrack = (MpdMCTrack*) MCTracks->UncheckedAt(imctrack);
+#endif
+      bool isUsed = (UsedMCTracks.count(imctrack));
+
+      // If motherId != 1 and mc track doesn't relate to any reco track - skip
+      if (mctrack->GetMotherId() != -1 && !isUsed) continue;
+
+      auto *track = mc_tracks->AddChannel();
+      track->Init(out_config->GetBranchConfig(str_mc_tracks_branch));
+
+      // Collect new Mc Ids
+      InitMcNewMcId[imctrack] = track->GetId();
+
+      int mc_eta_sign = (mctrack->GetRapidity() < 0) ? -1 : 1;
+
+      track->SetMomentum(mctrack->GetPx(), mctrack->GetPy(), mctrack->GetPz());
+      track->SetPid(int(mctrack->GetPdgCode()));
+      track->SetMass(float(mctrack->GetMass()));
+      track->SetField(int(mctrack->GetMotherId()), imother_id);
+      track->SetField(bool(IsCharged(mctrack->GetPdgCode())), icharge_mc);
+      track->SetField(int(mc_eta_sign), ietasign_mc);
+    } // End of the mc track loop
+
+    // reco-mc tracks matching
+    for (int itrack=0; itrack<Num_of_tpc_tracks; itrack++)
+    {
+      MpdTrack* mpdtrack = (MpdTrack*) MpdGlobalTracks->UncheckedAt(itrack);
+      const int tpc_track_id = itrack;
+      const int mc_track_id  = InitMcNewMcId[mpdtrack->GetID()];
+      tpc2mc_tracks->AddMatch(tpc_track_id, mc_track_id);
+    }
+
+    outTree->Fill();
+  
+  } // End of the event loop
+  
+  outFile->cd();
+  outTree->Print();
+  outTree->Write();
+  if (isDataHeaderDefined) dataHeader->Write("DataHeader");
+  out_config->Write("Configuration");
+  outFile->Close();
+
+  timer.Stop();
+  timer.Print();
+}
+
+
+float get_beamP(float sqrtSnn, float m_target, float m_beam)
+{
+  return sqrt( pow((pow(sqrtSnn,2) - pow(m_target,2) - pow(m_beam,2))/(2*m_target), 2) - pow(m_beam,2) );
+}
+
+float GetFHCalPhi(int iModule)
+{
+  const int Nmodules = 45;
+  int xAxisSwitch = (iModule < Nmodules) ? 1 : -1;
+  int module = (iModule < Nmodules) ? iModule : iModule - Nmodules;
+  float x, y, phi = -999.;
+  if (module >= 0 && module <= 4)
+  {
+    y = 45.;
+    x = (module - 2) * 15.;
+    phi = TMath::ATan2(y, x * xAxisSwitch);
+  }
+  else if ((module >= 5) && (module <= 39))
+  {
+    y = (3 - (module + 2) / 7) * 15.;
+    x = (3 - (module + 2) % 7) * 15.;
+    phi = TMath::ATan2(y, x * xAxisSwitch);
+  }
+  else if ((module >= 40) && (module <= 44))
+  {
+    y = -45.;
+    x = (module - 42) * 15.;
+    phi = TMath::ATan2(y, x * xAxisSwitch);
+  }
+  return phi;
+}
+
+TVector3 GetFHCalPos(int iModule)
+{
+  const int Nmodules = 45;
+  int xAxisSwitch = (iModule < Nmodules) ? 1 : -1;
+  int module = (iModule < Nmodules) ? iModule : iModule - Nmodules;
+  float x, y, z;
+  z = (iModule < Nmodules) ? -320. : 320.;
+  if (module >= 0 && module <= 4)
+  {
+    y = 45.;
+    x = (module - 2) * 15.;
+  }
+  else if ((module >= 5) && (module <= 39))
+  {
+    y = (3 - (module + 2) / 7) * 15.;
+    x = (3 - (module + 2) % 7) * 15.;
+  }
+  else if ((module >= 40) && (module <= 44))
+  {
+    y = -45.;
+    x = (module - 42) * 15.;
+  }
+  TVector3 vec(x * xAxisSwitch, y, z);
+
+  return vec;
+}
+
+Bool_t IsCharged(Int_t pdg)
+{
+  auto particle = (TParticlePDG *)TDatabasePDG::Instance()->GetParticle(pdg);
+  if (!particle)
+    return false;
+  return ( particle->Charge() != 0 );
+}
+
+Float_t GetRapidity(Float_t p, Float_t pz, Int_t pid_type)
+{
+  Float_t mass; // in GeV/c^2
+  if (pid_type == 0) mass = 0.13957; // pions
+  if (pid_type == 1) mass = 0.493677; // kaons
+  if (pid_type == 2) mass = 0.938272; // protons
+
+  if (pid_type != 0 && pid_type != 1 && pid_type != 2) return -999.;
+
+  Float_t E = TMath::Sqrt(p*p + mass*mass);
+
+  return 0.5 * TMath::Log((E + pz)/(E-pz));
+}
+
+Float_t GetRapidityPDG(Float_t p, Float_t pz, Int_t pdg)
+{
+  Float_t mass;
+
+  auto pdgParticle = (TParticlePDG*)TDatabasePDG::Instance()->GetParticle(pdg);
+  if (!pdgParticle) return -999.;
+  
+  mass = pdgParticle->Mass();
+  Float_t E = TMath::Sqrt(p*p + mass*mass);
+
+  return 0.5 * TMath::Log((E + pz)/(E - pz));
+}

File diff suppressed because it is too large
+ 1075 - 0
mpd-analysis-config.yml


+ 654 - 0
mpd-correlation-plain.yml

@@ -0,0 +1,654 @@
+_detectors: &detectors_mc
+  - name: mc_hadrons_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_hadrons_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_proton_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_proton_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_kaon_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_kaon_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_akaon_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_akaon_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_pion_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_pion_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_apion_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_apion_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mc_TPC_EP_L
+    tags: [ qn_vector ]
+    correction-step: plain
+  - name: mc_TPC_EP_R
+    tags: [ qn_vector ]
+    correction-step: plain
+  - name: mc_FHCal_L
+    tags: [ FHCal_qn_vector ]
+    correction-step: plain
+  - name: mc_FHCal_R
+    tags: [ FHCal_qn_vector ]
+    correction-step: plain
+  - name: mc_RP
+    correction-step: plain
+    tags: [ FHCal_qn_vector ]
+
+_detectors: &detectors_reco
+  - name: reco_hadrons_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_hadrons_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_proton_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_proton_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_kaon_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_kaon_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_akaon_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_akaon_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_pion_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_pion_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_apion_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_apion_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: reco_TPC_EP_L
+    tags: [ qn_vector ]
+    correction-step: plain
+  - name: reco_TPC_EP_R
+    tags: [ qn_vector ]
+    correction-step: plain
+  - name: reco_FHCal_L
+    tags: [ FHCal_qn_vector ]
+    correction-step: plain
+  - name: reco_FHCal_R
+    tags: [ FHCal_qn_vector ]
+    correction-step: plain
+  - name: reco_RP
+    correction-step: plain
+    tags: [ FHCal_qn_vector ]
+
+_detectors: &detectors_mcreco
+  - name: mcreco_hadrons_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_hadrons_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_proton_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_proton_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_kaon_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_kaon_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_akaon_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_akaon_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_pion_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_pion_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_apion_L
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_apion_R
+    tags: [ un_vector ]
+    correction-step: plain
+  - name: mcreco_TPC_EP_L
+    tags: [ qn_vector ]
+    correction-step: plain
+  - name: mcreco_TPC_EP_R
+    tags: [ qn_vector ]
+    correction-step: plain
+  - name: mcreco_FHCal_L
+    tags: [ FHCal_qn_vector ]
+    correction-step: plain
+  - name: mcreco_FHCal_R
+    tags: [ FHCal_qn_vector ]
+    correction-step: plain
+  - name: mcreco_RP
+    correction-step: plain
+    tags: [ FHCal_qn_vector ]
+
+_axes:
+  - &centrality
+    name: McEvent_B
+    bin-edges: [0., 4.06, 5.84, 7.17, 8.27, 9.26, 10.17, 11.04, 11.99, 14.]
+
+_components:
+  - &v2_sp_components
+    [ x2,y2 ]
+  - &v2_ep_components
+    [ cos2,sin2 ]
+  - &fhcal_sp_components
+    [ x1,y1 ]
+  - &fhcal_ep_components
+    [ cos1,sin1 ]
+  - &v1_sp_components
+    [ x1,y1 ]
+  - &v1_ep_components
+    [ cos1,sin1 ]
+
+_tasks_mc:
+  # u2 x Q2
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/EP"
+    axes: [ *centrality ]
+  # u2 x Q1 x Q1
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/EP"
+    axes: [ *centrality ]
+  # u1 x Q1
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mc
+        components: *v1_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v1_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v1/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mc
+        components: *v1_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v1_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v1/uQ/EP"
+    axes: [ *centrality ]
+  # Q2 x Q2
+  - args:
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v2/QQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mc
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v2/QQ/EP"
+    axes: [ *centrality ]
+  # Q1 x Q1
+  - args:
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v1/QQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mc
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v1/QQ/EP"
+    axes: [ *centrality ]
+
+_tasks_reco:
+  # u2 x Q2
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/EP"
+    axes: [ *centrality ]
+  # u2 x Q1 x Q1
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/EP"
+    axes: [ *centrality ]
+  # u1 x Q1
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_reco
+        components: *v1_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v1_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v1/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_reco
+        components: *v1_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v1_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v1/uQ/EP"
+    axes: [ *centrality ]
+  # Q2 x Q2
+  - args:
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v2/QQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_reco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v2/QQ/EP"
+    axes: [ *centrality ]
+  # Q1 x Q1
+  - args:
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v1/QQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_reco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v1/QQ/EP"
+    axes: [ *centrality ]
+
+_tasks_mcreco:
+  # u2 x Q2
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/EP"
+    axes: [ *centrality ]
+  # u2 x Q1 x Q1
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v2/uQ/EP"
+    axes: [ *centrality ]
+  # u1 x Q1
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v1_sp_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v1_sp_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v1/uQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ un_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v1_ep_components
+        correction-steps: [ plain ]
+        weight: sumw
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v1_ep_components
+        correction-steps: [ plain ]
+        weight: ones
+    n-samples: 50
+    weights-type: observable
+    folder: "/v1/uQ/EP"
+    axes: [ *centrality ]
+  # Q2 x Q2
+  - args:
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_sp_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v2/QQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *v2_ep_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v2/QQ/EP"
+    axes: [ *centrality ]
+  # Q1 x Q1
+  - args:
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_sp_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v1/QQ/SP"
+    axes: [ *centrality ]
+  - args:
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+      - query: { tags: { any-in: [ FHCal_qn_vector ] } }
+        query-list: *detectors_mcreco
+        components: *fhcal_ep_components
+        correction-steps: [ plain ]
+    n-samples: 50
+    weights-type: reference
+    folder: "/v1/QQ/EP"
+    axes: [ *centrality ]