Browse Source

Resolution & basic EP QA calculations were implemented.
All EP & centrality info is stored in mfEventInfo, main class to
provide mfEventInfo is mfProcessEvent.

PeterParfenov 2 years ago
parent
commit
259d103c90

+ 3 - 0
CMakeLists.txt

@@ -38,6 +38,9 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O2 -Wall")
 set(MPD_FLOW_INCLUDE_DIRECTORIES
   ${CMAKE_CURRENT_SOURCE_DIR}
   "${CMAKE_CURRENT_SOURCE_DIR}/utils"
+  "${CMAKE_CURRENT_SOURCE_DIR}/qvector"
+  "${CMAKE_CURRENT_SOURCE_DIR}/centralityreader"
+  "${CMAKE_CURRENT_SOURCE_DIR}/eventanalysis"
   "${CMAKE_CURRENT_SOURCE_DIR}/bin"
   ${ROOT_INCLUDE_DIRS}
   ${PICO_DST_INC}

+ 17 - 3
MpdFlow.LinkDef.h

@@ -9,10 +9,24 @@
 #pragma link C++ nestedclass;
 #pragma link C++ nestedtypedef;
 
-#pragma link C++ namespace mpdflow;
-#pragma link C++ function mpdflow::GetPtBin(float);
-#pragma link C++ function mpdflow::GetEtaBin(float);
+#pragma link C++ namespace mfdca;
+#pragma link C++ function mfdca::GetPtBin(float);
+#pragma link C++ function mfdca::GetEtaBin(float);
+
+#pragma link C++ namespace mfcent;
+#pragma link C++ function mfcent::CentB(float);
+#pragma link C++ function mfcent::isRefMult(float,float);
+
+#pragma link C++ namespace mfEp;
+#pragma link C++ function mfEp::GetFHCalPhi(int);
+#pragma link C++ namespace mfRes;
 
 #pragma link C++ class mfPicoDstReader+;
+#pragma link C++ class mfCentralityInfo+;
+#pragma link C++ class mfQvector+;
+#pragma link C++ class mfEventInfo+;
+#pragma link C++ class mfDcaReader+;
+#pragma link C++ class mfInputReader+;
+#pragma link C++ class mfProcessEvent+;
 
 #endif

+ 7 - 1
bin/CMakeLists.txt

@@ -3,4 +3,10 @@ add_executable(get-dca get-dca.cpp)
 target_link_libraries(get-dca MpdFlow ${MPD_FLOW_INCLUDE_LIBRARIES})
 
 add_executable(fit-dca fit-dca.cpp)
-target_link_libraries(fit-dca MpdFlow ${ROOT_LIBRARIES})
+target_link_libraries(fit-dca MpdFlow ${ROOT_LIBRARIES})
+
+add_executable(get-res get-res.cpp)
+target_link_libraries(get-res MpdFlow ${MPD_FLOW_INCLUDE_LIBRARIES})
+
+add_executable(get-qv get-qv.cpp)
+target_link_libraries(get-qv MpdFlow ${MPD_FLOW_INCLUDE_LIBRARIES})

+ 12 - 10
bin/fit-dca.cpp

@@ -9,12 +9,12 @@
 #include <TFile.h>
 #include <iostream>
 
-using mpdflow::etaBins;
-using mpdflow::Ndim;
-using mpdflow::NetaBins;
-using mpdflow::NptBins;
-using mpdflow::ptBins;
-using mpdflow::NdcaFitIter;
+using mfdca::etaBins;
+using mfdca::Ndim;
+using mfdca::NetaBins;
+using mfdca::NptBins;
+using mfdca::ptBins;
+using mfdca::NdcaFitIter;
 
 int main(int argc, char **argv)
 {
@@ -73,7 +73,7 @@ int main(int argc, char **argv)
       {
         h_dca[dim][ptbin][etabin] = (TH1F *)fi->Get(Form("h_dca_%i_%i_%i", dim, ptbin, etabin));
         dca_fit[dim][ptbin][etabin] = new TF1(Form("f_dca_%i_%i_%i", dim, ptbin, etabin), "gaus", -0.2, 0.2);
-        h_dca[dim][ptbin][etabin]->Fit(Form("f_dca_%i_%i_%i", dim, ptbin, etabin), "R");
+        h_dca[dim][ptbin][etabin]->Fit(Form("f_dca_%i_%i_%i", dim, ptbin, etabin), "RQ");
       }
     }
   }
@@ -115,10 +115,11 @@ int main(int argc, char **argv)
   for (Int_t i_dim = 0; i_dim < Ndim; i_dim++)
   {
     // Fit N iterations
-    for (int it=0; it<NdcaFitIter; it++)
+    h_sigma[i_dim]->Fit(f_sigma[i_dim], "R");
+    par = f_sigma[i_dim]->GetParameters();
+    for (int it=0; it<NdcaFitIter-1; it++)
     {
-      std::cout << "\tIteration " << it << std::endl;
-      if (it != 0) f_sigma[i_dim]->SetParameters(par);
+      f_sigma[i_dim]->SetParameters(par);
       h_sigma[i_dim]->Fit(f_sigma[i_dim], "R");
       par = f_sigma[i_dim]->GetParameters();
     }
@@ -150,4 +151,5 @@ int main(int argc, char **argv)
 
   timer.Stop();
   timer.Print();
+  return 0;
 }

+ 16 - 15
bin/get-dca.cpp

@@ -13,14 +13,14 @@
 #include <PicoDstRecoTrack.h>
 #include <PicoDstFHCal.h>
 
-using mpdflow::etaBins;
-using mpdflow::Ndim;
-using mpdflow::NetaBins;
-using mpdflow::NptBins;
-using mpdflow::ptBins;
+using mfdca::etaBins;
+using mfdca::Ndim;
+using mfdca::NetaBins;
+using mfdca::NptBins;
+using mfdca::ptBins;
 
-using mpdflow::GetEtaBin;
-using mpdflow::GetPtBin;
+using mfdca::GetEtaBin;
+using mfdca::GetPtBin;
 
 int main(int argc, char **argv)
 {
@@ -94,17 +94,17 @@ int main(int argc, char **argv)
     }
   }
 
-  PicoDstMCEvent *mcEvent = nullptr;
-  PicoDstRecoEvent *recoEvent = nullptr;
+  // PicoDstMCEvent *mcEvent = nullptr;
+  // PicoDstRecoEvent *recoEvent = nullptr;
   TClonesArray *recoTracks = nullptr;
-  TClonesArray *mcTracks = nullptr;
-  TClonesArray *fhcalmodules = nullptr;
+  // TClonesArray *mcTracks = nullptr;
+  // TClonesArray *fhcalmodules = nullptr;
 
