mpd_summer2019/real-flow/MpdCalculator.C

//!!!!!!!!!!!!!!!!!!!!!!
#include "MpdCalculator.h"

MpdCalculator::MpdCalculator(TString inFileName, TString outFileName, TString dcaFileName)
{

  this->inFileName = inFileName;
  this->outFileName = outFileName;

  inChain = new TChain("cbmsim_reduced");
  inChain->Add(inFileName.Data());

  outFile = new TFile(outFileName.Data(), "RECREATE");

  dcaFile = new TFile(dcaFileName.Data(), "READ");

  inChain->SetBranchAddress("b_mc", &b_mc);
  inChain->SetBranchAddress("centrality_tpc_mpd", &centrality_tpc_mpd);
  inChain->SetBranchAddress("n_tracks_mc", &n_tracks_mc);
  inChain->SetBranchAddress("n_tracks_mpd", &n_tracks_mpd);
  inChain->SetBranchAddress("k_tracks_mpd", &k_tracks_mpd);
  inChain->SetBranchAddress("signed_pt_mpd", signed_pt_mpd);
  inChain->SetBranchAddress("eta_mpd", eta_mpd);
  inChain->SetBranchAddress("n_hits_mpd", n_hits_mpd);
  inChain->SetBranchAddress("phi_mpd", phi_mpd);
  inChain->SetBranchAddress("theta_mpd", theta_mpd);
  inChain->SetBranchAddress("tof_beta_mpd", tof_beta_mpd);
  inChain->SetBranchAddress("id_from_mc_mpd", id_from_mc_mpd);
  inChain->SetBranchAddress("mother_ID_mc", mother_ID_mc);
  inChain->SetBranchAddress("PDG_code_mc", PDG_code_mc);
  inChain->SetBranchAddress("phiEP_mc", &phiEP_mc);
  inChain->SetBranchAddress("ZDC_energy_mpd", ZDC_energy_mpd);
  inChain->SetBranchAddress("pt_mc", pt_mc);
  inChain->SetBranchAddress("chi2_mpd", chi2_mpd);
  inChain->SetBranchAddress("eta_mc", eta_mc);
  inChain->SetBranchAddress("px_mc", px_mc);
  inChain->SetBranchAddress("py_mc", py_mc);
  inChain->SetBranchAddress("pz_mc", pz_mc);
  inChain->SetBranchAddress("energy_mc", energy_mc);
  inChain->SetBranchAddress("DCA_x_mpd", DCA_x_mpd);
  inChain->SetBranchAddress("DCA_y_mpd", DCA_y_mpd);
  inChain->SetBranchAddress("DCA_z_mpd", DCA_z_mpd);
  inChain->SetBranchAddress("pid_tpc_prob_pion_mpd", pid_tpc_prob_pion_mpd);
  inChain->SetBranchAddress("pid_tpc_prob_kaon_mpd", pid_tpc_prob_kaon_mpd);
  inChain->SetBranchAddress("pid_tpc_prob_proton_mpd", pid_tpc_prob_proton_mpd);
  inChain->SetBranchAddress("TOF_flag_mpd", TOF_flag_mpd);

  h_nhits_TPC = new TH1F("h_nhits_TPC", "number of hits in TPC distribution;N_{hits};N_{tracks};", 60, 0., 60.);
  h_nhits_TPC_after = new TH1F("h_nhits_TPC_after", "number of hits in TPC distribution after cuts;N_{hits};N_{tracks};", 60, 0., 60.);
  h_multiplicity = new TH1F("h_multiplicity", "multiplicity in TPC;N_{tracks} in TPC;N_{events};", 1600, 0., 1600.);
  h_multiplicity_before = new TH1F("h_multiplicity_before", "multiplicity in TPCbefore cuts;N_{tracks} in TPC;N_{events};", 1600, 0., 1600.);
  h_impact_parameter = new TH1F("h_impact_parameter", "impact parameter distribution;b,fm;N_{events};", 100, 0., 20.);
  h2_energyZDC_L_vs_energy_ZDC_R = new TH2F("h2_energyZDC_L_vs_energy_ZDC_R", "FHCal energy left vs right;E_{left};E_{right};", 300, 0., 70., 300, 0., 70.);
  p_b_vs_multiplicity = new TProfile("p_b_vs_multiplicity", "Impact parameter vs TPC multiplicity", 800, 0, 1600);
  h2_b_vs_multiplicity = new TH2F("h2_b_vs_multiplicity", "Impact parameter vs TPC multiplicity", 800, 0., 1600., 100, 0., 20.);
  h2_b_vs_multiplicity->Sumw2();
  p_b_vs_multiplicity->Sumw2();
  p_b_vs_energy = new TProfile("p_b_vs_energy", "Impact parameter vs FHCal energy", 200, 0, 70);
  h2_b_vs_energy = new TH2F("h2_b_vs_energy", "Impact parameter vs FHCal energy", 200, 0., 70., 200, 0., 20.);
  h2_b_vs_energy->Sumw2();
  p_b_vs_energy->Sumw2();
  h2_b_vs_centrality = new TH2F("h2_b_vs_centrality", "h2_b_vs_centrality", 100, 0, 100, 200, 0., 20.);

  char name[200];
  char title[200];

  for (Int_t dim = 0; dim < Ndim; dim++)
  {
    for (Int_t eta_bin = 0; eta_bin < NetaBins; eta_bin++)
    {
      f_pt_fit[dim][eta_bin] = (TF1 *)dcaFile->Get(Form("sigma_pt_fit%i%i", dim, eta_bin));
    }
    for (Int_t pt_bin = 0; pt_bin < NptBins; pt_bin++)
    {
      for (Int_t eta_bin = 0; eta_bin < NetaBins; eta_bin++)
      {
        //f_dca[dim][pt_bin][eta_bin] = (TF1*) dcaFile->Get(Form("dca_fit[%i][%i][%i]",dim,pt_bin,eta_bin));

        sprintf(name, "h_DCA_primary[%i][%i][%i]", dim, pt_bin, eta_bin);
        sprintf(title, "DCA_{%i} primary for %.2f < p_{t} < %.2f and %.2f < #eta < %.2f ", dim, ptBins[pt_bin], ptBins[pt_bin + 1], etaBins[eta_bin], etaBins[eta_bin + 1]);
        h_DCA_primary[dim][pt_bin][eta_bin] = new TH1F(name, title, 200, -50., 50.);

        sprintf(name, "h_DCA_secondary[%i][%i][%i]", dim, pt_bin, eta_bin);
        sprintf(title, "DCA_{%i} secondary for %.2f < p_{t} < %.2f and %.2f < #eta < %.2f", dim, ptBins[pt_bin], ptBins[pt_bin + 1], etaBins[eta_bin], etaBins[eta_bin + 1]);
        h_DCA_secondary[dim][pt_bin][eta_bin] = new TH1F(name, title, 200, -50., 50.);

        sprintf(name, "h_DCA_all[%i][%i][%i]", dim, pt_bin, eta_bin);
        sprintf(title, "DCA_{%i} all for %.2f < p_{t} < %.2f and %.2f < #eta < %.2f", dim, ptBins[pt_bin], ptBins[pt_bin + 1], etaBins[eta_bin], etaBins[eta_bin + 1]);
        h_DCA_all[dim][pt_bin][eta_bin] = new TH1F(name, title, 200, -50., 50.);
      }
    }
  }

  for (int sort = 0; sort < _N_SORTS; ++sort)
  {
    sprintf(name, "p_momenta_resolution[%i]", sort);
    sprintf(title, "#Delta p_{t} / p_{t} for %s", sorts_of_particles[sort].Data());
    p_momenta_resolution[sort] = new TProfile(name, title, NptBins, ptBins);
  }

  for (int harm = 0; harm < _N_HARM; ++harm)
  {
    sprintf(name, "p_flow_wrt_RP_vs_centrality[%i]", harm);
    sprintf(title, "v_{%i} wrt RP for ;centrailty, %;v_{%i}", harm + 1, harm + 1);
    p_flow_wrt_RP_vs_centrality[harm] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

    for (int cenralityBin = 0; cenralityBin < NcentralityBinsFlow; ++cenralityBin)
    {
      sprintf(name, "p_flow_wrt_RP_vs_pt[%i][%i]", cenralityBin, harm);
      sprintf(title, "v_{%i} wrt RP for %.2f < b < %.2f;p_{t};", harm + 1, centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
      p_flow_wrt_RP_vs_pt[cenralityBin][harm] = new TProfile(name, title, NptBins, ptBins);

      sprintf(name, "p_flow_wrt_RP_vs_eta[%i][%i]", cenralityBin, harm);
      sprintf(title, "v_{%i} wrt RP for %.2f < b < %.2f;#eta;", harm + 1, centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
      p_flow_wrt_RP_vs_eta[cenralityBin][harm] = new TProfile(name, title, NetaBins, etaBins);

      sprintf(name, "p_flow_wrt_RP_vs_rapidity[%i][%i]", cenralityBin, harm);
      sprintf(title, "v_{%i} wrt RP for %.2f < b < %.2f;#y;", harm + 1, centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
      p_flow_wrt_RP_vs_rapidity[cenralityBin][harm] = new TProfile(name, title, NrapidityBins, rapidityBins);
    }

    for (int _harm = 0; _harm < _N_HARM; ++_harm)
    {
      for (int method = 0; method < _N_METHOD; ++method)
      {
        sprintf(name, "p_true_Res_vs_b[%i][%i][%i]", harm, _harm, method);
        sprintf(title, "%s%i%s%i%s%s%s", "<cos", harm + 1, "(#Psi_{", _harm + 1, "}^{FULL} - #Psi_{RP})> using ",
                methods_names[method].Data(), ";b,fm;");
        p_true_Res_vs_b[harm][_harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

        sprintf(name, "%s%i%s%i%s%i%s", "p_Res2Psi_vs_b[", harm, "][", _harm, "][", method, "]");
        sprintf(title, "%s%i%s%i%s%s%s%i%s%s%s", "cos(", harm + 1, "(#Psi_{", _harm + 1, ",", methods_names[method].Data(), "}^{R} - #Psi_{",
                _harm + 1, ",", methods_names[method].Data(), "}^{L});b,fm;");
        p_Res2Psi_vs_b[harm][_harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

        sprintf(name, "p_flow_wrt_full_vs_centrality[%i][%i][%i]", harm, _harm, method);
        sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} ;centrality, %; ", harm + 1, _harm + 1, methods_names[method].Data());
        p_flow_wrt_full_vs_centrality[harm][_harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

        sprintf(name, "p_flow_wrt_full_vs_centrality_divided[%i][%i][%i]", harm, _harm, method);
        sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} divided ;centrality, %; ", harm + 1, _harm + 1, methods_names[method].Data());
        p_flow_wrt_full_vs_centrality_divided[harm][_harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

        for (int cenralityBin = 0; cenralityBin < NcentralityBinsFlow; ++cenralityBin)
        {
          sprintf(name, "p_flow_wrt_full_vs_pt[%i][%i][%i][%i]", cenralityBin, harm, _harm, method);
          sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} for %.2f < b < %.2f;p_{t};", harm + 1, _harm + 1, methods_names[method].Data(),
                  centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
          p_flow_wrt_full_vs_pt[cenralityBin][harm][_harm][method] = new TProfile(name, title, NptBins, ptBins);

          sprintf(name, "p_flow_wrt_full_vs_pt_divided[%i][%i][%i][%i]", cenralityBin, harm, _harm, method);
          sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} divided for %.2f < b < %.2f;p_{t};", harm + 1, _harm + 1, methods_names[method].Data(), centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
          p_flow_wrt_full_vs_pt_divided[cenralityBin][harm][_harm][method] = new TProfile(name, title, NptBins, ptBins);

          sprintf(name, "p_flow_wrt_full_vs_eta[%i][%i][%i][%i]", cenralityBin, harm, _harm, method);
          sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} for %.2f < b < %.2f;#eta;", harm + 1, _harm + 1, methods_names[method].Data(),
                  centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
          p_flow_wrt_full_vs_eta[cenralityBin][harm][_harm][method] = new TProfile(name, title, NetaBins, etaBins);

          sprintf(name, "p_flow_wrt_full_vs_eta_divided[%i][%i][%i][%i]", cenralityBin, harm, _harm, method);
          sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} divided for %.2f < b < %.2f;#eta;", harm + 1, _harm + 1, methods_names[method].Data(),
                  centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
          p_flow_wrt_full_vs_eta_divided[cenralityBin][harm][_harm][method] = new TProfile(name, title, NetaBins, etaBins);

          sprintf(name, "p_flow_wrt_full_vs_rapidity[%i][%i][%i][%i]", cenralityBin, harm, _harm, method);
          sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} for %.2f < b < %.2f;#y;", harm + 1, _harm + 1, methods_names[method].Data(),
                  centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
          p_flow_wrt_full_vs_rapidity[cenralityBin][harm][_harm][method] = new TProfile(name, title, NrapidityBins, rapidityBins);

          sprintf(name, "p_flow_wrt_full_vs_rapidity_divided[%i][%i][%i][%i]", cenralityBin, harm, _harm, method);
          sprintf(title, "v_{%i} wrt #Psi_{%i,%s}^{FULL} divided for %.2f < b < %.2f;y;", harm + 1, _harm + 1, methods_names[method].Data(),
                  centralityBinsFlow[cenralityBin], centralityBinsFlow[cenralityBin + 1]);
          p_flow_wrt_full_vs_rapidity_divided[cenralityBin][harm][_harm][method] = new TProfile(name, title, NrapidityBins, rapidityBins);
        }
      }
    }
  }

  for (int arm = 0; arm < _N_ARM; ++arm)
  {

    sprintf(name, "%s%i%s", "h_energy_ZDC_total[", arm, "]");
    sprintf(title, "%s%s%s", "FHCal energy distribution for ", arm_names[arm].Data(), " arm;E_{total};N_{events};");
    h_energy_ZDC_total[arm] = new TH1F(name, title, 100, 0., 22.);

    sprintf(name, "%s%i%s", "h2_energy_ZDC_vs_multiplicity[", arm, "]");
    sprintf(title, "%s%s%s", "Total FHCal energy in ", arm_names[arm].Data(), " arm vs TPC multiplicity;E_{total};multiplicity;");
    h2_energy_ZDC_vs_multiplicity[arm] = new TH2F(name, title, 300, 0., 70., 200, 0., 800.);

    for (int harm = 0; harm < _N_HARM; ++harm)
    {
      for (int method = 0; method < _N_METHOD; ++method)
      {
        sprintf(name, "%s%i%s%i%s%i%s", "p_qx_vs_b[", arm, "][", harm, "][", method, "]");
        sprintf(title, "%s%i%s%s%s%s", "Q_{x}^{", harm + 1, ",", arm_names[arm].Data(), "} using ", methods_names[method].Data());
        p_qx_vs_b[arm][harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

        sprintf(name, "%s%i%s%i%s%i%s", "p_qy_vs_b[", arm, "][", harm, "][", method, "]");
        sprintf(title, "%s%i%s%s%s%s", "Q_{y}^{", harm + 1, ",", arm_names[arm].Data(), "} using ", methods_names[method].Data());
        p_qy_vs_b[arm][harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);

        for (int _harm = 0; _harm < _N_HARM; ++_harm)
        {
          sprintf(name, "p_true_Res_half_vs_b[%i][%i][%i][%i]", arm, harm, _harm, method);
          sprintf(title, "<cos(%i(#Psi_{%i,%s}^{%s} - #Psi_{RP} ))> ;b,fm;", harm + 1, _harm + 1, methods_names[method].Data(), arm_names[arm].Data());
          p_true_Res_half_vs_b[arm][harm][_harm][method] = new TProfile(name, title, NcentralityBinsRes, centralityBinsRes);
        }
      }
    }
  }

  for (int centralityBin = 0; centralityBin < NcentralityBinsRes; ++centralityBin)
  {
    for (int sort = 0; sort < _N_SORTS; ++sort)
    {
      sprintf(name, "%s%i%s%i%s", "h_pt[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "P_{t} distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};N_{tracks};");
      h_pt[centralityBin][sort] = new TH1F(name, title, 100, 0., 3.5);

      sprintf(name, "%s%i%s%i%s", "h_eta[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#eta distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#eta;N_{tracks};");
      h_eta[centralityBin][sort] = new TH1F(name, title, 100, -2., 2.);

      sprintf(name, "%s%i%s%i%s", "h_phi[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#phi distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#phi;N_{tracks};");
      h_phi[centralityBin][sort] = new TH1F(name, title, 100, -4., 4.);

      sprintf(name, "%s%i%s%i%s", "h2_pt_vs_eta[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "p_{t} vs #eta distribution for ", sorts_of_particles[sort].Data(), " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};#eta;");
      h2_pt_vs_eta[centralityBin][sort] = new TH2F(name, title, 100, 0., 3.5, 100., -1.5, 1.5);

      sprintf(name, "%s%i%s%i%s", "h2_pt_vs_phi[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "p_{t} vs #phi distribution for ", sorts_of_particles[sort].Data(), " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};#phi;");
      h2_pt_vs_phi[centralityBin][sort] = new TH2F(name, title, 100, 0., 3.5, 100., -4., 4.);

      sprintf(name, "%s%i%s%i%s", "h2_phi_vs_eta[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#phi vs #eta distribution for ", sorts_of_particles[sort].Data(), " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#phi;#eta;");
      h2_phi_vs_eta[centralityBin][sort] = new TH2F(name, title, 100, -4., 4., 100., -1.5, 1.5); //

      sprintf(name, "%s%i%s%i%s", "h_pt_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "P_{t} distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};N_{tracks};");
      h_pt_after[centralityBin][sort] = new TH1F(name, title, 100, 0., 3.5);

      sprintf(name, "%s%i%s%i%s", "h_eta_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#eta distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#eta;N_{tracks};");
      h_eta_after[centralityBin][sort] = new TH1F(name, title, 100, -2., 2.);

      sprintf(name, "%s%i%s%i%s", "h_phi_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#phi distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#phi;N_{tracks};");
      h_phi_after[centralityBin][sort] = new TH1F(name, title, 100, -4., 4.);

      sprintf(name, "%s%i%s%i%s", "h2_pt_vs_eta_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "p_{t} vs #eta distribution for ", sorts_of_particles[sort].Data(), " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};#eta;");
      h2_pt_vs_eta_after[centralityBin][sort] = new TH2F(name, title, 100, 0., 3.5, 100., -1.5, 1.5);

      sprintf(name, "%s%i%s%i%s", "h2_pt_vs_phi_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "p_{t} vs #phi distribution for ", sorts_of_particles[sort].Data(), " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};#phi;");
      h2_pt_vs_phi_after[centralityBin][sort] = new TH2F(name, title, 100, 0., 3.5, 100., -4., 4.);

      sprintf(name, "%s%i%s%i%s", "h2_phi_vs_eta_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#phi vs #eta distribution for ", sorts_of_particles[sort].Data(), " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#phi;#eta;");
      h2_phi_vs_eta_after[centralityBin][sort] = new TH2F(name, title, 100, -4., 4., 100., -1.5, 1.5); //

      sprintf(name, "%s%i%s%i%s", "h_pt_mc[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "P_{t} distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};N_{tracks};");
      h_pt_mc[centralityBin][sort] = new TH1F(name, title, 100, 0., 3.5);

      sprintf(name, "%s%i%s%i%s", "h_eta_mc[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#eta distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#eta;N_{tracks};");
      h_eta_mc[centralityBin][sort] = new TH1F(name, title, 100, -2., 2.);

      sprintf(name, "%s%i%s%i%s", "h_phi_mc[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#phi distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#phi;N_{tracks};");
      h_phi_mc[centralityBin][sort] = new TH1F(name, title, 100, -4., 4.);

      sprintf(name, "%s%i%s%i%s", "h_pt_mc_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "P_{t} distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";p_{t};N_{tracks};");
      h_pt_mc_after[centralityBin][sort] = new TH1F(name, title, 100, 0., 3.5);

      sprintf(name, "%s%i%s%i%s", "h_eta_mc_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#eta distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#eta;N_{tracks};");
      h_eta_mc_after[centralityBin][sort] = new TH1F(name, title, 100, -2., 2.);

      sprintf(name, "%s%i%s%i%s", "h_phi_mc_after[", centralityBin, "][", sort, "]");
      sprintf(title, "%s%s%s%.2f%s%.2f%s", "#phi distribution for ", sorts_of_particles[sort].Data(),
              " for ", centralityBinsRes[centralityBin], " < b < ", centralityBinsRes[centralityBin + 1], ";#phi;N_{tracks};");
      h_phi_mc_after[centralityBin][sort] = new TH1F(name, title, 100, -4., 4.);
    }
  }
}

void MpdCalculator::FillTPC(Int_t cenrality_bin, EPParticle *event_plane_buffer, Int_t ep_particle_count, Float_t weightB)
{

  Double_t qx, qy;
  Double_t psi_N_HALF;
  phiEP_mc = ATan2(Sin(phiEP_mc), Cos(phiEP_mc));

  for (Int_t harm = 0; harm < _N_HARM; ++harm)
  {
    Double_t psi_N_R = GetPsiHalfTpc(event_plane_buffer, ep_particle_count, 1, harm + 1, qx, qy);
    Double_t psi_N_L = GetPsiHalfTpc(event_plane_buffer, ep_particle_count, -1, harm + 1, qx, qy);

    Double_t psi_N_FULL = GetPsiFullTpc(event_plane_buffer, ep_particle_count, harm + 1);
    for (Int_t _harm = 0; _harm < _N_HARM; ++_harm)
    {
      Double_t psi_N_R = GetPsiHalfTpc(event_plane_buffer, ep_particle_count, 1, _harm + 1, qx, qy);
      Double_t psi_N_L = GetPsiHalfTpc(event_plane_buffer, ep_particle_count, -1, _harm + 1, qx, qy);
      // << "Filling with centrality = " << centralityBinsRes[cenrality_bin]+0.1 << "Res2 = " << Cos((harm+1)*(psi_N_R - psi_N_L)) << endl;
      p_Res2Psi_vs_b[harm][_harm][0]->Fill(centralityBinsRes[cenrality_bin] + 0.1, Cos((harm + 1) * (psi_N_R - psi_N_L)), weightB);

      Double_t psi_N_FULL = GetPsiFullTpc(event_plane_buffer, ep_particle_count, harm + 1);
      p_true_Res_vs_b[harm][_harm][0]->Fill(centralityBinsRes[cenrality_bin] + 0.1, Cos((harm + 1) * (psi_N_FULL - phiEP_mc)), weightB);

      for (Int_t arm = 0; arm < _N_ARM; ++arm)
      {
        Double_t psi_N_HALF = GetPsiHalfTpc(event_plane_buffer, ep_particle_count, -2 * arm + 1, _harm + 1, qx, qy);
        p_true_Res_half_vs_b[arm][harm][_harm][0]->Fill(centralityBinsRes[cenrality_bin] + 0.1, Cos((harm + 1) * (psi_N_HALF - phiEP_mc)), weightB);
      }
    }

    for (Int_t arm = 0; arm < _N_ARM; ++arm)
    {
      psi_N_HALF = GetPsiHalfTpc(event_plane_buffer, ep_particle_count, -2 * arm + 1, harm + 1, qx, qy);
      p_qx_vs_b[arm][harm][0]->Fill(centralityBinsRes[cenrality_bin] + 0.1, qx, weightB);
      p_qy_vs_b[arm][harm][0]->Fill(centralityBinsRes[cenrality_bin] + 0.1, qy, weightB);
    }
  }
}

void MpdCalculator::FillZDC(Int_t cenrality_bin, Float_t *ZDC_energy_mpd, Float_t weightB)
{
  Float_t values[10], absvalues[10];
  Double_t qx, qy;
  Double_t psi_N_HALF;
  phiEP_mc = ATan2(Sin(phiEP_mc), Cos(phiEP_mc));

  for (Int_t harm = 0; harm < _N_HARM; ++harm)
  {
    for (Int_t _harm = 0; _harm < _N_HARM; ++_harm)
    {
      Double_t psi_N_R = GetPsiHalfZdc(ZDC_energy_mpd, 0, _harm + 1, qx, qy);
      Double_t psi_N_L = GetPsiHalfZdc(ZDC_energy_mpd, 1, _harm + 1, qx, qy);
      //cout << "Filling with centrality = " << centralityBinsRes[cenrality_bin]+0.1 << "Res2 = " << Cos((harm+1)*(psi_N_R - psi_N_L)) << endl;
      p_Res2Psi_vs_b[harm][_harm][1]->Fill(centralityBinsRes[cenrality_bin] + 0.1, Cos((harm + 1) * (psi_N_R - psi_N_L)),weightB);

      Double_t psi_N_FULL = GetPsiFullZdc(ZDC_energy_mpd, _harm + 1);
      p_true_Res_vs_b[harm][_harm][1]->Fill(centralityBinsRes[cenrality_bin] + 0.1, Cos((harm + 1) * (psi_N_FULL - phiEP_mc)),weightB);

      for (Int_t arm = 0; arm < _N_ARM; ++arm)
      {
        Double_t psi_N_HALF = GetPsiHalfZdc(ZDC_energy_mpd, arm, _harm + 1, qx, qy);
        p_true_Res_half_vs_b[arm][harm][_harm][1]->Fill(centralityBinsRes[cenrality_bin] + 0.1, Cos((harm + 1) * (psi_N_HALF - phiEP_mc)),weightB);
      }
    }

    for (Int_t arm = 0; arm < _N_ARM; ++arm)
    {
      psi_N_HALF = GetPsiHalfZdc(ZDC_energy_mpd, arm, harm + 1, qx, qy);
      p_qx_vs_b[arm][harm][1]->Fill(centralityBinsRes[cenrality_bin] + 0.1, qx,weightB);
      p_qy_vs_b[arm][harm][1]->Fill(centralityBinsRes[cenrality_bin], qy, weightB);
    }
  }
}

void MpdCalculator::CalculateResolutions(Int_t nevents)
{
  const Float_t dB = 0.25;
  Float_t				B,wB;
  EPParticle *event_plane_buffer = new EPParticle[10000];
  if (nevents == 0)
    nevents = inChain->GetEntries();
  for (Int_t event = 0; event < nevents; ++event)
  {
    inChain->GetEntry(event); //reading all the branches for current event
    cout << "EVENT # " << event << endl;
    B = b_mc;
    wB = 2*TMath::Pi()*B*dB;
    centrality_tpc_mpd = 100 - centrality_tpc_mpd; //FIXED CENTRALITY!!!!!!!!!!!!!!!!!!!!!!

    h_impact_parameter->Fill(b_mc,wB); //filling impact parameter hitsogram with no cuts on events
    //Int_t multiplicity = GetMultiplicityTPC();
    if (n_tracks_mpd == 0)
      continue;																 //skipping empty events
    h_multiplicity_before->Fill(n_tracks_mpd,wB); //filling multiplicity TPC before cuts hitsogram
    p_b_vs_multiplicity->Fill(k_tracks_mpd, b_mc,wB);
    h2_b_vs_centrality->Fill(centrality_tpc_mpd, b_mc,wB);
    h2_b_vs_multiplicity->Fill(k_tracks_mpd, b_mc,wB);

    int cenrality_bin_res = GetCentralityBinRes(centrality_tpc_mpd); //getting the resolution bin in whic the event is
    if (cenrality_bin_res == -1)
      continue; // just in case...

    for (Long64_t track = 0; track < n_tracks_mc; ++track)
    {
      //~ _      _  _                                              _             __                                 _
      //~ | |    (_)| |                                            | |           / _|                               | |
      //~ | |__   _ | |_  ___    __ _  _ __  __ _  _ __ ___   ___  | |__    ___ | |_  ___   _ __  ___    ___  _   _ | |_  ___
      //~ | '_ \ | || __|/ _ \  / _` || '__|/ _` || '_ ` _ \ / __| | '_ \  / _ \|  _|/ _ \ | '__|/ _ \  / __|| | | || __|/ __|
      //~ | | | || || |_| (_) || (_| || |  | (_| || | | | | |\__ \ | |_) ||  __/| | | (_) || |  |  __/ | (__ | |_| || |_ \__ \
      //~ |_| |_||_| \__|\___/  \__, ||_|   \__,_||_| |_| |_||___/ |_.__/  \___||_|  \___/ |_|   \___|  \___| \__,_| \__||___/
      //~ __/ |
      //~ |___/
      int sort = 0;
      h_phi_mc[cenrality_bin_res][0]->Fill(ATan2(px_mc[track], py_mc[track]),wB);
      h_eta_mc[cenrality_bin_res][0]->Fill(eta_mc[track],wB);
      h_pt_mc[cenrality_bin_res][0]->Fill(pt_mc[track],wB);

      switch (PDG_code_mc[track])
      {
      case 211:
        sort = 1;
        break;
      case 2212:
        sort = 2;
        break;
      case 321:
        sort = 3;
        break;
      default:
        break;
      }

      if ((sort == 1) || (sort == 2) || (sort == 3))
      {
        h_phi_mc[cenrality_bin_res][sort]->Fill(ATan2(px_mc[track], py_mc[track]),wB);
        h_eta_mc[cenrality_bin_res][sort]->Fill(eta_mc[track],wB);
        h_pt_mc[cenrality_bin_res][sort]->Fill(pt_mc[track],wB);
      }

      //~ | |
      //~ ___  _   _ | |_  ___
      //~ / __|| | | || __|/ __|
      //~ | (__ | |_| || |_ \__ \
      //~ \___| \__,_| \__||___/

      Int_t pt_bin = GetPtBin(pt_mc[track]);
      Int_t eta_bin = GetEtaBin(eta_mc[track]);
      if ((eta_bin == -1) || (pt_bin == -1))
        continue;
      if (mother_ID_mc[track] > -1)
        continue;

      //~ | |    (_)     | |                                                    / _|| |                           | |
      //~ | |__   _  ___ | |_  ___    __ _  _ __  __ _  _ __ ___   ___    __ _ | |_ | |_  ___  _ __    ___  _   _ | |_  ___
      //~ | '_ \ | |/ __|| __|/ _ \  / _` || '__|/ _` || '_ ` _ \ / __|  / _` ||  _|| __|/ _ \| '__|  / __|| | | || __|/ __|
      //~ | | | || |\__ \| |_| (_) || (_| || |  | (_| || | | | | |\__ \ | (_| || |  | |_|  __/| |    | (__ | |_| || |_ \__ \
      //~ |_| |_||_||___/ \__|\___/  \__, ||_|   \__,_||_| |_| |_||___/  \__,_||_|   \__|\___||_|     \___| \__,_| \__||___/
      //~ __/ |
      //~ |___/
      sort = 0;
      h_phi_mc_after[cenrality_bin_res][0]->Fill(ATan2(px_mc[track], py_mc[track]),wB);
      h_eta_mc_after[cenrality_bin_res][0]->Fill(eta_mc[track],wB);
      h_pt_mc_after[cenrality_bin_res][0]->Fill(pt_mc[track],wB);

      switch (PDG_code_mc[track])
      {
      case 211:
        sort = 1;
        break;
      case 2212:
        sort = 2;
        break;
      case 321:
        sort = 3;
        break;
      default:
        break;
      }

      if ((sort == 1) || (sort == 2) || (sort == 3))
      {
        h_phi_mc_after[cenrality_bin_res][sort]->Fill(ATan2(px_mc[track], py_mc[track]),wB);
        h_eta_mc_after[cenrality_bin_res][sort]->Fill(eta_mc[track],wB);
        h_pt_mc_after[cenrality_bin_res][sort]->Fill(pt_mc[track],wB);
      }
    }

    int ep_particle_count = 0, multiplicity_tpc_after = 0;
    for (Long64_t track = 0; track < n_tracks_mpd; ++track) //loop over mpdtracks
    {
      //~ _      _  _                                              _             __                                 _
      //~ | |    (_)| |                                            | |           / _|                               | |
      //~ | |__   _ | |_  ___    __ _  _ __  __ _  _ __ ___   ___  | |__    ___ | |_  ___   _ __  ___    ___  _   _ | |_  ___
      //~ | '_ \ | || __|/ _ \  / _` || '__|/ _` || '_ ` _ \ / __| | '_ \  / _ \|  _|/ _ \ | '__|/ _ \  / __|| | | || __|/ __|
      //~ | | | || || |_| (_) || (_| || |  | (_| || | | | | |\__ \ | |_) ||  __/| | | (_) || |  |  __/ | (__ | |_| || |_ \__ \
      //~ |_| |_||_| \__|\___/  \__, ||_|   \__,_||_| |_| |_||___/ |_.__/  \___||_|  \___/ |_|   \___|  \___| \__,_| \__||___/
      //~ __/ |
      //~ |___/
      Int_t pt_bin = GetPtBin(Abs(signed_pt_mpd[track]));
      Int_t eta_bin = GetEtaBin(eta_mpd[track]);

      h_DCA_all[0][pt_bin][eta_bin]->Fill(DCA_x_mpd[track],wB);
      h_DCA_all[1][pt_bin][eta_bin]->Fill(DCA_y_mpd[track],wB);
      h_DCA_all[2][pt_bin][eta_bin]->Fill(DCA_z_mpd[track],wB);
      if (mother_ID_mc[id_from_mc_mpd[track]] == -1)
      {
        h_DCA_primary[0][pt_bin][eta_bin]->Fill(DCA_x_mpd[track],wB);
        h_DCA_primary[1][pt_bin][eta_bin]->Fill(DCA_y_mpd[track],wB);
        h_DCA_primary[2][pt_bin][eta_bin]->Fill(DCA_z_mpd[track],wB);
      }
      else if (mother_ID_mc[id_from_mc_mpd[track]] > -1)
      {
        h_DCA_secondary[0][pt_bin][eta_bin]->Fill(DCA_x_mpd[track],wB);
        h_DCA_secondary[1][pt_bin][eta_bin]->Fill(DCA_y_mpd[track],wB);
        h_DCA_secondary[2][pt_bin][eta_bin]->Fill(DCA_z_mpd[track],wB);
      }

      h_nhits_TPC->Fill(n_hits_mpd[track],wB); //number of hits in TPC per track, before eta, nhits and pt cuts, but after mother ID cut
      int sort = 0;
      h_phi[cenrality_bin_res][0]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
      h_eta[cenrality_bin_res][0]->Fill(eta_mpd[track],wB);
      h_pt[cenrality_bin_res][0]->Fill(Abs(signed_pt_mpd[track]),wB);
      h2_pt_vs_eta[cenrality_bin_res][0]->Fill(Abs(signed_pt_mpd[track]), eta_mpd[track],wB);
      h2_pt_vs_phi[cenrality_bin_res][0]->Fill(Abs(signed_pt_mpd[track]), ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
      h2_phi_vs_eta[cenrality_bin_res][0]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])), eta_mpd[track],wB);
      p_momenta_resolution[0]->Fill(Abs(signed_pt_mpd[track]), Abs(Abs(signed_pt_mpd[track]) - pt_mc[id_from_mc_mpd[track]]) / Abs(signed_pt_mpd[track]),wB);

      switch (PDG_code_mc[id_from_mc_mpd[track]])
      {
      case 211:
        sort = 1;
        break;
      case 2212:
        sort = 2;
        break;
      case 321:
        sort = 3;
        break;
      default:
        break;
      }

      /*if (pid_tpc_prob_pion_mpd[track] > 0.9 && 
        pid_tpc_prob_kaon_mpd[track] < 0.1 &&
        pid_tpc_prob_proton_mpd[track] < 0.1) sort = 1;

        if (pid_tpc_prob_pion_mpd[track] < 0.1 && 
        pid_tpc_prob_kaon_mpd[track] > 0.9 &&
        pid_tpc_prob_proton_mpd[track] < 0.1) sort = 3;

        if (pid_tpc_prob_pion_mpd[track] < 0.1 && 
        pid_tpc_prob_kaon_mpd[track] < 0.1 &&
        pid_tpc_prob_proton_mpd[track] > 0.9) sort = 2;

        if (signed_pt_mpd[track] >= 0 ) continue; // Only positively charged particles being selected
        if (TOF_flag_mpd[track] == 0) continue; // Only tracks with TOF hits being selected
        if (TOF_flag_mpd[track] == 4) continue; // Only tracks with TOF hits being selected*/

      if ((sort == 1) || (sort == 2) || (sort == 3))
      {
        h_phi[cenrality_bin_res][sort]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
        h_eta[cenrality_bin_res][sort]->Fill(eta_mpd[track],wB);
        h_pt[cenrality_bin_res][sort]->Fill(Abs(signed_pt_mpd[track]),wB);
        h2_pt_vs_eta[cenrality_bin_res][sort]->Fill(Abs(signed_pt_mpd[track]), eta_mpd[track],wB);
        h2_pt_vs_phi[cenrality_bin_res][sort]->Fill(Abs(signed_pt_mpd[track]), ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
        h2_phi_vs_eta[cenrality_bin_res][sort]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])), eta_mpd[track],wB);
        p_momenta_resolution[sort]->Fill(Abs(signed_pt_mpd[track]), Abs(Abs(signed_pt_mpd[track]) - pt_mc[id_from_mc_mpd[track]]) / Abs(signed_pt_mpd[track]),wB);
      }

      //~ | |
      //~ ___  _   _ | |_  ___
      //~ / __|| | | || __|/ __|
      //~ | (__ | |_| || |_ \__ \
      //~ \___| \__,_| \__||___/

      //if (id_from_mc_mpd[track] == -1) continue; //equivalent to mother id cut
      pt_bin = GetPtBin(Abs(signed_pt_mpd[track]));
      eta_bin = GetEtaBin(eta_mpd[track]);
      if ((eta_bin == -1) || (pt_bin == -1))
        continue;
      if (n_hits_mpd[track] < Cut_No_Of_hits_min)
        continue; //n hits in TPC per track cut

      //if (TMath::Abs(DCA_x_mpd[track]) >= f_dca[0][pt_bin][eta_bin]->GetParameter(2)*2) continue;
      //if (TMath::Abs(DCA_y_mpd[track]) >= f_dca[1][pt_bin][eta_bin]->GetParameter(2)*2) continue;
      //if (TMath::Abs(DCA_z_mpd[track]) >= f_dca[2][pt_bin][eta_bin]->GetParameter(2)*2) continue;
      TF1 sigma_fit_res_X = *f_pt_fit[0][eta_bin];
      TF1 sigma_fit_res_Y = *f_pt_fit[1][eta_bin];
      TF1 sigma_fit_res_Z = *f_pt_fit[2][eta_bin];

      if (TMath::Abs(DCA_x_mpd[track]) >= sigma_fit_res_X(Abs(signed_pt_mpd[track])) * 2)
        continue;
      if (TMath::Abs(DCA_y_mpd[track]) >= sigma_fit_res_Y(Abs(signed_pt_mpd[track])) * 2)
        continue;
      if (TMath::Abs(DCA_z_mpd[track]) >= sigma_fit_res_Z(Abs(signed_pt_mpd[track])) * 2)
        continue;
      //if (mother_ID_mc[id_from_mc_mpd[track]] > -1) continue;

      //~ | |    (_)     | |                                                    / _|| |                           | |
      //~ | |__   _  ___ | |_  ___    __ _  _ __  __ _  _ __ ___   ___    __ _ | |_ | |_  ___  _ __    ___  _   _ | |_  ___
      //~ | '_ \ | |/ __|| __|/ _ \  / _` || '__|/ _` || '_ ` _ \ / __|  / _` ||  _|| __|/ _ \| '__|  / __|| | | || __|/ __|
      //~ | | | || |\__ \| |_| (_) || (_| || |  | (_| || | | | | |\__ \ | (_| || |  | |_|  __/| |    | (__ | |_| || |_ \__ \
      //~ |_| |_||_||___/ \__|\___/  \__, ||_|   \__,_||_| |_| |_||___/  \__,_||_|   \__|\___||_|     \___| \__,_| \__||___/
      //~ __/ |
      //~ |___/

      h_nhits_TPC_after->Fill(n_hits_mpd[track],wB); //n hits in TPC after eta, nhits, pt and mother ID cuts
      multiplicity_tpc_after++;										//multiplicity in TPC after eta, nhits, pt and mother ID cuts
      sort = 0;
      //filling some distributions after track cuts, depending on sort of the particle from MC data
      h_phi_after[cenrality_bin_res][0]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
      h_eta_after[cenrality_bin_res][0]->Fill(eta_mpd[track],wB);
      h_pt_after[cenrality_bin_res][0]->Fill(Abs(signed_pt_mpd[track]),wB);
      h2_pt_vs_eta_after[cenrality_bin_res][0]->Fill(Abs(signed_pt_mpd[track]), eta_mpd[track],wB);
      h2_pt_vs_phi_after[cenrality_bin_res][0]->Fill(Abs(signed_pt_mpd[track]), ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
      h2_phi_vs_eta_after[cenrality_bin_res][0]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])), eta_mpd[track],wB);
      //~ p_momenta_resolution[0]->Fill(Abs(signed_pt_mpd[track]),Abs(Abs(signed_pt_mpd[track]) - pt_mc[id_from_mc_mpd[track]])
      //~ /Abs(signed_pt_mpd[track]));

      switch (PDG_code_mc[id_from_mc_mpd[track]])
      {
      case 211:
        sort = 1;
        break;
      case 2212:
        sort = 2;
        break;
      case 321:
        sort = 3;
        break;
      default:
        break;
      }

      /*if (pid_tpc_prob_pion_mpd[track] > 0.9 && 
        pid_tpc_prob_kaon_mpd[track] < 0.1 &&
        pid_tpc_prob_proton_mpd[track] < 0.1) sort = 1;

        if (pid_tpc_prob_pion_mpd[track] < 0.1 && 
        pid_tpc_prob_kaon_mpd[track] > 0.9 &&
        pid_tpc_prob_proton_mpd[track] < 0.1) sort = 3;

        if (pid_tpc_prob_pion_mpd[track] < 0.1 && 
        pid_tpc_prob_kaon_mpd[track] < 0.1 &&
        pid_tpc_prob_proton_mpd[track] > 0.9) sort = 2;

        if (signed_pt_mpd[track] >= 0 ) continue; // Only positively charged particles being selected
        if (TOF_flag_mpd[track] == 0) continue; // Only tracks with TOF hits being selected
        if (TOF_flag_mpd[track] == 4) continue; // Only tracks with TOF hits being selected*/

      if ((sort == 1) || (sort == 2) || (sort == 3))
      {
        h_phi_after[cenrality_bin_res][sort]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
        h_eta_after[cenrality_bin_res][sort]->Fill(eta_mpd[track],wB);
        h_pt_after[cenrality_bin_res][sort]->Fill(Abs(signed_pt_mpd[track]),wB);
        h2_pt_vs_eta_after[cenrality_bin_res][sort]->Fill(Abs(signed_pt_mpd[track]), eta_mpd[track],wB);
        h2_pt_vs_phi_after[cenrality_bin_res][sort]->Fill(Abs(signed_pt_mpd[track]), ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])),wB);
        h2_phi_vs_eta_after[cenrality_bin_res][sort]->Fill(ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track])), eta_mpd[track],wB);
        //~ p_momenta_resolution[sort]->Fill(Abs(signed_pt_mpd[track]),Abs(Abs(signed_pt_mpd[track]) - pt_mc[id_from_mc_mpd[track]])
        //~ /Abs(signed_pt_mpd[track]));
      }

      //~ if (Abs(eta_mpd[track]) > Cut_Eta_Min) //and finally, the resolution gap eta cut, then filling the buffer of EP particles:
      //~ {
      event_plane_buffer[ep_particle_count].Pt = Abs(signed_pt_mpd[track]);
      event_plane_buffer[ep_particle_count].Eta = eta_mpd[track];
      event_plane_buffer[ep_particle_count].Phi = ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track]));
      ep_particle_count++;
      //~ }
    }
    h_multiplicity->Fill(multiplicity_tpc_after,wB); //filling multiplicity TPC after cuts hitsogram, now that the loop over mpd tracks is over

    for (int arm = 0; arm < _N_ARM; ++arm)
    {
      h_energy_ZDC_total[arm]->Fill(GetTotalEnergy(ZDC_energy_mpd, arm),wB);																		 //filling the total ZDC energy for given event
      h2_energy_ZDC_vs_multiplicity[arm]->Fill(GetTotalEnergy(ZDC_energy_mpd, arm), multiplicity_tpc_after,wB); //filling total ZDC energy vs TPC
      //multiplicity after cuts for a given event
    }
    p_b_vs_energy->Fill(GetTotalEnergy(ZDC_energy_mpd, 0) + GetTotalEnergy(ZDC_energy_mpd, 1), b_mc,wB);
    h2_b_vs_energy->Fill(GetTotalEnergy(ZDC_energy_mpd, 0) + GetTotalEnergy(ZDC_energy_mpd, 1), b_mc,wB);
    h2_energyZDC_L_vs_energy_ZDC_R->Fill(GetTotalEnergy(ZDC_energy_mpd, 0), GetTotalEnergy(ZDC_energy_mpd, 1),wB); //filling ZDC energy L vs R for a given event

    int multiplicity_R = 0, multiplicity_L = 0;
    for (int ep_particle = 0; ep_particle < ep_particle_count; ++ep_particle) //TPC subevent multiplicity after cuts
    {
      if (event_plane_buffer[ep_particle].Eta > 0)
        multiplicity_R++;
      else if (event_plane_buffer[ep_particle].Eta < 0)
        multiplicity_L++;
    }

    //calculate the EP resolution using TPC if subevents have enough multipllicity
    if ((multiplicity_L >= 4) && (multiplicity_R >= 4))
      FillTPC(cenrality_bin_res, event_plane_buffer, ep_particle_count,wB);

    //if both ZDCs have some signal then calculate EP resolutions using ZDC
    if ((GetTotalEnergy(ZDC_energy_mpd, 0) != 0) && (GetTotalEnergy(ZDC_energy_mpd, 1) != 0))
      FillZDC(cenrality_bin_res, ZDC_energy_mpd,wB);
  }

  return;
}

