Просмотр исходного кода

Created dedicated class for flattening corrections

PeterParfenov лет назад: 3
Родитель
Сommit
471de53a22
7 измененных файлов с 203 добавлено и 224 удалено
  1. 1 0
      MpdFlow.LinkDef.h
  2. 13 75
      bin/get-flow.cpp
  3. 25 55
      bin/get-qv.cpp
  4. 13 94
      bin/get-res.cpp
  5. 116 0
      eventanalysis/mfFlattening.cxx
  6. 30 0
      eventanalysis/mfFlattening.h
  7. 5 0
      eventanalysis/mfRecentering.cxx

+ 1 - 0
MpdFlow.LinkDef.h

@@ -34,6 +34,7 @@
 #pragma link C++ class mfQvector+;
 #pragma link C++ class mfEventInfo+;
 #pragma link C++ class mfRecentering+;
+#pragma link C++ class mfFlattening+;
 #pragma link C++ class mfDcaReader+;
 #pragma link C++ class mfInputReader+;
 #pragma link C++ class mfRecoTrackSelector+;

+ 13 - 75
bin/get-flow.cpp

@@ -4,6 +4,7 @@
 #include "mfInputReader.h"
 #include "mfSelector.h"
 #include "mfRecentering.h"
+#include "mfFlattening.h"
 
 #include <TStopwatch.h>
 #include <TProfile2D.h>
@@ -144,49 +145,8 @@ int main(int argc, char **argv)
   recentering->ReadRecenteringFile(recenteringFileName);
 
   // Read flattening
-  TProfile *tmp;
-  std::map<double, double> vCosFlat[NsubEvents][Nflattening];
-  std::map<double, double> vSinFlat[NsubEvents][Nflattening];
-  TFile *fiFlat;
-  if (flatteningFileName != "")
-  {
-    fiFlat = new TFile(flatteningFileName.Data(), "read");
-    fiFlat->cd();
-    for (int i = 0; i < NsubEvents; i++)
-    {
-      for (int j = 0; j < Nflattening; j++)
-      {
-        tmp = (TProfile *)fiFlat->Get(Form("pCosAvg_%i_%i", i, j));
-        for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
-        {
-          vCosFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
-        }
-        delete tmp;
-
-        tmp = (TProfile *)fiFlat->Get(Form("pSinAvg_%i_%i", i, j));
-        for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
-        {
-          vSinFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
-        }
-        delete tmp;
-      }
-    }
-    fiFlat->Close();
-  }
-  else
-  {
-    for (int i = 0; i < NsubEvents; i++)
-    {
-      for (int j = 0; j < Nflattening; j++)
-      {
-        for (int ibin = 0; ibin < NresBin; ibin++)
-        {
-          vCosFlat[i][j][ibin * 100. / NresBin+5.] = 0.;
-          vSinFlat[i][j][ibin * 100. / NresBin+5.] = 0.;
-        }
-      }
-    }
-  }
+  mfFlattening *flattening = new mfFlattening();
+  flattening->ReadFlatteningFile(flatteningFileName);
 
   // Init output profiles
   TProfile2D *pv2PtTpc[NetaGaps][Npid];
@@ -326,16 +286,7 @@ int main(int argc, char **argv)
       continue;
     procEvent->SetInputReader(inputReader);
 
-    // std::vector<double> vResTpc;
-    // for (int iGap=0; iGap<NetaGaps; iGap++)
-    // {
-    //   vResTpc.push_back(fResTpc[iGap]->Eval(cent));
-    // }
-
-    // std::cout << "\tNtracks = " << recoTracks->GetEntriesFast() << std::endl;
-
     mfQvector qv;
-    double dPsi = 0.;
 
     mfEventInfo evInfoTpcRaw = procEvent->ProcessRecoTpcEP(recoTracks, 2.);
     vPsiTpcR.clear();
@@ -343,32 +294,23 @@ int main(int argc, char **argv)
     PsiFHCal = 0.;
     if (evInfoTpcRaw.GetQvectors().size() == 0)
       continue;