-  inChain->SetBranchAddress("mcevent.", &mcEvent);
-  inChain->SetBranchAddress("recoevent.", &recoEvent);
-  inChain->SetBranchAddress("mctracks", &mcTracks);
+  // inChain->SetBranchAddress("mcevent.", &mcEvent);
+  // inChain->SetBranchAddress("recoevent.", &recoEvent);
+  // inChain->SetBranchAddress("mctracks", &mcTracks);
   inChain->SetBranchAddress("recotracks", &recoTracks);
-  inChain->SetBranchAddress("FHCalModules", &fhcalmodules);
+  // inChain->SetBranchAddress("FHCalModules", &fhcalmodules);
 
   // Start event loop
   int n_entries = inChain->GetEntries();
@@ -150,4 +150,5 @@ int main(int argc, char **argv)
 
   timer.Stop();
   timer.Print();
+  return 0;
 }

+ 239 - 0
bin/get-qv.cpp

@@ -0,0 +1,239 @@
+#include "mfPicoDstReader.h"
+#include "mfNamespace.h"
+#include "mfProcessEvent.h"
+#include "mfInputReader.h"
+
+#include <TStopwatch.h>
+#include <TH2F.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <iostream>
+#include <TMath.h>
+
+#include <PicoDstMCEvent.h>
+#include <PicoDstRecoEvent.h>
+#include <PicoDstMCTrack.h>
+#include <PicoDstRecoTrack.h>
+#include <PicoDstFHCal.h>
+
+using mfdca::Ndim;
+using mfEp::EpNames;
+using mfEp::NetaGaps;
+using mfEp::NsubEvents;
+using mfEp::RangeAngle;
+using mfEp::TpcEtaGap;
+using mfEp::TpcLEp;
+using mfEp::TpcREp;
+using mfEp::EpNames::kFHCal;
+using mfEp::EpNames::kFHCal_L;
+using mfEp::EpNames::kFHCal_R;
+using mfRes::NresBin;
+
+using TMath::Cos;
+using TMath::Pi;
+using TMath::Sin;
+
+int main(int argc, char **argv)
+{
+  TString iFileName, oFileName, dcaFileName;
+
+  if (argc < 7)
+  {
+    std::cerr << "./get-qv -i INPUTFILE/LIST -o OUTPUT.root -dca DCA_FIT.root" << std::endl;
+    return 1;
+  }
+  for (int i = 1; i < argc; i++)
+  {
+    if (std::string(argv[i]) != "-i" &&
+        std::string(argv[i]) != "-o" &&
+        std::string(argv[i]) != "-dca")
+    {
+      std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
+      return 2;
+    }
+    else
+    {
+      if (std::string(argv[i]) == "-i" && i != argc - 1)
+      {
+        iFileName = argv[++i];
+        continue;
+      }
+      if (std::string(argv[i]) == "-i" && i == argc - 1)
+      {
+        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
+        return 1;
+      }
+      if (std::string(argv[i]) == "-o" && i != argc - 1)
+      {
+        oFileName = argv[++i];
+        continue;
+      }
+      if (std::string(argv[i]) == "-o" && i == argc - 1)
+      {
+        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
+        return 1;
+      }
+      if (std::string(argv[i]) == "-dca" && i != argc - 1)
+      {
+        dcaFileName = argv[++i];
+        continue;
+      }
+      if (std::string(argv[i]) == "-dca" && i == argc - 1)
+      {
+        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
+        return 1;
+      }
+    }
+  }
+
+  TStopwatch timer;
+  timer.Start();
+
+  mfPicoDstReader *input = new mfPicoDstReader();
+  TChain *inChain = (TChain *)input->Open(iFileName);
+
+  if (!inChain)
+  {
+    std::cerr << "[ERROR]: nothing to read. Maybe the input TTree(s) is(are) empty?" << std::endl;
+    return 10;
+  }
+
+  // Init output profiles
+  TH2F *hQxVsCent[NsubEvents];
+  TH2F *hQyVsCent[NsubEvents];
+  TH2F *hPsiEpVsCent[NsubEvents];
+  TH2F *hQxRLVsCent[NetaGaps + 1];
+  TH2F *hQyRLVsCent[NetaGaps + 1];
+  TH2F *hPsiEpRLVsCent[NetaGaps + 1];
+  TH2F *hPsiRpRVsCent[NetaGaps + 1];
+  TH2F *hPsiRpLVsCent[NetaGaps + 1];
+  for (int i = 0; i < NsubEvents; i++)
+  {
+    hQxVsCent[i] = new TH2F(Form("hQxVsCent_%i", i),
+                            Form("hQxVsCent_%i", i),
+                            600, -1.1, 1.1, NresBin, 0., 100.);
+    hQyVsCent[i] = new TH2F(Form("hQyVsCent_%i", i),
+                            Form("hQyVsCent_%i", i),
+                            600, -1.1, 1.1, NresBin, 0., 100.);
+    hPsiEpVsCent[i] = new TH2F(Form("hPsiEpVsCent_%i", i),
+                               Form("hPsiEpVsCent_%i", i),
+                               360, 0., 2 * Pi(), NresBin, 0., 100.);
+  }
+  for (int i = 0; i < NetaGaps + 1; i++)
+  {
+    hQxRLVsCent[i] = new TH2F(Form("hQxRLVsCent_%i", i),
+                              Form("hQxRLVsCent_%i", i),
+                              600, -1.1, 1.1, NresBin, 0., 100.);
+    hQyRLVsCent[i] = new TH2F(Form("hQyRLVsCent_%i", i),
+                              Form("hQyRLVsCent_%i", i),
+                              600, -1.1, 1.1, NresBin, 0., 100.);
+    hPsiEpRLVsCent[i] = new TH2F(Form("hPsiEpRLVsCent_%i", i),
+                                 Form("hPsiEpRLVsCent_%i", i),
+                                 360, -Pi(), Pi(), NresBin, 0., 100.);
+    hPsiRpRVsCent[i] = new TH2F(Form("hPsiRpRVsCent_%i", i),
+                                Form("hPsiRpRVsCent_%i", i),
+                                360, -Pi(), Pi(), NresBin, 0., 100.);
+    hPsiRpLVsCent[i] = new TH2F(Form("hPsiRpLVsCent_%i", i),
+                                Form("hPsiRpLVsCent_%i", i),
+                                360, -Pi(), Pi(), NresBin, 0., 100.);
+  }
+
+  PicoDstMCEvent *mcEvent = nullptr;
+  // PicoDstRecoEvent *recoEvent = nullptr;
+  TClonesArray *recoTracks = nullptr;
+  // TClonesArray *mcTracks = nullptr;
+  TClonesArray *fhcalmodules = nullptr;
+
+  inChain->SetBranchAddress("mcevent.", &mcEvent);
+  // inChain->SetBranchAddress("recoevent.", &recoEvent);
+  // inChain->SetBranchAddress("mctracks", &mcTracks);
+  inChain->SetBranchAddress("recotracks", &recoTracks);
+  inChain->SetBranchAddress("FHCalModules", &fhcalmodules);
+
+  // Start event loop
+  int n_entries = inChain->GetEntries();
+  mfInputReader *inputReader = new mfInputReader();
+  inputReader->SetDcaFitFile(dcaFileName);
+  for (int iEv = 0; iEv < n_entries; iEv++)
+  {
+    if (iEv % 1000 == 0)
+      std::cout << "Event [" << iEv << "/" << n_entries << "]" << std::endl;
+    inChain->GetEntry(iEv);
+
+    mfProcessEvent *procEvent = new mfProcessEvent();
+    float cent = procEvent->GetCentB(mcEvent);
+    procEvent->SetInputReader(inputReader);
+
+    // std::cout << "\tNtracks = " << recoTracks->GetEntriesFast() << std::endl;
+
+    mfEventInfo evInfo = procEvent->ProcessRecoEP(recoTracks, fhcalmodules);
+    if (evInfo.GetQvectors().size() == 0)
+      continue;
+
+    for (int iGap = 0; iGap < NetaGaps; iGap++)
+    {
+      // Ntracks >= 2 for TpcR EP
+      if (evInfo.GetEpMult(TpcREp[iGap]) < 2)
+        continue;
+      // Ntracks >= 2 for TpcL EP
+      if (evInfo.GetEpMult(TpcLEp[iGap]) < 2)
+        continue;
+
+      hQxVsCent[TpcREp[iGap]]->Fill(evInfo.GetQvector(TpcREp[iGap]).X(), cent);
+      hQyVsCent[TpcREp[iGap]]->Fill(evInfo.GetQvector(TpcREp[iGap]).Y(), cent);
+      hQxVsCent[TpcLEp[iGap]]->Fill(evInfo.GetQvector(TpcLEp[iGap]).X(), cent);
+      hQyVsCent[TpcLEp[iGap]]->Fill(evInfo.GetQvector(TpcLEp[iGap]).Y(), cent);
+
+      hPsiEpVsCent[TpcREp[iGap]]->Fill(evInfo.GetQvector(TpcREp[iGap]).GetPhi(), cent);
+      hPsiEpVsCent[TpcLEp[iGap]]->Fill(evInfo.GetQvector(TpcLEp[iGap]).GetPhi(), cent);
+
+      hQxRLVsCent[iGap]->Fill(evInfo.GetQvector(TpcREp[iGap]).X() - evInfo.GetQvector(TpcLEp[iGap]).X(), cent);
+      hQyRLVsCent[iGap]->Fill(evInfo.GetQvector(TpcREp[iGap]).Y() - evInfo.GetQvector(TpcLEp[iGap]).Y(), cent);
+
+      hPsiEpRLVsCent[iGap]->Fill(RangeAngle(evInfo.GetQvector(TpcREp[iGap]).GetPhi() - evInfo.GetQvector(TpcLEp[iGap]).GetPhi()), cent);
+      hPsiRpRVsCent[iGap]->Fill(RangeAngle(evInfo.GetQvector(TpcREp[iGap]).GetPhi() - mcEvent->GetPhiRP()), cent);
+      hPsiRpLVsCent[iGap]->Fill(RangeAngle(evInfo.GetQvector(TpcLEp[iGap]).GetPhi() - mcEvent->GetPhiRP()), cent);
+    }
+
+    hQxVsCent[kFHCal_R]->Fill(evInfo.GetQvector(kFHCal_R).X(), cent);
+    hQyVsCent[kFHCal_R]->Fill(evInfo.GetQvector(kFHCal_R).Y(), cent);
+    hQxVsCent[kFHCal_L]->Fill(evInfo.GetQvector(kFHCal_L).X(), cent);
+    hQyVsCent[kFHCal_L]->Fill(evInfo.GetQvector(kFHCal_L).Y(), cent);
+    hQxVsCent[kFHCal]->Fill(evInfo.GetQvector(kFHCal).X(), cent);
+    hQyVsCent[kFHCal]->Fill(evInfo.GetQvector(kFHCal).Y(), cent);
+
+    hPsiEpVsCent[kFHCal_R]->Fill(evInfo.GetQvector(kFHCal_R).GetPhi(), cent);
+    hPsiEpVsCent[kFHCal_L]->Fill(evInfo.GetQvector(kFHCal_L).GetPhi(), cent);
+    hPsiEpVsCent[kFHCal]->Fill(evInfo.GetQvector(kFHCal).GetPhi(), cent);
+
+    hQxRLVsCent[NetaGaps]->Fill(evInfo.GetQvector(kFHCal_R).X() - evInfo.GetQvector(kFHCal_L).X(), cent);
+    hQyRLVsCent[NetaGaps]->Fill(evInfo.GetQvector(kFHCal_R).Y() - evInfo.GetQvector(kFHCal_L).Y(), cent);
+
+    hPsiEpRLVsCent[NetaGaps]->Fill(RangeAngle(evInfo.GetQvector(kFHCal_R).GetPhi() - evInfo.GetQvector(kFHCal_L).GetPhi()), cent);
+    hPsiRpRVsCent[NetaGaps]->Fill(RangeAngle(evInfo.GetQvector(kFHCal_R).GetPhi() - mcEvent->GetPhiRP()), cent);
+    hPsiRpLVsCent[NetaGaps]->Fill(RangeAngle(evInfo.GetQvector(kFHCal_L).GetPhi() - mcEvent->GetPhiRP()), cent);
+  }
+
+  // Init output file
+  TFile *fo = new TFile(oFileName.Data(), "recreate");
+  fo->cd();
+  for (int i = 0; i < NsubEvents; i++)
+  {
+    hQxVsCent[i]->Write();
+    hQyVsCent[i]->Write();
+    hPsiEpVsCent[i]->Write();
+  }
+  for (int i = 0; i < NetaGaps + 1; i++)
+  {
+    hQxRLVsCent[i]->Write();
+    hQyRLVsCent[i]->Write();
+    hPsiEpRLVsCent[i]->Write();
+    hPsiRpRVsCent[i]->Write();
+    hPsiRpLVsCent[i]->Write();
+  }
+  fo->Close();
+
+  timer.Stop();
+  timer.Print();
+  return 0;
+}