//~ __  _
//~ / _|| |
//~ | |_ | |  ___ __      __
//~ |  _|| | / _ \\ \ /\ / /
//~ | |  | || (_) |\ V  V /
//~ |_|  |_| \___/  \_/\_/

void MpdCalculator::CalculateFlow(Int_t nevents, TString fitFile)
{
  TFile *f = new TFile(fitFile.Data());							//opening file with fit data
  TF1 *resolution_fit[_N_HARM][_N_HARM][_N_METHOD]; //declaring the fit-function array
  for (int harm = 0; harm < _N_HARM; ++harm)
  {
    for (int _harm = 0; _harm < _N_HARM; ++_harm)
    {
      for (int method = 0; method < _N_METHOD; ++method)
      {
        char name[200];
        sprintf(name, "resolution_fit[%i][%i][%i]", harm, _harm, method);
        resolution_fit[harm][_harm][method] = (TF1 *)f->Get(name); //reading fit-functions from file
      }
    }
  }

  const Float_t dB = 0.25;
  Float_t				B,wB;
  if (nevents == 0)
    nevents = inChain->GetEntries();
  for (Int_t event = 0; event < nevents; ++event) //loop over events
  {
    inChain->GetEntry(event); //reading all branches for a given event
    cout << "FLOW EVENT # " << event << endl;
    B = b_mc;
    wB = 2*TMath::Pi()*B*dB;
    //Int_t multiplicity = GetMultiplicityTPC();
    centrality_tpc_mpd = 100 - centrality_tpc_mpd;

    int cenrality_bin_flow = GetCentralityBinFlow(centrality_tpc_mpd); //getting b bin of a given event for flow measurements
    //if (cenrality_bin_flow == -1) continue;
    int cenrality_bin_res = GetCentralityBinRes(centrality_tpc_mpd); //getting the resolution bin in whic the event is
    if (cenrality_bin_res == -1)
      continue; // just in case...
    //skip the event if the ZDC signal is zero
    if ((GetTotalEnergy(ZDC_energy_mpd, 0) == 0) || (GetTotalEnergy(ZDC_energy_mpd, 1) == 0))
      continue;

    Double_t phi_EP = ATan2(Sin(phiEP_mc), Cos(phiEP_mc)); //unfold the generated event plane

    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //RECONSTRUCTED
    const Float_t p_mass2[] = {0.88035, 0.2437, 0.01948};
    const Int_t pdg_codes[] = {2212, 321, 211}; //0 - PROTON, 1 - KAON, 2 -PION
    Bool_t isPion, isKaon, isProton;

    //		______		 _		 _____
    //		| ___ \		| |		|  __ \
    //		| |__| |	| |		| |  \ |
    //		|	____/		| |		| |  | |
    //		| |				| |		| |__/ |
    //		|_|				|_|		|_____/
    Int_t sort = 0; // 0 -proton; 1 - kaon; 2 - pion;

    for (Long64_t track = 0; track < n_tracks_mpd; ++track) //loop over reconstructed tracks to fill flow buffer
    {
      isPion = kFALSE;
      isKaon = kFALSE;
      isProton = kFALSE;
      //if (id_from_mc_mpd[track] == -1) continue; //equivalent to mother ID cut
      if (n_hits_mpd[track] < Cut_No_Of_hits_min)
        continue; //n hits in TPC cut
      if (PDG_code_mc[id_from_mc_mpd[track]] != pdg_codes[sort])
        continue;

      Float_t Pt = Abs(signed_pt_mpd[track]);
      Float_t Eta = eta_mpd[track];
      Float_t Phi = ATan2(Sin(phi_mpd[track]), Cos(phi_mpd[track]));
      //~ Float_t Rapidity = 0.5*TMath::Log((energy_mc[id_from_mc_mpd[track]] + pz_mc[id_from_mc_mpd[track]])
      //~ /(energy_mc[id_from_mc_mpd[track]] - pz_mc[id_from_mc_mpd[track]]));

      Int_t pt_bin = GetPtBin(Pt);
      Int_t eta_bin = GetEtaBin(Eta);
      if ((eta_bin == -1) || (pt_bin == -1))
        continue;

      //if (TMath::Abs(DCA_x_mpd[track]) >= f_dca[0][pt_bin][eta_bin]->GetParameter(2)*2) continue;
      //if (TMath::Abs(DCA_y_mpd[track]) >= f_dca[1][pt_bin][eta_bin]->GetParameter(2)*2) continue;
      //if (TMath::Abs(DCA_z_mpd[track]) >= f_dca[2][pt_bin][eta_bin]->GetParameter(2)*2) continue;

      TF1 sigma_fit_X = *f_pt_fit[0][eta_bin];
      TF1 sigma_fit_Y = *f_pt_fit[1][eta_bin];
      TF1 sigma_fit_Z = *f_pt_fit[2][eta_bin];

      if (TMath::Abs(DCA_x_mpd[track]) >= sigma_fit_X(Pt) * 2)
        continue;
      if (TMath::Abs(DCA_y_mpd[track]) >= sigma_fit_Y(Pt) * 2)
        continue;
      if (TMath::Abs(DCA_z_mpd[track]) >= sigma_fit_Z(Pt) * 2)
        continue;

      /*if (pid_tpc_prob_pion_mpd[track] > 0.9 && 
        pid_tpc_prob_kaon_mpd[track] < 0.1 &&
        pid_tpc_prob_proton_mpd[track] < 0.1) isPion = kTRUE;

        if (pid_tpc_prob_pion_mpd[track] < 0.1 && 
        pid_tpc_prob_kaon_mpd[track] > 0.9 &&
        pid_tpc_prob_proton_mpd[track] < 0.1) isKaon = kTRUE;

        if (pid_tpc_prob_pion_mpd[track] < 0.1 && 
        pid_tpc_prob_kaon_mpd[track] < 0.1 &&
        pid_tpc_prob_proton_mpd[track] > 0.9) isProton = kTRUE;

        if (signed_pt_mpd[track] >= 0 ) continue; // Only positively charged particles being selected
        if (TOF_flag_mpd[track] == 0) continue; // Only tracks with TOF hits being selected
        if (TOF_flag_mpd[track] == 4) continue; // Only tracks with TOF hits being selected

        if (sort == 0 && !isProton)	continue;
        if (sort == 1 && !isKaon) 	continue;
        if (sort == 2 && !isPion)		continue;*/

      Float_t Rapidity = TMath::Log((TMath::Sqrt(p_mass2[sort] + Pt * Pt * TMath::CosH(Eta) * TMath::CosH(Eta)) + Pt * TMath::SinH(Eta)) / (TMath::Sqrt(p_mass2[sort] + Pt * Pt)));

      Double_t Psi_1_FULL_ZDC = GetPsiFullZdc(ZDC_energy_mpd, 1); //getting event plane using ZDC for 1st harmonic

      Float_t Pt_MC = Abs(pt_mc[id_from_mc_mpd[track]]);
      Float_t Eta_MC = eta_mc[id_from_mc_mpd[track]];
      Float_t Phi_MC = ATan2(py_mc[id_from_mc_mpd[track]], px_mc[id_from_mc_mpd[track]]);
      Float_t Rapidity_MC = TMath::Log((TMath::Sqrt(p_mass2[sort] + Pt_MC * Pt_MC * TMath::CosH(Eta_MC) * TMath::CosH(Eta_MC)) + Pt_MC * TMath::SinH(Eta_MC)) / (TMath::Sqrt(p_mass2[sort] + Pt_MC * Pt_MC)));

      int sign;
      if (Eta < 0)
        sign = -1;
      else
        sign = 1;
      TF1 res11 = *resolution_fit[0][0][1];
      Float_t mid_bin_cent = (centralityBinsRes[cenrality_bin_res] + centralityBinsRes[cenrality_bin_res + 1]) / 2;
      if (res11(mid_bin_cent) != 0)
      {
        if (cenrality_bin_flow != -1)
        {
          if ((Abs(Rapidity) > 0.2) && (Abs(Rapidity) < Cut_Eta_Max))
          {
            p_flow_wrt_full_vs_pt_divided[cenrality_bin_flow][0][0][1]->Fill(Pt, sign * Cos(Psi_1_FULL_ZDC - Phi) / res11(mid_bin_cent),wB);
            p_flow_wrt_RP_vs_pt[cenrality_bin_flow][0]->Fill(Pt, sign * Cos(phi_EP - Phi_MC),wB);
          }
          //if (Pt > Cut_Pt_Min)
          p_flow_wrt_full_vs_eta_divided[cenrality_bin_flow][0][0][1]->Fill(Eta, Cos(Psi_1_FULL_ZDC - Phi) / res11(mid_bin_cent),wB);
          p_flow_wrt_RP_vs_eta[cenrality_bin_flow][0]->Fill(Eta, Cos(phi_EP - Phi_MC),wB);
          //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min))
          p_flow_wrt_full_vs_rapidity_divided[cenrality_bin_flow][0][0][1]->Fill(Rapidity, Cos(Psi_1_FULL_ZDC - Phi) / res11(mid_bin_cent),wB);
          p_flow_wrt_RP_vs_rapidity[cenrality_bin_flow][0]->Fill(Rapidity, Cos(phi_EP - Phi_MC),wB);
        }
        //cout << "filling with centralityBinsRes[cenrality_bin_res] = " << centralityBinsRes[cenrality_bin_res] << " Cos(Psi_1_FULL_ZDC - Phi)/res11(mid_bin_cent) = " << Cos(Psi_1_FULL_ZDC - Phi)/res11(mid_bin_cent) << endl;
        //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min))
        p_flow_wrt_full_vs_centrality_divided[0][0][1]->Fill(centralityBinsRes[cenrality_bin_res] + 0.1, sign * Cos(Psi_1_FULL_ZDC - Phi) / res11(mid_bin_cent),wB);
        p_flow_wrt_RP_vs_centrality[0]->Fill(centralityBinsRes[cenrality_bin_res] + 0.1, sign * Cos(phi_EP - Phi_MC),wB);
      }
      TF1 res21 = *resolution_fit[1][0][1];
      if (res21(mid_bin_cent) != 0)
      {
        if (cenrality_bin_flow != -1)
        {
          //if (Abs(Eta) < Cut_Eta_Max)
          p_flow_wrt_full_vs_pt_divided[cenrality_bin_flow][1][0][1]->Fill(Pt, Cos(2. * (Psi_1_FULL_ZDC - Phi)) / res21(mid_bin_cent),wB);
          p_flow_wrt_RP_vs_pt[cenrality_bin_flow][1]->Fill(Pt, Cos(2. * (phi_EP - Phi_MC)),wB);
          //if (Pt > Cut_Pt_Min)
          p_flow_wrt_full_vs_eta_divided[cenrality_bin_flow][1][0][1]->Fill(Eta, Cos(2. * (Psi_1_FULL_ZDC - Phi)) / res21(mid_bin_cent),wB);
          p_flow_wrt_RP_vs_eta[cenrality_bin_flow][1]->Fill(Eta, Cos(2. * (phi_EP - Phi_MC)),wB);
          //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min))
          p_flow_wrt_full_vs_rapidity_divided[cenrality_bin_flow][1][0][1]->Fill(Rapidity, Cos(2. * (Psi_1_FULL_ZDC - Phi)) / res21(mid_bin_cent),wB);
          p_flow_wrt_RP_vs_rapidity[cenrality_bin_flow][1]->Fill(Rapidity, Cos(2. * (phi_EP - Phi_MC)),wB);
        }
        //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min))
        p_flow_wrt_full_vs_centrality_divided[1][0][1]->Fill(centralityBinsRes[cenrality_bin_res] + 0.1, Cos(2. * (Psi_1_FULL_ZDC - Phi)) / res21(mid_bin_cent),wB);
        p_flow_wrt_RP_vs_centrality[1]->Fill(centralityBinsRes[cenrality_bin_res] + 0.1, Cos(2. * (phi_EP - Phi_MC)),wB);
      }
    }
    ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    //GENERATED
    /*		for (Long64_t track = 0; track < n_tracks_mc; ++track)
                {
                if (PDG_code_mc[track] != pdg_codes[sort]) continue;

                Float_t Pt = Abs(pt_mc[track]);
                Float_t Eta = eta_mc[track];
                Float_t Phi = ATan2(py_mc[track],px_mc[track]);
                Float_t Rapidity = .5*TMath::Log((energy_mc[track] + pz_mc[track])/(energy_mc[track] - pz_mc[track]));

                Int_t pt_bin = GetPtBin(Pt);
                Int_t eta_bin = GetEtaBin(Eta);
                if ( (eta_bin == -1) || (pt_bin == -1 )) continue;
                if (mother_ID_mc[track] > -1 ) continue;

                int sign;
                if (Eta < 0) sign = -1;
                else sign = 1;

    //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min)) 
    //{
    p_flow_wrt_RP_vs_centrality[0]->Fill(centralityBinsRes[cenrality_bin_res] + 0.1, sign*Cos(phi_EP - Phi));
    //cout << "Filling with centralityBinsRes[cenrality_bin_res] + 0.1 = " << centralityBinsRes[cenrality_bin_res] + 0.1 << "Cos(phi_EP - Phi) = " << Cos(phi_EP - Phi) << endl;
    p_flow_wrt_RP_vs_centrality[1]->Fill(centralityBinsRes[cenrality_bin_res] + 0.1, Cos(2.*(phi_EP - Phi)));
    //}

    if (cenrality_bin_flow != -1)
    {
    if ((Abs(Eta) > 0.2) && (Abs(Eta) < Cut_Eta_Max))
    p_flow_wrt_RP_vs_pt[cenrality_bin_flow][0]->Fill(Pt , sign*Cos(phi_EP - Phi));
    //if (Abs(Eta) < Cut_Eta_Max)	
    p_flow_wrt_RP_vs_pt[cenrality_bin_flow][1]->Fill(Pt , Cos(2.*(phi_EP-Phi)));

    //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min))
    p_flow_wrt_RP_vs_rapidity[cenrality_bin_flow][0]->Fill(Rapidity , Cos(phi_EP	- Phi));
    //if ((Abs(Eta) < Cut_Eta_Max) && (Pt > Cut_Pt_Min))	
    p_flow_wrt_RP_vs_rapidity[cenrality_bin_flow][1]->Fill(Rapidity , Cos(2.*(phi_EP - Phi)));
    //if (Pt > Cut_Pt_Min)	
    p_flow_wrt_RP_vs_eta[cenrality_bin_flow][0]->Fill(Eta , Cos(phi_EP - Phi));
    //if (Pt > Cut_Pt_Min)	
    p_flow_wrt_RP_vs_eta[cenrality_bin_flow][1]->Fill(Eta , Cos(2.*(phi_EP - Phi)));
    }
    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    }*/
  }
}