-    mfEventInfo evInfoTpc = recentering->ProcessRecentering(evInfoTpcRaw, cent);
+    mfEventInfo evInfoTpcRec = recentering->ProcessRecentering(evInfoTpcRaw, cent);
+    mfEventInfo evInfoTpc = flattening->ProcessFlattening(evInfoTpcRec, cent);
 
     for (int iGap = 0; iGap < NetaGaps; iGap++)
     {
       // Ntracks >= 4 for TpcR and TpcL EP
       if (selector->isGoodTpcEp(evInfoTpc.GetEpMult(TpcREp[iGap]), evInfoTpc.GetEpMult(TpcLEp[iGap])))
       {
-        qv.SetHarmonic(2.);
+        qv.SetHarmonic(evInfoTpc.GetQvector(TpcREp[iGap]).GetHarmonic());
         qv.Set(evInfoTpc.GetQvector(TpcREp[iGap]).X(),
                evInfoTpc.GetQvector(TpcREp[iGap]).Y());
-        dPsi = 0.;
-        for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-        {
-          dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[TpcREp[iGap]][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[TpcREp[iGap]][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-        }
-        vPsiTpcR.push_back(qv.GetPhi() + dPsi);
+        vPsiTpcR.push_back(qv.GetPhi());
 
-        qv.SetHarmonic(2.);
+        qv.SetHarmonic(evInfoTpc.GetQvector(TpcLEp[iGap]).GetHarmonic());
         qv.Set(evInfoTpc.GetQvector(TpcLEp[iGap]).X(),
                evInfoTpc.GetQvector(TpcLEp[iGap]).Y());
-        dPsi = 0.;
-        for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-        {
-          dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[TpcLEp[iGap]][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[TpcLEp[iGap]][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-        }
-        vPsiTpcL.push_back(qv.GetPhi() + dPsi);
+        vPsiTpcL.push_back(qv.GetPhi());
       }
       else
       {
@@ -378,19 +320,15 @@ int main(int argc, char **argv)
     }
 
     mfEventInfo evInfoFHCalRaw = procEvent->ProcessRecoFHCalEP(fhcalmodules, 1.);
-    mfEventInfo evInfoFHCal = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
+    mfEventInfo evInfoFHCalRec = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
+    mfEventInfo evInfoFHCal = flattening->ProcessFlattening(evInfoFHCalRec, cent);
     // FHCal energy > 0 in FHCal_L and FHCal_R
     if (selector->isGoodFHCalEp(evInfoFHCal.GetEpWeight(EpNames::kFHCal_R), evInfoFHCal.GetEpWeight(EpNames::kFHCal_R)))
     {
-      qv.SetHarmonic(1.);
+      qv.SetHarmonic(evInfoFHCal.GetQvector(EpNames::kFHCal).GetHarmonic());
       qv.Set(evInfoFHCal.GetQvector(EpNames::kFHCal).X(),
              evInfoFHCal.GetQvector(EpNames::kFHCal).Y());
-      dPsi = 0.;
-      for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-      {
-        dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[EpNames::kFHCal][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[EpNames::kFHCal][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-      }
-      PsiFHCal = qv.GetPhi() + dPsi;
+      PsiFHCal = qv.GetPhi();
     }
     else
       PsiFHCal = -999.;

+ 25 - 55
bin/get-qv.cpp

@@ -3,6 +3,7 @@
 #include "mfProcessEvent.h"
 #include "mfInputReader.h"
 #include "mfRecentering.h"
+#include "mfFlattening.h"
 
 #include <vector>
 #include <map>
@@ -129,31 +130,8 @@ int main(int argc, char **argv)
   recentering->ReadRecenteringFile(recenteringFileName);
 
   // Read flattening
-  TProfile *tmp;
-  std::map<double, double> vCosFlat[NsubEvents][Nflattening];
-  std::map<double, double> vSinFlat[NsubEvents][Nflattening];
-  TFile *fiFlat = new TFile(flatteningFileName.Data(), "read");
-  fiFlat->cd();
-  for (int i = 0; i < NsubEvents; i++)
-  {
-    for (int j = 0; j < Nflattening; j++)
-    {
-      tmp = (TProfile *)fiFlat->Get(Form("pCosAvg_%i_%i", i, j));
-      for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
-      {
-        vCosFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
-      }
-      delete tmp;
-
-      tmp = (TProfile *)fiFlat->Get(Form("pSinAvg_%i_%i", i, j));
-      for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
-      {
-        vSinFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
-      }
-      delete tmp;
-    }
-  }
-  fiFlat->Close();
+  mfFlattening *flattening = new mfFlattening();
+  flattening->ReadFlatteningFile(flatteningFileName);
 
   // Init output profiles
   TH2F *hQxVsCent[3][NsubEvents];
@@ -231,7 +209,8 @@ int main(int argc, char **argv)
     mfEventInfo evInfoTpcRaw = procEvent->ProcessRecoTpcEP(recoTracks, 2.);
     if (evInfoTpcRaw.GetQvectors().size() == 0)
       continue;
-    mfEventInfo evInfoTpc = recentering->ProcessRecentering(evInfoTpcRaw, cent);
+    mfEventInfo evInfoTpcRec = recentering->ProcessRecentering(evInfoTpcRaw, cent);
+    mfEventInfo evInfoTpcFlat = flattening->ProcessFlattening(evInfoTpcRec, cent);
 
     std::vector<double> qx[3], qy[3], psi[3];
     std::vector<mfQvector> qv[3];
@@ -241,46 +220,37 @@ int main(int argc, char **argv)
 
       qx[0].push_back(evInfoTpcRaw.GetQvector(iSub).X());
       qy[0].push_back(evInfoTpcRaw.GetQvector(iSub).Y());
-      qx[1].push_back(evInfoTpc.GetQvector(iSub).X());
-      qy[1].push_back(evInfoTpc.GetQvector(iSub).Y());
-      qx[2].push_back(0.);
-      qy[2].push_back(0.);
-      qv[0].push_back(mfQvector(evInfoTpc.GetQvector(iSub).GetHarmonic(), qx[0].at(iSub), qy[0].at(iSub)));
-      qv[1].push_back(mfQvector(evInfoTpc.GetQvector(iSub).GetHarmonic(), qx[1].at(iSub), qy[1].at(iSub)));
-      qv[2].push_back(mfQvector(evInfoTpc.GetQvector(iSub).GetHarmonic(), qx[2].at(iSub), qy[2].at(iSub)));
+      qx[1].push_back(evInfoTpcRec.GetQvector(iSub).X());
+      qy[1].push_back(evInfoTpcRec.GetQvector(iSub).Y());
+      qx[2].push_back(evInfoTpcFlat.GetQvector(iSub).X());
+      qy[2].push_back(evInfoTpcFlat.GetQvector(iSub).Y());
+      qv[0].push_back(mfQvector(evInfoTpcRaw.GetQvector(iSub).GetHarmonic(), qx[0].at(iSub), qy[0].at(iSub)));
+      qv[1].push_back(mfQvector(evInfoTpcRec.GetQvector(iSub).GetHarmonic(), qx[1].at(iSub), qy[1].at(iSub)));
+      qv[2].push_back(mfQvector(evInfoTpcFlat.GetQvector(iSub).GetHarmonic(), qx[2].at(iSub), qy[2].at(iSub)));
       psi[0].push_back(qv[0].at(iSub).GetPhi());
       psi[1].push_back(qv[1].at(iSub).GetPhi());
-      double dPsi = 0.;
-      for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-      {
-        dPsi += 2 / ((iFlat + 1.) * qv[1].at(iSub).GetHarmonic()) * (-vSinFlat[iSub][iFlat].at(cent) * Cos((iFlat + 1.) * qv[1].at(iSub).GetHarmonic() * qv[1].at(iSub).GetPhi()) + vCosFlat[iSub][iFlat].at(cent) * Sin((iFlat + 1.) * qv[1].at(iSub).GetHarmonic() * qv[1].at(iSub).GetPhi()));
-      }
-      psi[2].push_back(qv[1].at(iSub).GetPhi() + dPsi);
+      psi[2].push_back(qv[2].at(iSub).GetPhi());
     }
 
     mfEventInfo evInfoFHCalRaw = procEvent->ProcessRecoFHCalEP(fhcalmodules, 1.);
-    mfEventInfo evInfoFHCal = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
+    mfEventInfo evInfoFHCalRec = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
+    mfEventInfo evInfoFHCalFlat = flattening->ProcessFlattening(evInfoFHCalRec, cent);
     for (int iSub = NsubEventsTpc; iSub < NsubEvents; iSub++)
     {
       // std::cout << "iSub = " << iSub << " cent = " << cent << std::endl;
 
       qx[0].push_back(evInfoFHCalRaw.GetQvector(iSub).X());
       qy[0].push_back(evInfoFHCalRaw.GetQvector(iSub).Y());
-      qx[1].push_back(evInfoFHCal.GetQvector(iSub).X());
-      qy[1].push_back(evInfoFHCal.GetQvector(iSub).Y());
-      qx[2].push_back(0.);
-      qy[2].push_back(0.);
-      qv[0].push_back(mfQvector(evInfoFHCal.GetQvector(iSub).GetHarmonic(), qx[0].at(iSub), qy[0].at(iSub)));
-      qv[1].push_back(mfQvector(evInfoFHCal.GetQvector(iSub).GetHarmonic(), qx[1].at(iSub), qy[1].at(iSub)));
-      qv[2].push_back(mfQvector(evInfoFHCal.GetQvector(iSub).GetHarmonic(), qx[2].at(iSub), qy[2].at(iSub)));
+      qx[1].push_back(evInfoFHCalRec.GetQvector(iSub).X());
+      qy[1].push_back(evInfoFHCalRec.GetQvector(iSub).Y());
+      qx[2].push_back(evInfoFHCalFlat.GetQvector(iSub).X());
+      qy[2].push_back(evInfoFHCalFlat.GetQvector(iSub).Y());
+      qv[0].push_back(mfQvector(evInfoFHCalRaw.GetQvector(iSub).GetHarmonic(), qx[0].at(iSub), qy[0].at(iSub)));
+      qv[1].push_back(mfQvector(evInfoFHCalRec.GetQvector(iSub).GetHarmonic(), qx[1].at(iSub), qy[1].at(iSub)));
+      qv[2].push_back(mfQvector(evInfoFHCalFlat.GetQvector(iSub).GetHarmonic(), qx[2].at(iSub), qy[2].at(iSub)));
       psi[0].push_back(qv[0].at(iSub).GetPhi());
       psi[1].push_back(qv[1].at(iSub).GetPhi());
-      double dPsi = 0.;
-      for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-      {
-        dPsi += 2 / ((iFlat + 1.) * qv[1].at(iSub).GetHarmonic()) * (-vSinFlat[iSub][iFlat].at(cent) * Cos((iFlat + 1.) * qv[1].at(iSub).GetHarmonic() * qv[1].at(iSub).GetPhi()) + vCosFlat[iSub][iFlat].at(cent) * Sin((iFlat + 1.) * qv[1].at(iSub).GetHarmonic() * qv[1].at(iSub).GetPhi()));
-      }
-      psi[2].push_back(qv[1].at(iSub).GetPhi() + dPsi);
+      psi[2].push_back(qv[2].at(iSub).GetPhi());
     }
 
     for (int i = 0; i < 3; i++)
@@ -288,10 +258,10 @@ int main(int argc, char **argv)
       for (int iGap = 0; iGap < NetaGaps; iGap++)
       {
         // Ntracks >= 2 for TpcR EP
-        if (evInfoTpc.GetEpMult(TpcREp[iGap]) < 2)
+        if (evInfoTpcFlat.GetEpMult(TpcREp[iGap]) < 2)
           continue;
         // Ntracks >= 2 for TpcL EP
-        if (evInfoTpc.GetEpMult(TpcLEp[iGap]) < 2)
+        if (evInfoTpcFlat.GetEpMult(TpcLEp[iGap]) < 2)
           continue;
 
         hQxVsCent[i][TpcREp[iGap]]->Fill(qv[i].at(TpcREp[iGap]).X(), cent);

+ 13 - 94
bin/get-res.cpp

@@ -4,6 +4,7 @@
 #include "mfInputReader.h"
 #include "mfSelector.h"
 #include "mfRecentering.h"
+#include "mfFlattening.h"
 
 #include <TStopwatch.h>
 #include <TProfile.h>
@@ -125,49 +126,8 @@ int main(int argc, char **argv)
   recentering->ReadRecenteringFile(recenteringFileName);
 
   // Read flattening
-  TProfile *tmp;
-  std::map<double, double> vCosFlat[NsubEvents][Nflattening];
-  std::map<double, double> vSinFlat[NsubEvents][Nflattening];
-  TFile *fiFlat;
-  if (flatteningFileName != "")
-  {
-    fiFlat = new TFile(flatteningFileName.Data(), "read");
-    fiFlat->cd();
-    for (int i = 0; i < NsubEvents; i++)
-    {
-      for (int j = 0; j < Nflattening; j++)
-      {
-        tmp = (TProfile *)fiFlat->Get(Form("pCosAvg_%i_%i", i, j));
-        for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
-        {
-          vCosFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
-        }
-        delete tmp;
-
-        tmp = (TProfile *)fiFlat->Get(Form("pSinAvg_%i_%i", i, j));
-        for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
-        {
-          vSinFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
-        }
-        delete tmp;
-      }
-    }
-    fiFlat->Close();
-  }
-  else
-  {
-    for (int i = 0; i < NsubEvents; i++)
-    {
-      for (int j = 0; j < Nflattening; j++)
-      {
-        for (int ibin = 0; ibin < NresBin; ibin++)
-        {
-          vCosFlat[i][j][ibin * 100. / NresBin+5.] = 0.;
-          vSinFlat[i][j][ibin * 100. / NresBin+5.] = 0.;
-        }
-      }
-    }
-  }
+  mfFlattening *flattening = new mfFlattening();
+  flattening->ReadFlatteningFile(flatteningFileName);
 
   // Init output profiles
   TProfile *pResRecoSqrTpc[NetaGaps];
@@ -270,35 +230,19 @@ int main(int argc, char **argv)
     mfEventInfo evInfoTpcRaw = procEvent->ProcessRecoTpcEP(recoTracks, 2.);
     if (evInfoTpcRaw.GetQvectors().size() == 0)
       continue;
-    mfEventInfo evInfoTpc = recentering->ProcessRecentering(evInfoTpcRaw, cent);
+    mfEventInfo evInfoTpcRec = recentering->ProcessRecentering(evInfoTpcRaw, cent);
+    mfEventInfo evInfoTpcFlat = flattening->ProcessFlattening(evInfoTpcRec, cent);
 
     mfQvector qv;
     double dPsi = 0.;
     for (int iGap = 0; iGap < NetaGaps; iGap++)
     {
       // Ntracks >= 4 for TpcR and TpcL EP
-      if (!selector->isGoodTpcEp(evInfoTpc.GetEpMult(TpcREp[iGap]), evInfoTpc.GetEpMult(TpcLEp[iGap])))
+      if (!selector->isGoodTpcEp(evInfoTpcRec.GetEpMult(TpcREp[iGap]), evInfoTpcRec.GetEpMult(TpcLEp[iGap])))
         continue;
 
-      qv.SetHarmonic(2.);
-      qv.Set(evInfoTpc.GetQvector(TpcREp[iGap]).X(),
-             evInfoTpc.GetQvector(TpcREp[iGap]).Y());
-      dPsi = 0.;
-      for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-      {
-        dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[TpcREp[iGap]][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[TpcREp[iGap]][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-      }
-      double PsiR = qv.GetPhi() + dPsi;
-
-      qv.SetHarmonic(2.);
-      qv.Set(evInfoTpc.GetQvector(TpcLEp[iGap]).X(),
-             evInfoTpc.GetQvector(TpcLEp[iGap]).Y());
-      dPsi = 0.;
-      for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-      {
-        dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[TpcLEp[iGap]][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[TpcLEp[iGap]][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-      }
-      double PsiL = qv.GetPhi() + dPsi;
+      double PsiR = evInfoTpcFlat.GetQvector(TpcREp[iGap]).GetPhi();
+      double PsiL = evInfoTpcFlat.GetQvector(TpcLEp[iGap]).GetPhi();
 
       double resReco = Cos(2. * (PsiR - PsiL));
       double resRpR = Cos(2. * (PsiR - mcEvent->GetPhiRP()));
@@ -315,37 +259,12 @@ int main(int argc, char **argv)
     // FHCal energy > 0 in FHCal_L and FHCal_R
     if (!selector->isGoodFHCalEp(evInfoFHCalRaw.GetEpWeight(EpNames::kFHCal_R), evInfoFHCalRaw.GetEpWeight(EpNames::kFHCal_R)))
       continue;
-    mfEventInfo evInfoFHCal = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
-
-    qv.SetHarmonic(1.);
-    qv.Set(evInfoFHCal.GetQvector(EpNames::kFHCal_R).X(),
-           evInfoFHCal.GetQvector(EpNames::kFHCal_R).Y());
-    dPsi = 0.;
-    for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-    {
-      dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[EpNames::kFHCal_R][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[EpNames::kFHCal_R][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-    }
-    double PsiFHCalR = qv.GetPhi() + dPsi;
-
-    qv.SetHarmonic(1.);
-    qv.Set(evInfoFHCal.GetQvector(EpNames::kFHCal_L).X(),
-           evInfoFHCal.GetQvector(EpNames::kFHCal_L).Y());
-    dPsi = 0.;
-    for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-    {
-      dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[EpNames::kFHCal_L][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[EpNames::kFHCal_L][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-    }
-    double PsiFHCalL = qv.GetPhi() + dPsi;
+    mfEventInfo evInfoFHCalRec = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
+    mfEventInfo evInfoFHCalFlat = flattening->ProcessFlattening(evInfoFHCalRec, cent);
 
-    qv.SetHarmonic(1.);
-    qv.Set(evInfoFHCal.GetQvector(EpNames::kFHCal).X(),
-           evInfoFHCal.GetQvector(EpNames::kFHCal).Y());
-    dPsi = 0.;
-    for (int iFlat = 0; iFlat < Nflattening; iFlat++)
-    {
-      dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[EpNames::kFHCal][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[EpNames::kFHCal][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
-    }
-    double PsiFHCal = qv.GetPhi() + dPsi;
+    double PsiFHCalR = evInfoFHCalFlat.GetQvector(EpNames::kFHCal_R).GetPhi();
+    double PsiFHCalL = evInfoFHCalFlat.GetQvector(EpNames::kFHCal_L).GetPhi();
+    double PsiFHCal = evInfoFHCalFlat.GetQvector(EpNames::kFHCal).GetPhi();
 
     double resEpFHCal = Cos(1. * (PsiFHCalR - PsiFHCalL));
     double resV1RpFHCal = Cos(1. * (PsiFHCal - mcEvent->GetPhiRP()));

+ 116 - 0
eventanalysis/mfFlattening.cxx

@@ -0,0 +1,116 @@
+#include "mfFlattening.h"
+#include "mfQvector.h"
+
+#include <TProfile.h>
+#include <TFile.h>
+#include <TMath.h>
+
+using TMath::Sin;
+using TMath::Cos;
+
+using mfRes::NresBin;
+
+ClassImp(mfFlattening);
+
+mfFlattening::mfFlattening(/* args */)
+{
+}
+
+mfFlattening::~mfFlattening()
+{
+  for (int i = 0; i < NsubEvents; i++)
+  {
+    for (int j = 0; j < Nflattening; j++)
+    {
+      vCosFlat[i][j].clear();
+      vSinFlat[i][j].clear();
+    }
+  }
+}
+
+void mfFlattening::ReadFlatteningFile(TString str)
+{
+  TFile *fiFlat;
+  TProfile *tmp;
+  if (str != "")
+  {
+    fiFlat = new TFile(str.Data(), "read");
+    fiFlat->cd();
+    for (int i = 0; i < NsubEvents; i++)
+    {
+      for (int j = 0; j < Nflattening; j++)
+      {
+        tmp = (TProfile *)fiFlat->Get(Form("pCosAvg_%i_%i", i, j));
+        for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
+        {
+          vCosFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
+        }
+        delete tmp;
+
+        tmp = (TProfile *)fiFlat->Get(Form("pSinAvg_%i_%i", i, j));
+        for (int ibin = 0; ibin < tmp->GetNbinsX(); ibin++)
+        {
+          vSinFlat[i][j][tmp->GetBinCenter(ibin + 1)] = tmp->GetBinContent(ibin + 1);
+        }
+        delete tmp;
+      }
+    }
+    fiFlat->Close();
+  }
+  else
+  {
+    for (int i = 0; i < NsubEvents; i++)
+    {
+      for (int j = 0; j < Nflattening; j++)
+      {
+        for (int ibin = 0; ibin < NresBin; ibin++)
+        {
+          vCosFlat[i][j][ibin * 100. / NresBin + 5.] = 0.;
+          vSinFlat[i][j][ibin * 100. / NresBin + 5.] = 0.;
+        }
+      }
+    }
+  }
+}
+
+mfEventInfo mfFlattening::ProcessFlattening(mfEventInfo evt, float cent)
+{
+  mfEventInfo result;
+
+  result.SetCentralityRange(evt.GetCentralityRange().first, evt.GetCentralityRange().second);
+  // result.SetQv(evt.GetQvectors());
+  result.SetEpMult(evt.GetEpMult());
+  result.SetEpWeight(evt.GetEpWeight());
+  result.SetMultiplicity(evt.GetMultiplicity());
+
+  mfQvector qv;
+  double dPsi = 0., Psi = 0., Qmod = 0.;;
+
+  for (int iSub = 0; iSub < NsubEvents; iSub++)
+  {
+    // Recenter Q-vector if it's not empty ( i.e. Qv = (0.,0.) )
+    if (evt.GetQvector(iSub).X() != 0. && evt.GetQvector(iSub).Y() != 0.)
+    {
+      qv.SetHarmonic(evt.GetQvector(iSub).GetHarmonic());
+      qv.Set(evt.GetQvector(iSub).X(), evt.GetQvector(iSub).Y());
+      dPsi = 0.;
+      for (int iFlat = 0; iFlat < Nflattening; iFlat++)
+      {
+        dPsi += 2 / ((iFlat + 1.) * qv.GetHarmonic()) * (-vSinFlat[iSub][iFlat].at(cent) * Cos((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()) + vCosFlat[iSub][iFlat].at(cent) * Sin((iFlat + 1.) * qv.GetHarmonic() * qv.GetPhi()));
+      }
+      Psi = qv.GetPhi() + dPsi;
+      Qmod = qv.Mod();
+
+      qv.Set(0.,0.);
+      qv.Set(Qmod*Cos(qv.GetHarmonic()*Psi), Qmod*Sin(qv.GetHarmonic()*Psi));
+    }
+    else
+    {
+      qv.Set(evt.GetQvector(iSub).X(), evt.GetQvector(iSub).Y());
+      qv.SetHarmonic(evt.GetQvector(iSub).GetHarmonic());
+    }
+    result.AddQv(qv);
+  }
+
+  return result;
+}

+ 30 - 0
eventanalysis/mfFlattening.h

@@ -0,0 +1,30 @@
+#ifndef MPD_FLOW_FLATTENING_H
+#define MPD_FLOW_FLATTENING_H
+
+#include <map>
+#include <Rtypes.h>
+#include <TString.h>
+
+#include "mfNamespace.h"
+#include "mfEventInfo.h"
+
+using mfEp::Nflattening;
+using mfEp::NsubEvents;
+
+class mfFlattening
+{
+private:
+  std::map<double, double> vCosFlat[NsubEvents][Nflattening];
+  std::map<double, double> vSinFlat[NsubEvents][Nflattening];
+public:
+  mfFlattening(/* args */);
+  virtual ~mfFlattening();
+
+  void ReadFlatteningFile(TString str);
+
+  mfEventInfo ProcessFlattening(mfEventInfo evt, float cent);
+
+  ClassDef(mfFlattening, 0);
+};
+
+#endif

+ 5 - 0
eventanalysis/mfRecentering.cxx

@@ -13,6 +13,11 @@ mfRecentering::mfRecentering(/* args */)
 
 mfRecentering::~mfRecentering()
 {
+  for (int i = 0; i < NsubEvents; i++)
+  {
+    vQxRec[i].clear();
+    vQyRec[i].clear();
+  }
 }
 
 void mfRecentering::ReadRecenteringFile(TString str)