+ 228 - 0
bin/get-res.cpp

@@ -0,0 +1,228 @@
+#include "mfPicoDstReader.h"
+#include "mfNamespace.h"
+#include "mfProcessEvent.h"
+#include "mfInputReader.h"
+
+#include <TStopwatch.h>
+#include <TProfile.h>
+#include <TClonesArray.h>
+#include <TFile.h>
+#include <iostream>
+#include <TMath.h>
+
+#include <PicoDstMCEvent.h>
+#include <PicoDstRecoEvent.h>
+#include <PicoDstMCTrack.h>
+#include <PicoDstRecoTrack.h>
+#include <PicoDstFHCal.h>
+
+using mfdca::Ndim;
+using mfEp::EpNames;
+using mfEp::NetaGaps;
+using mfEp::TpcEtaGap;
+using mfEp::TpcLEp;
+using mfEp::TpcREp;
+using mfEp::EpNames::kFHCal;
+using mfEp::EpNames::kFHCal_L;
+using mfEp::EpNames::kFHCal_R;
+using mfRes::NresBin;
+
+using TMath::ATan2;
+using TMath::Cos;
+using TMath::Pi;
+using TMath::Sin;
+
+int main(int argc, char **argv)
+{
+  TString iFileName, oFileName, dcaFileName;
+
+  if (argc < 7)
+  {
+    std::cerr << "./get-res -i INPUTFILE/LIST -o OUTPUT.root -dca DCA_FIT.root" << std::endl;
+    return 1;
+  }
+  for (int i = 1; i < argc; i++)
+  {
+    if (std::string(argv[i]) != "-i" &&
+        std::string(argv[i]) != "-o" &&
+        std::string(argv[i]) != "-dca")
+    {
+      std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
+      return 2;
+    }
+    else
+    {
+      if (std::string(argv[i]) == "-i" && i != argc - 1)
+      {
+        iFileName = argv[++i];
+        continue;
+      }
+      if (std::string(argv[i]) == "-i" && i == argc - 1)
+      {
+        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
+        return 1;
+      }
+      if (std::string(argv[i]) == "-o" && i != argc - 1)
+      {
+        oFileName = argv[++i];
+        continue;
+      }
+      if (std::string(argv[i]) == "-o" && i == argc - 1)
+      {
+        std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
+        return 1;
+      }
+      if (std::string(argv[i]) == "-dca" && i != argc - 1)
+      {
+        dcaFileName = argv[++i];
+        continue;
+      }
+      if (std::string(argv[i]) == "-dca" && i == argc - 1)
+      {
+        std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
+        return 1;
+      }
+    }
+  }
+
+  TStopwatch timer;
+  timer.Start();
+
+  mfPicoDstReader *input = new mfPicoDstReader();
+  TChain *inChain = (TChain *)input->Open(iFileName);
+
+  if (!inChain)
+  {
+    std::cerr << "[ERROR]: nothing to read. Maybe the input TTree(s) is(are) empty?" << std::endl;
+    return 10;
+  }
+
+  // Init output profiles
+  TProfile *pResRecoSqrTpc[NetaGaps];
+  TProfile *pResRpTpc[NetaGaps];
+  TProfile *pResRpTpcR[NetaGaps];
+  TProfile *pResRpTpcL[NetaGaps];
+  TProfile *pResRecoSqrFHCal;
+  TProfile *pResRpFHCal;
+  TProfile *pResRpFHCalR;
+  TProfile *pResRpFHCalL;
+  for (int i = 0; i < NetaGaps; i++)
+  {
+    pResRecoSqrTpc[i] = new TProfile(Form("pResRecoSqrtTpc_gap%i", i),
+                                     Form("pResRecoSqrtTpc_gap%i", i),
+                                     NresBin, 0., 100.);
+    pResRpTpc[i] = new TProfile(Form("pResRpTpc_gap%i", i),
+                                Form("pResRpTpc_gap%i", i),
+                                NresBin, 0, 100);
+    pResRpTpcR[i] = new TProfile(Form("pResRpTpcR_gap%i", i),
+                                 Form("pResRpTpcR_gap%i", i),
+                                 NresBin, 0, 100);
+    pResRpTpcL[i] = new TProfile(Form("pResRpTpcL_gap%i", i),
+                                 Form("pResRpTpcL_gap%i", i),
+                                 NresBin, 0, 100);
+  }
+  pResRecoSqrFHCal = new TProfile(Form("pResRecoSqrFHCal"),
+                                  Form("pResRecoSqrFHCal"),
+                                  NresBin, 0., 100.);
+  pResRpFHCal = new TProfile(Form("pResRpFHCal"),
+                             Form("pResRpFHCal"),
+                             NresBin, 0., 100.);
+  pResRpFHCalR = new TProfile(Form("pResRpFHCalR"),
+                              Form("pResRpFHCalR"),
+                              NresBin, 0., 100.);
+  pResRpFHCalL = new TProfile(Form("pResRpFHCalL"),
+                              Form("pResRpFHCalL"),
+                              NresBin, 0., 100.);
+
+  PicoDstMCEvent *mcEvent = nullptr;
+  // PicoDstRecoEvent *recoEvent = nullptr;
+  TClonesArray *recoTracks = nullptr;
+  // TClonesArray *mcTracks = nullptr;
+  TClonesArray *fhcalmodules = nullptr;
+
+  inChain->SetBranchAddress("mcevent.", &mcEvent);
+  // inChain->SetBranchAddress("recoevent.", &recoEvent);
+  // inChain->SetBranchAddress("mctracks", &mcTracks);
+  inChain->SetBranchAddress("recotracks", &recoTracks);
+  inChain->SetBranchAddress("FHCalModules", &fhcalmodules);
+
+  // Start event loop
+  int n_entries = inChain->GetEntries();
+  mfInputReader *inputReader = new mfInputReader();
+  inputReader->SetDcaFitFile(dcaFileName);
+  for (int iEv = 0; iEv < n_entries; iEv++)
+  {
+    if (iEv % 1000 == 0)
+      std::cout << "Event [" << iEv << "/" << n_entries << "]" << std::endl;
+    inChain->GetEntry(iEv);
+
+    mfProcessEvent *procEvent = new mfProcessEvent();
+    float cent = procEvent->GetCentB(mcEvent);
+    procEvent->SetInputReader(inputReader);
+
+    // std::cout << "\tNtracks = " << recoTracks->GetEntriesFast() << std::endl;
+
+    mfEventInfo evInfo = procEvent->ProcessRecoEP(recoTracks, fhcalmodules);
+    if (evInfo.GetQvectors().size() == 0)
+      continue;
+
+    for (int iGap = 0; iGap < NetaGaps; iGap++)
+    {
+      // Ntracks >= 2 for TpcR EP
+      if (evInfo.GetEpMult(TpcREp[iGap]) < 2)
+        continue;
+      // Ntracks >= 2 for TpcL EP
+      if (evInfo.GetEpMult(TpcLEp[iGap]) < 2)
+        continue;
+
+      double PsiR = evInfo.GetQvector(TpcREp[iGap]).GetPhi();
+      double PsiL = evInfo.GetQvector(TpcLEp[iGap]).GetPhi();
+      double resReco = Cos(2. * (PsiR - PsiL));
+      double resRpR = Cos(2. * (PsiR - mcEvent->GetPhiRP()));
+      double resRpL = Cos(2. * (PsiL - mcEvent->GetPhiRP()));
+
+      pResRecoSqrTpc[iGap]->Fill(cent, resReco);
+      pResRpTpcR[iGap]->Fill(cent, resRpR);
+      pResRpTpcL[iGap]->Fill(cent, resRpL);
+      pResRpTpc[iGap]->Fill(cent, resRpR);
+      pResRpTpc[iGap]->Fill(cent, resRpL);
+    }
+
+    double PsiFHCalR = evInfo.GetQvector(kFHCal_R).GetPhi();
+    double PsiFHCalL = evInfo.GetQvector(kFHCal_L).GetPhi();
+    double PsiFHCal = evInfo.GetQvector(kFHCal).GetPhi();
+    // double qx = evInfo.GetQvector(kFHCal_R).X() + evInfo.GetQvector(kFHCal_L).X();
+    // double qy = evInfo.GetQvector(kFHCal_R).Y() + evInfo.GetQvector(kFHCal_L).Y();
+    // double PsiFHCal = ATan2(qy, qx) / evInfo.GetQvector(kFHCal).GetHarmonic();
+    double resEpFHCal = Cos(1. * (PsiFHCalR - PsiFHCalL));
+    double resRpFHCal = Cos(1. * (PsiFHCal - mcEvent->GetPhiRP()));
+    double resRpFHCalR= Cos(1. * (PsiFHCalR - mcEvent->GetPhiRP()));
+    double resRpFHCalL= Cos(1. * (PsiFHCalL - mcEvent->GetPhiRP()));
+
+    pResRecoSqrFHCal->Fill(cent, resEpFHCal);
+    pResRpFHCal->Fill(cent, resRpFHCal);
+    pResRpFHCalR->Fill(cent, resRpFHCalR);
+    pResRpFHCalL->Fill(cent, resRpFHCalL);
+  }
+
+  // Init output file
+  TFile *fo = new TFile(oFileName.Data(), "recreate");
+  fo->cd();
+  for (int iGap = 0; iGap < NetaGaps; iGap++)
+  {
+    pResRecoSqrTpc[iGap]->Write();
+    pResRpTpcR[iGap]->Write();
+    pResRpTpcL[iGap]->Write();
+    pResRpTpc[iGap]->Write();
+    pResRpTpc[iGap]->Write();
+  }
+  pResRecoSqrFHCal->Write();
+  pResRpFHCal->Write();
+  pResRpFHCalR->Write();
+  pResRpFHCalL->Write();
+  fo->Close();
+
+  timer.Stop();
+  timer.Print();
+  return 0;
+}