void MpdCalculator::GetQsTpc(EPParticle *event_plane_buffer, Int_t buffer_size, Int_t sign, Int_t harm, Double_t &Qx, Double_t &Qy)
{
  Double_t Qcos = 0., Qsin = 0., total_momenta = 0;

  total_momenta = GetTotalMomenta(event_plane_buffer, buffer_size, sign);

  Int_t w_sign = 0;
  if (harm == 2)
    w_sign = 1;
  else
    w_sign = sign;

  for (int ep_particle = 0; ep_particle < buffer_size; ++ep_particle)
  {
    if (!(event_plane_buffer[ep_particle].Eta * sign > 0))
      continue;
    Qcos += w_sign * event_plane_buffer[ep_particle].Pt / total_momenta * Cos(harm * event_plane_buffer[ep_particle].Phi);
    Qsin += w_sign * event_plane_buffer[ep_particle].Pt / total_momenta * Sin(harm * event_plane_buffer[ep_particle].Phi);
  }

  Qx = Qcos;
  Qy = Qsin;
}

void MpdCalculator::GetQsZdc(Float_t *zdc_energy, Int_t zdc_ID, Int_t harm, Double_t &Qx, Double_t &Qy)
{
  Double_t *phi_angle_of_modules = GetAngles();
  Double_t Qcos = 0, Qsin = 0;
  Double_t w_sign = 0;
  Double_t total_energy = GetTotalEnergy(zdc_energy, zdc_ID);

  if (harm == 2)
    w_sign = 1;
  else if (harm == 1)
  {
    if (zdc_ID == 0)
      w_sign = 1;
    else if (zdc_ID == 1)
      w_sign = -1;
  }

  for (int module = _N_MODULES_TOTAL / 2 * zdc_ID; module < _N_MODULES_TOTAL / 2 * (zdc_ID + 1); ++module)
  {
    if ((module == 22) || (module == 67))
      continue;
    //		if ((i ==15)||(i ==21)||(i==23)||(i ==29)||(i==60)||
    //				(i==66)||(i==68)||(i==74)) continue;
    Qcos += w_sign * zdc_energy[module] / total_energy * Cos(harm * phi_angle_of_modules[module]);
    Qsin += w_sign * zdc_energy[module] / total_energy * Sin(harm * phi_angle_of_modules[module]);
  }

  Qx = Qcos;
  Qy = Qsin;
}