+ 19 - 0
centralityreader/mfCentralityInfo.cxx

@@ -0,0 +1,19 @@
+#include "mfCentralityInfo.h"
+
+ClassImp(mfCentralityInfo);
+
+mfCentralityInfo::mfCentralityInfo(/* args */)
+{
+}
+
+mfCentralityInfo::~mfCentralityInfo()
+{
+}
+
+void mfCentralityInfo::Set(Float_t centMin, Float_t centMax, Float_t multMin, Float_t multMax)
+{
+  fCentMin = centMin;
+  fCentMax = centMax;
+  fMultMin = multMin;
+  fMultMax = multMax;
+}

+ 36 - 0
centralityreader/mfCentralityInfo.h

@@ -0,0 +1,36 @@
+#ifndef MPD_FLOW_CENTRALITYINFO_H
+#define MPD_FLOW_CENTRALITYINFO_H
+
+#include <iostream>
+#include <Rtypes.h>
+
+class mfCentralityInfo
+{
+private:
+  Float_t fCentMin;
+  Float_t fCentMax;
+  Float_t fMultMin;
+  Float_t fMultMax;
+
+public:
+  mfCentralityInfo(/* args */);
+  virtual ~mfCentralityInfo();
+
+  // Setters
+  void Set(Float_t centMin, Float_t centMax, Float_t multMin, Float_t multMax);
+  void SetCentMin(Float_t _a) { fCentMin = _a; }
+  void SetCentMax(Float_t _a) { fCentMax = _a; }
+  void SetMultMin(Float_t _a) { fMultMin = _a; }
+  void SetMultMax(Float_t _a) { fMultMax = _a; }
+
+  // Getters
+  Float_t GetCentMin() const { return fCentMin; }
+  Float_t GetCentMax() const { return fCentMax; }
+  Float_t GetCentMean() const { return 0.5 * (fCentMax - fCentMin); }
+  Float_t GetMultMin() const { return fMultMin; }
+  Float_t GetMultMax() const { return fMultMax; }
+
+  ClassDef(mfCentralityInfo, 0);
+};
+
+#endif