Double_t MpdCalculator::GetPsiHalfTpc(EPParticle *event_plane_buffer, Int_t buffer_size, Int_t sign, Int_t harm, Double_t &Qx, Double_t &Qy) //if sign == 1 then its positive pseudorapidity and all weights are positive,
{
  GetQsTpc(event_plane_buffer, buffer_size, sign, harm, Qx, Qy);
  Double_t PsiEP = (1 / (Double_t)harm) * ATan2(Qx, Qy);
  return PsiEP;
}

Double_t MpdCalculator::GetPsiHalfZdc(Float_t *zdc_energy, Int_t zdc_ID, Int_t n, Double_t &qx, Double_t &qy)
{
  Double_t Qcos = 0., Qsin = 0.;
  GetQsZdc(zdc_energy, zdc_ID, n, Qcos, Qsin);

  qx = Qcos;
  qy = Qsin;
  Double_t PsiEP = (1 / (Double_t)n) * ATan2(Qsin, Qcos);
  return PsiEP;
}

Double_t MpdCalculator::GetPsiFullZdc(Float_t *zdc_energy, Int_t n)
{
  Double_t QcosR = 0., QsinR = 0.;
  GetQsZdc(zdc_energy, 0, n, QcosR, QsinR);

  Double_t QcosL = 0., QsinL = 0.;
  GetQsZdc(zdc_energy, 1, n, QcosL, QsinL);

  Double_t psiEP = ATan2(QsinR + QsinL, QcosR + QcosL) / (Double_t)n; // (-pi/n,pi/n]
  return psiEP;
}