+ 18 - 0
eventanalysis/mfEventInfo.cxx

@@ -0,0 +1,18 @@
+#include "mfEventInfo.h"
+
+ClassImp(mfEventInfo);
+
+mfEventInfo::mfEventInfo(/* args */)
+{
+}
+
+mfEventInfo::~mfEventInfo()
+{
+  fQvectors.clear();
+}
+
+void mfEventInfo::SetCentralityRange(Float_t centMin, Float_t centMax)
+{
+  fCentRange.first = centMin;
+  fCentRange.second = centMax;
+}

+ 42 - 0
eventanalysis/mfEventInfo.h

@@ -0,0 +1,42 @@
+#ifndef MPD_FLOW_EVENTINFO_H
+#define MPD_FLOW_EVENTINFO_H
+
+#include <iostream>
+#include <vector>
+#include <Rtypes.h>
+#include <TString.h>
+
+#include "mfQvector.h"
+
+class mfEventInfo
+{
+private:
+  std::pair<Float_t,Float_t> fCentRange;
+  Int_t fMult;
+  std::vector<mfQvector> fQvectors;
+  std::vector<int> fEpMult;
+public:
+  mfEventInfo(/* args */);
+  virtual ~mfEventInfo();
+
+  // Setters
+  void SetCentralityRange(Float_t centMin, Float_t centMax);
+  void AddQv(mfQvector qv) { fQvectors.push_back(qv); }
+  void AddEpMult(int _a) { fEpMult.push_back(_a); }
+  void SetQv(std::vector<mfQvector> _a) { fQvectors = _a; }
+  void SetEpMult(std::vector<int> _a) { fEpMult = _a; }
+  void SetMultiplicity(Int_t _a) { fMult = _a; }
+
+  // Getters
+  std::pair<Float_t,Float_t> GetCentralityRange() const { return fCentRange; }
+  Float_t GetCentralityMean() const { return 0.5*(fCentRange.second - fCentRange.first); }
+  std::vector<mfQvector> GetQvectors() const { return fQvectors; }
+  std::vector<int> GetEpMult() const { return fEpMult; }
+  mfQvector GetQvector(Int_t i) const { return fQvectors.at(i); }
+  int GetEpMult(int i) const { return fEpMult.at(i); }
+  Int_t GetMultiplicity() const { return fMult; }
+
+  ClassDef(mfEventInfo,0);
+};
+
+#endif