Double_t MpdCalculator::GetTotalMomenta(EPParticle *event_plane_buffer, Int_t buffer_size, Int_t sign)
{
  Double_t total_momenta = 0.;
  for (int ep_particle = 0; ep_particle < buffer_size; ++ep_particle)
  {
    if (!(event_plane_buffer[ep_particle].Eta * sign > 0))
      continue;
    total_momenta += event_plane_buffer[ep_particle].Pt;
  }
  return total_momenta;
}

Double_t MpdCalculator::GetTotalEnergy(Float_t *zdc_energy, Int_t zdc_ID)
{
  Double_t total_energy = 0.;
  for (int i = _N_MODULES_TOTAL / 2 * zdc_ID; i < _N_MODULES_TOTAL / 2 * (zdc_ID + 1); ++i)
  {
    if ((i == 22) || (i == 67))
      continue;
    //		if ((i ==15)||(i ==21)||(i==23)||(i ==29)||(i==60)||
    //				(i==66)||(i==68)||(i==74)) continue;
    total_energy += zdc_energy[i];
  }
  return total_energy;
}

Double_t MpdCalculator::GetPsiFullTpc(EPParticle *event_plane_buffer, Int_t buffer_size, Int_t harm)
{
  Double_t QcosR = 0., QsinR = 0.;
  GetQsTpc(event_plane_buffer, buffer_size, 1, harm, QcosR, QsinR);

  Double_t QcosL = 0., QsinL = 0.;
  GetQsTpc(event_plane_buffer, buffer_size, -1, harm, QcosL, QsinL);

  Double_t psiEP = ATan2(QsinR + QsinL, QcosR + QcosL) / (Double_t)harm; // (-pi/n,pi/n]
  return psiEP;
}