+ 174 - 0
eventanalysis/mfProcessEvent.cxx

@@ -0,0 +1,174 @@
+#include <map>
+#include <iostream>
+
+#include <TString.h>
+#include <TMath.h>
+
+#include "mfProcessEvent.h"
+#include "mfQvector.h"
+
+using mfcent::CentB;
+using mfcent::isRefMult;
+using mfdca::NsigmaDcaCut;
+using mfEp::EpNames;
+using mfEp::GetFHCalPhi;
+using mfEp::NetaGaps;
+using mfEp::Nmodules;
+using mfEp::NsubEvents;
+using mfEp::NsubEventsTpc;
+using mfEp::TpcEtaGap;
+using mfEp::TpcLEp;
+using mfEp::TpcREp;
+using mfEp::EpNames::kFHCal;
+using mfEp::EpNames::kFHCal_L;
+using mfEp::EpNames::kFHCal_R;
+
+using TMath::Abs;
+using TMath::Cos;
+using TMath::Sin;
+
+ClassImp(mfProcessEvent);
+
+mfProcessEvent::mfProcessEvent(/* args */)
+{
+}
+
+mfProcessEvent::~mfProcessEvent()
+{
+}
+
+Float_t mfProcessEvent::GetCentB(PicoDstMCEvent *const &mcEvent)
+{
+  return CentB(mcEvent->GetB());
+}
+
+mfEventInfo mfProcessEvent::ProcessRecoEP(TClonesArray *const &recoTracks, TClonesArray *const &fhcalModules)
+{
+  int mult = 0;
+  std::vector<mfQvector> Qvectors;
+  std::vector<double> Qx, Qy, Wt;
+  std::vector<int> EpMult;
+
+  if (!recoTracks)
+    return fEventInfo;
+  if (recoTracks->GetEntriesFast() == 0)
+    return fEventInfo;
+  if (!fInputReader)
+    return fEventInfo;
+
+  for (int i = 0; i < NsubEvents; i++)
+  {
+    Qx.push_back(0.);
+    Qy.push_back(0.);
+    Wt.push_back(0.);
+    EpMult.push_back(0);
+  }
+
+  for (int iTrack = 0; iTrack < recoTracks->GetEntriesFast(); iTrack++)
+  {
+    auto recoTrack = (PicoDstRecoTrack *)recoTracks->UncheckedAt(iTrack);
+    if (!recoTrack)
+      continue;
+    if (isRefMult(recoTrack->GetPt(), recoTrack->GetEta()))
+    {
+      mult++;
+    }
+
+    // Primary track selection
+    if (!fInputReader->isPrimary(recoTrack))
+      continue;
+
+    for (int iGap = 0; iGap < NetaGaps; iGap++)
+    {
+      // Tpc left
+      if (recoTrack->GetEta() < -1. * TpcEtaGap[iGap])
+      {
+        Qx.at(TpcLEp[iGap]) += recoTrack->GetPt() * Cos(2. * recoTrack->GetPhi());
+        Qy.at(TpcLEp[iGap]) += recoTrack->GetPt() * Sin(2. * recoTrack->GetPhi());
+        Wt.at(TpcLEp[iGap]) += recoTrack->GetPt();
+        EpMult.at(TpcLEp[iGap]) += 1;
+      }
+      // Tpc right
+      if (recoTrack->GetEta() > TpcEtaGap[iGap])
+      {
+        Qx.at(TpcREp[iGap]) += recoTrack->GetPt() * Cos(2. * recoTrack->GetPhi());
+        Qy.at(TpcREp[iGap]) += recoTrack->GetPt() * Sin(2. * recoTrack->GetPhi());
+        Wt.at(TpcREp[iGap]) += recoTrack->GetPt();
+        EpMult.at(TpcREp[iGap]) += 1;
+      }
+    }
+  }
+  for (int i = 0; i < NsubEventsTpc; i++)
+  {
+    if (Wt.at(i) == 0 || EpMult.at(i) == 0)
+    {
+      // std::cerr << "No tracks for EP" << std::endl;
+      return fEventInfo;
+    }
+    else
+    {
+      Qx.at(i) /= Wt.at(i);
+      Qy.at(i) /= Wt.at(i);
+      Qvectors.push_back(mfQvector(2., Qx.at(i), Qy.at(i)));
+      // std::cout << "Qx[" << i << "] = " << Qx.at(i) << ", "
+      //           << "Qy[" << i << "] = " << Qy.at(i) << std::endl;
+    }
+  }
+
+  for (int iMod = 0; iMod < fhcalModules->GetEntriesFast(); iMod++)
+  {
+    auto fhcalModule = (PicoDstFHCal *)fhcalModules->UncheckedAt(iMod);
+    if (!fhcalModule)
+      continue;
+    if ((iMod == 22) || (iMod == 67))
+      continue;
+    double phi = GetFHCalPhi(iMod);
+
+    // std::cout << "\tModule " << iMod
+    //           << ": phi " << phi
+    //           << ", Energy " << fhcalModule->GetEnergy() << std::endl;
+    
+    int w_sign = (iMod < Nmodules) ? 1 : -1;
+    // for left FHCal modules
+    if (iMod < Nmodules)
+    {
+      Qx.at(kFHCal_L) += w_sign * fhcalModule->GetEnergy() * Cos(1. * phi);
+      Qy.at(kFHCal_L) += w_sign * fhcalModule->GetEnergy() * Sin(1. * phi);
+      Wt.at(kFHCal_L) += fhcalModule->GetEnergy();
+    }
+    // for right FHCal modules
+    if (iMod >= Nmodules)
+    {
+      Qx.at(kFHCal_R) += w_sign * fhcalModule->GetEnergy() * Cos(1. * phi);
+      Qy.at(kFHCal_R) += w_sign * fhcalModule->GetEnergy() * Sin(1. * phi);
+      Wt.at(kFHCal_R) += fhcalModule->GetEnergy();
+    }
+    // Total FHCal modules
+    Qx.at(kFHCal) += w_sign * fhcalModule->GetEnergy() * Cos(1. * phi);
+    Qy.at(kFHCal) += w_sign * fhcalModule->GetEnergy() * Sin(1. * phi);
+    Wt.at(kFHCal) += fhcalModule->GetEnergy();
+  }
+
+  for (int i = NsubEventsTpc; i < NsubEvents; i++)
+  {
+    if (Wt.at(i) == 0)
+    {
+      // std::cerr << "No tracks for EP" << std::endl;
+      return fEventInfo;
+    }
+    else
+    {
+      Qx.at(i) /= Wt.at(i);
+      Qy.at(i) /= Wt.at(i);
+      Qvectors.push_back(mfQvector(1., Qx.at(i), Qy.at(i)));
+      // std::cout << "Qx[" << i << "] = " << Qx.at(i) << ", "
+      //           << "Qy[" << i << "] = " << Qy.at(i) << std::endl;
+    }
+  }
+
+  fEventInfo.SetMultiplicity(mult);
+  fEventInfo.SetQv(Qvectors);
+  fEventInfo.SetEpMult(EpMult);
+
+  return fEventInfo;
+}

+ 36 - 0
eventanalysis/mfProcessEvent.h

@@ -0,0 +1,36 @@
+#ifndef MPD_FLOW_PROCESSEVENT_H
+#define MPD_FLOW_PROCESSEVENT_H
+
+#include <Rtypes.h>
+#include <TString.h>
+#include <TClonesArray.h>
+
+#include <PicoDstMCEvent.h>
+#include <PicoDstRecoEvent.h>
+#include <PicoDstMCTrack.h>
+#include <PicoDstRecoTrack.h>
+#include <PicoDstFHCal.h>
+
+#include "mfEventInfo.h"
+#include "mfNamespace.h"
+#include "mfInputReader.h"
+
+class mfProcessEvent
+{
+private:
+  mfEventInfo fEventInfo;
+  mfInputReader *fInputReader;
+public:
+  mfProcessEvent(/* args */);
+  virtual ~mfProcessEvent();
+
+  void SetInputReader(mfInputReader *const& _a) { fInputReader = _a; }
+
+  Float_t GetCentB(PicoDstMCEvent *const& mcEvent);
+  mfEventInfo ProcessRecoEP(TClonesArray *const& recoTracks, TClonesArray *const& fhcalModules);
+  Bool_t isPrimary(PicoDstRecoTrack *const& recoTrack);
+
+  ClassDef(mfProcessEvent,0);
+};
+
+#endif

+ 18 - 0
qvector/mfQvector.cxx

@@ -0,0 +1,18 @@
+#include "mfQvector.h"
+
+ClassImp(mfQvector);
+
+mfQvector::mfQvector(/* args */)
+{
+}
+
+mfQvector::mfQvector(Float_t harm, Float_t x, Float_t y)
+{
+  fHarmonic = harm;
+  fX = x;
+  fY = y;
+}
+
+mfQvector::~mfQvector()
+{
+}

+ 26 - 0
qvector/mfQvector.h

@@ -0,0 +1,26 @@
+#ifndef MPD_FLOW_QVECTOR_H
+#define MPD_FLOW_QVECTOR_H
+
+#include <Rtypes.h>
+#include <TVector2.h>
+
+class mfQvector : public TVector2
+{
+private:
+  Float_t fHarmonic;
+public:
+  mfQvector(/* args */);
+  mfQvector(Float_t harm, Float_t x, Float_t y);
+  virtual ~mfQvector();
+
+  // Setters
+  void SetHarmonic(Float_t _a) { fHarmonic = _a; }
+
+  // Getters
+  Float_t GetHarmonic() const { return fHarmonic; }
+  Float_t GetPhi() const { return (this->Phi()/fHarmonic); }
+
+  ClassDef(mfQvector,0);
+};
+
+#endif

+ 43 - 0
utils/mfDcaReader.cxx

@@ -0,0 +1,43 @@
+#include "mfDcaReader.h"
+#include <TFile.h>
+#include <TMath.h>
+
+using mfdca::Ndim;
+using mfdca::NsigmaDcaCut;
+
+using TMath::Abs;
+
+ClassImp(mfDcaReader);
+
+mfDcaReader::mfDcaReader(/* args */)
+{
+}
+
+mfDcaReader::~mfDcaReader()
+{
+}
+
+void mfDcaReader::SetDcaFitFile(TString str)
+{
+  fDcaFileName = str;
+  TFile *fi = new TFile(fDcaFileName.Data(),"READ");
+  for (int i = 0; i < Ndim; i++)
+  {
+    fdca[i] = (TF2 *)fi->Get(Form("f_sigma%i", i));
+  }
+  fi->Close();
+}
+
+Bool_t mfDcaReader::isPrimary(PicoDstRecoTrack *recoTrack)
+{
+  if (fDcaFileName == "") return false;
+  if (!recoTrack) return false;
+  if (Abs(recoTrack->GetDCAx()) > fdca[0]->Eval(recoTrack->GetPt(),recoTrack->GetEta())*NsigmaDcaCut)
+  return false;
+  if (Abs(recoTrack->GetDCAy()) > fdca[1]->Eval(recoTrack->GetPt(),recoTrack->GetEta())*NsigmaDcaCut)
+  return false;
+  if (Abs(recoTrack->GetDCAz()) > fdca[2]->Eval(recoTrack->GetPt(),recoTrack->GetEta())*NsigmaDcaCut)
+  return false;
+
+  return true;
+}

+ 27 - 0
utils/mfDcaReader.h