//~ _  _
//~ (_)| |
//~ __      __ _ __  _ | |_  ___
//~ \ \ /\ / /| '__|| || __|/ _ \
//~ \ V  V / | |   | || |_|  __/
//~ \_/\_/  |_|   |_| \__|\___|

void MpdCalculator::Write()
{
  outFile->cd();

  h_nhits_TPC->Write();
  h_nhits_TPC_after->Write();
  h_impact_parameter->Write();
  h_multiplicity->Write();
  h_multiplicity_before->Write();
  h2_energyZDC_L_vs_energy_ZDC_R->Write();
  p_b_vs_multiplicity->Write();
  p_b_vs_energy->Write();
  h2_b_vs_energy->Write();
  h2_b_vs_centrality->Write();
  h2_b_vs_multiplicity->Write();

  for (Int_t dim = 0; dim < Ndim; dim++)
  {
    for (Int_t pt_bin = 0; pt_bin < NptBins; pt_bin++)
    {
      for (Int_t eta_bin = 0; eta_bin < NetaBins; eta_bin++)
      {
        h_DCA_all[dim][pt_bin][eta_bin]->Write();
        h_DCA_primary[dim][pt_bin][eta_bin]->Write();
        h_DCA_secondary[dim][pt_bin][eta_bin]->Write();
      }
    }
  }

  for (int sort = 0; sort < _N_SORTS; ++sort)
  {
    p_momenta_resolution[sort]->Write();
  }

  for (int harm = 0; harm < _N_HARM; ++harm)
  {
    p_flow_wrt_RP_vs_centrality[harm]->Write();
    for (int centralityBin = 0; centralityBin < NcentralityBinsFlow; ++centralityBin)
    {
      p_flow_wrt_RP_vs_pt[centralityBin][harm]->Write();
      p_flow_wrt_RP_vs_eta[centralityBin][harm]->Write();
      p_flow_wrt_RP_vs_rapidity[centralityBin][harm]->Write();
    }

    for (int method = 0; method < _N_METHOD; ++method)
    {
      for (int _harm = 0; _harm < _N_HARM; ++_harm)
      {
        p_flow_wrt_full_vs_centrality[harm][_harm][method]->Write();
        p_flow_wrt_full_vs_centrality_divided[harm][_harm][method]->Write();
        p_Res2Psi_vs_b[harm][_harm][method]->Write();
        p_true_Res_vs_b[harm][_harm][method]->Write();
        for (int centralityBin = 0; centralityBin < NcentralityBinsFlow; ++centralityBin)
        {
          p_flow_wrt_full_vs_pt[centralityBin][harm][_harm][method]->Write();
          p_flow_wrt_full_vs_eta[centralityBin][harm][_harm][method]->Write();
          p_flow_wrt_full_vs_rapidity[centralityBin][harm][_harm][method]->Write();
          p_flow_wrt_full_vs_pt_divided[centralityBin][harm][_harm][method]->Write();
          p_flow_wrt_full_vs_eta_divided[centralityBin][harm][_harm][method]->Write();
          p_flow_wrt_full_vs_rapidity_divided[centralityBin][harm][_harm][method]->Write();
        }
      }
    }
  }

  for (int arm = 0; arm < _N_ARM; ++arm)
  {
    h_energy_ZDC_total[arm]->Write();
    h2_energy_ZDC_vs_multiplicity[arm]->Write();

    for (int harm = 0; harm < _N_HARM; ++harm)
    {
      for (int method = 0; method < _N_METHOD; ++method)
      {
        p_qx_vs_b[arm][harm][method]->Write();
        p_qy_vs_b[arm][harm][method]->Write();
        for (int _harm = 0; _harm < _N_HARM; ++_harm)
        {
          p_true_Res_half_vs_b[arm][harm][_harm][method]->Write();
        }
      }
    }
  }

  for (int centralityBin = 0; centralityBin < NcentralityBinsRes; ++centralityBin)
  {
    for (int sort = 0; sort < _N_SORTS; ++sort)
    {
      h_phi[centralityBin][sort]->Write();
      h_pt[centralityBin][sort]->Write();
      h_eta[centralityBin][sort]->Write();
      h2_pt_vs_eta[centralityBin][sort]->Write();
      h2_pt_vs_phi[centralityBin][sort]->Write();
      h2_phi_vs_eta[centralityBin][sort]->Write();
      h_phi_after[centralityBin][sort]->Write();
      h_pt_after[centralityBin][sort]->Write();
      h_eta_after[centralityBin][sort]->Write();
      h2_pt_vs_eta_after[centralityBin][sort]->Write();
      h2_pt_vs_phi_after[centralityBin][sort]->Write();
      h2_phi_vs_eta_after[centralityBin][sort]->Write();

      h_phi_mc[centralityBin][sort]->Write();
      h_pt_mc[centralityBin][sort]->Write();
      h_eta_mc[centralityBin][sort]->Write();
      h_phi_mc_after[centralityBin][sort]->Write();
      h_pt_mc_after[centralityBin][sort]->Write();
      h_eta_mc_after[centralityBin][sort]->Write();
    }
  }
  outFile->Close();
}