@@ -0,0 +1,27 @@
+#ifndef MPD_FLOW_DCAEADER_H
+#define MPD_FLOW_DCAEADER_H
+
+#include <TF2.h>
+#include <TString.h>
+#include "mfNamespace.h"
+
+#include <PicoDstRecoTrack.h>
+
+using mfdca::Ndim;
+
+class mfDcaReader
+{
+private:
+  TString fDcaFileName;
+  TF2 *fdca[Ndim];
+public:
+  mfDcaReader(/* args */);
+  virtual ~mfDcaReader();
+
+  void SetDcaFitFile(TString str);
+  Bool_t isPrimary(PicoDstRecoTrack *recoTrack);
+
+  ClassDef(mfDcaReader,0);
+};
+
+#endif

+ 11 - 0
utils/mfInputReader.cxx

@@ -0,0 +1,11 @@
+#include "mfInputReader.h"
+
+ClassImp(mfInputReader);
+
+mfInputReader::mfInputReader(/* args */)
+{
+}
+
+mfInputReader::~mfInputReader()
+{
+}

+ 21 - 0
utils/mfInputReader.h

@@ -0,0 +1,21 @@
+#ifndef MPD_FLOW_INPUTREADER_H
+#define MPD_FLOW_INPUTREADER_H
+
+#include <TF2.h>
+#include <TString.h>
+#include "mfNamespace.h"
+#include "mfDcaReader.h"
+
+
+class mfInputReader : public mfDcaReader
+{
+private:
+  /*args*/
+public:
+  mfInputReader(/* args */);
+  virtual ~mfInputReader();
+
+  ClassDef(mfInputReader,0);
+};
+
+#endif

+ 80 - 2
utils/mfNamespace.cxx

@@ -1,6 +1,12 @@
 #include "mfNamespace.h"
+#include <TMath.h>
+#include <TVector2.h>
 
-int mpdflow::GetPtBin(float pt)
+using TMath::ATan2;
+using TMath::Cos;
+using TMath::Sin;
+
+int mfdca::GetPtBin(float pt)
 {
   int pt_bin = -1;
   for (int i = 0; i < NptBins; ++i)
@@ -9,11 +15,83 @@ int mpdflow::GetPtBin(float pt)
   return pt_bin;
 }
 
-int mpdflow::GetEtaBin(float eta)
+int mfdca::GetEtaBin(float 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;
+}
+
+float mfcent::CentB(float bimp)
+{
+  // Hard coded centrality defenition
+  // based on the impact parameter
+  float fcent;
+  if (bimp < 4.18)
+    fcent = 0; // 0-10%
+  else if (bimp < 6.01)
+    fcent = 10; //10-20%
+  else if (bimp < 7.37)
+    fcent = 20; //20-30%
+  else if (bimp < 8.52)
+    fcent = 30; //30-40%
+  else if (bimp < 9.57)
+    fcent = 40; //40-50%
+  else if (bimp < 10.55)
+    fcent = 50; //50-60%
+  else if (bimp < 11.46)
+    fcent = 60; //60-70%
+  else if (bimp < 12.31)
+    fcent = 70; //70-80%
+  else
+    fcent = -1;
+
+  return fcent + 5.;
+}
+
+bool mfcent::isRefMult(float pt, float eta)
+{
+  bool result = (pt > ptRefMultCut && abs(eta) < etaRefMultCut) ? true : false;
+  return result;
+}
+
+double mfEp::GetFHCalPhi(int iModule)
+{
+  int xAxisSwitch = (iModule < Nmodules) ? 1 : -1;
+  int module = (iModule < Nmodules) ? iModule : iModule - Nmodules;
+  double x, y, phi=-999.;
+  if (module >= 0 && module <= 4)
+  {
+    y = 45.;
+    x = (module - 2) * 15.;
+    phi = ATan2(y, x * xAxisSwitch);
+  }
+  else if ((module >= 5) && (module <= 39))
+  {
+    y = (3-(module+2)/7)*15.;
+    x = (3-(module+2)%7)*15.;
+    phi = ATan2(y,x*xAxisSwitch);
+  }
+  else if ((module>=40) && (module<=44))
+  {
+    y = -45.;
+    x = (module-42)*15.;
+    phi = ATan2(y,x*xAxisSwitch);
+  }
+
+  return phi;
+}
+
+float mfEp::RangeAngle_02Pi(float phi)
+{
+  TVector2 vect;
+  vect.Set(Cos(phi),Sin(phi));
+  return vect.Phi();
+}
+
+float mfEp::RangeAngle(float phi)
+{
+  return ATan2(Sin(phi),Cos(phi));
 }

+ 47 - 2
utils/mfNamespace.h

@@ -1,7 +1,9 @@
 #ifndef MPD_FLOW_NAMESPACE_H
 #define MPD_FLOW_NAMESPACE_H
 
-namespace mpdflow
+#include <TString.h>
+
+namespace mfdca
 {
   const float ptBins[] = {0., 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.};
   const int NptBins = 12;
@@ -13,8 +15,51 @@ namespace mpdflow
 
   const int NdcaFitIter = 3;
 
+  const float NsigmaDcaCut = 2.;
+
   int GetPtBin(float pt);
   int GetEtaBin(float eta);
-}; // namespace mpdflow
+}; // namespace mfdca
+
+namespace mfcent
+{
+  const float ptRefMultCut = 0.15;
+  const float etaRefMultCut = 0.5;
+
+  float CentB(float bimp);
+  bool isRefMult(float pt, float eta);
+}; // namespace mfcent
+
+namespace mfEp
+{
+  const int NetaGaps = 3;
+  enum EpNames
+  {
+    kTpc_R_EtaGap1,
+    kTpc_R_EtaGap2,
+    kTpc_R_EtaGap3,
+    kTpc_L_EtaGap1,
+    kTpc_L_EtaGap2,
+    kTpc_L_EtaGap3,
+    kFHCal_R,
+    kFHCal_L,
+    kFHCal
+  };
+  const int NsubEvents = kFHCal + 1;
+  const int NsubEventsTpc = kTpc_L_EtaGap3 + 1;
+  const float TpcEtaGap[NetaGaps] = {0.05, 0.1, 0.25};
+  const int TpcREp[NetaGaps] = {kTpc_R_EtaGap1, kTpc_R_EtaGap2, kTpc_R_EtaGap3};
+  const int TpcLEp[NetaGaps] = {kTpc_L_EtaGap1, kTpc_L_EtaGap2, kTpc_L_EtaGap3};
+
+  const int Nmodules = 45;
+  double GetFHCalPhi(int iModule);
+  float RangeAngle_02Pi(float phi); // make it from 0 to 2*Pi
+  float RangeAngle(float phi); // make it from -Pi to Pi
+}; // namespace mfTpcEp
+
+namespace mfRes
+{
+  const int NresBin = 10;
+}; // namespace mfRes
 
 #endif