Int_t MpdCalculator::GetCentralityBinRes(Int_t centrality)
{
  int centrality_bin = -1;
  for (int c_bin = 0; c_bin < NcentralityBinsRes; ++c_bin)
  {
    if ((centrality > centralityBinsRes[c_bin]) && (centrality <= centralityBinsRes[c_bin + 1]))
      centrality_bin = c_bin;
  }
  return centrality_bin;
}

Int_t MpdCalculator::GetCentralityBinFlow(Int_t centrality)
{
  int centrality_bin = -1;
  for (int c_bin = 0; c_bin < NcentralityBinsFlow; ++c_bin)
    if ((centrality > centralityBinsFlow[c_bin]) && (centrality <= centralityBinsFlow[c_bin + 1]))
      centrality_bin = c_bin;
  if (centrality_bin == 1)
    centrality_bin = -1;
  else if (centrality_bin == 2)
    centrality_bin = 1; ////////////////////PODGON!!!!
  return centrality_bin;
}

Int_t MpdCalculator::GetMultiplicityTPC() //should be called in loop over events
{
  Int_t multiplicity = 0;
  for (Long64_t track = 0; track < n_tracks_mpd; ++track) //loop over mpdtracks, cut on mother ID will be everywhere
  {
    if (id_from_mc_mpd[track] == -1)
      continue;			//equivalent to mother id cut
    multiplicity++; //multiplicity in TPC before eta, nhits and pt cuts, but after mother ID cut
  }
  return multiplicity;
}

int MpdCalculator::GetPtBin(Float_t pt)
{
  int pt_bin = -1;
  for (int i = 0; i < NptBins; ++i)
    if ((pt > ptBins[i]) && (pt <= ptBins[i + 1]))
      pt_bin = i;
  return pt_bin;
}

int MpdCalculator::GetEtaBin(Float_t eta)
{
  int eta_bin = -1;
  for (int i = 0; i < NetaBins; ++i)
    if ((eta > etaBins[i]) && (eta <= etaBins[i + 1]))
      eta_bin = i;
  return eta_bin;
}