StFlowMephi/FemtoDstAnalyzer_FlowChargedHadrons.C

/**
 * \brief Example of how to read a file (list of files) using StFemtoEvent classes
 *
 * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
 * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
 * inFile.list) as an input, and preform physics analysis
 *
 * \author Grigory Nigmatkulov
 * \date May 29, 2018
 */

// This is needed for calling standalone classes
#define _VANILLA_ROOT_

// C++ headers
#include <iostream>
#include <vector>
#include <map>

// ROOT headers
#include "TROOT.h"
#include "TFile.h"
#include "TChain.h"
#include "TVector2.h"
#include "TVector3.h"
#include "TTree.h"
#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TProfile2D.h"
#include "TStopwatch.h"

// FemtoDst headers
#include "StFemtoDstReader.h"
#include "StFemtoDst.h"
#include "StFemtoEvent.h"
#include "StFemtoTrack.h"

// Load libraries (for ROOT_VERSTION_CODE >= 393215)
#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
R__LOAD_LIBRARY(/ mnt / pool / rhic / 4 / parfenovpeter / STAR / build / libStFemtoDst.so)
#endif

#include "Constants.h"

// inFile - is a name of name.FemtoDst.root file or a name
//          of a name.lis(t) files, that contains a list of
//          name1.FemtoDst.root, name2.FemtoDst.root, ... files

//Used functions (see them below)
Bool_t isGoodEvent(StFemtoEvent *const &event);
Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
Double_t GetWeight(StFemtoTrack *const &track);
TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
Double_t AngleShift(Double_t Psi, Double_t order);
Int_t GetVzBin(Double_t vtxZ);
Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
Bool_t isBadRun(StFemtoEvent *const &event);
Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip);

//_________________
void FemtoDstAnalyzer_FlowChargedHadrons(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
                                         const Char_t *outFile = "./oShiftTest.root",
                                         const Char_t *reCentFile = "",
                                         const Char_t *shiftFile = "",
                                         const Char_t *resFile = "",
                                         const Char_t *gainBBC = "",
                                         const Char_t *reCentFileBBC = "",
                                         const Char_t *shiftFileBBC = "",
                                         const Char_t *reCentFileZDC = "",
                                         const Char_t *shiftFileZDC = "",
                                         const Char_t *resFitFile = "")
{

  std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  TStopwatch timer;
  timer.Start();

#if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  gSystem->Load("../libStFemtoDst.so");
#endif

  StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  femtoReader->Init();

  // This is a way if you want to spead up IO
  std::cout << "Explicit read status for some branches" << std::endl;
  femtoReader->SetStatus("*", 0);
  femtoReader->SetStatus("Event", 1);
  femtoReader->SetStatus("Track", 1);
  std::cout << "Status has been set" << std::endl;

  std::cout << "Now I know what to read, Master!" << std::endl;

  if (!femtoReader->chain())
  {
    std::cout << "No chain has been found." << std::endl;
  }
  Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  std::cout << "eventsInTree: " << eventsInTree << std::endl;
  Long64_t events2read = femtoReader->chain()->GetEntries();

  std::cout << "Number of events to read: " << events2read << std::endl;

  //Read ReCentering <Q> vectors
  TFile *fiReCent = new TFile(reCentFile, "read");

  std::map<Double_t, Double_t> p_Qx_FULL_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_Qy_FULL_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_Qx_FULL_SP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_Qy_FULL_SP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_Qx_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qy_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qx_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qy_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qx_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qy_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qx_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Qy_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];

  TProfile2D *prof;
  double binCont;

  std::cout << "Reading file for TPC recentering..." << std::endl;
  for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
      prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_EP%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Qx_FULL_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_EP%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Qy_FULL_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_SP%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Qx_FULL_SP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_SP%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Qy_FULL_SP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;

      for (int iEtaGap = 0; iEtaGap < NEtaGaps; iEtaGap++)
      {
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qx_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qy_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qx_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qy_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;

        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qx_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qy_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qx_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Qy_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
      }
    }
  }

  //Read Shift coefficients
  TFile *fiShift = new TFile(shiftFile, "read");

  std::map<Double_t, Double_t> p_Cos2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Sin2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Cos2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Sin2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Cos2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Cos2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Cos2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Cos2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];

  std::map<Double_t, Double_t> p_Cos3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Sin3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Cos3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Sin3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  std::map<Double_t, Double_t> p_Cos3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Cos3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Cos3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Cos3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  std::map<Double_t, Double_t> p_Sin3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];

  std::cout << "Reading file for TPC flattening..." << std::endl;
  for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  {
    for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
    {
      std::cout << "iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
      prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Cos2_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Sin2_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Cos2_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Sin2_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;

      prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Cos3_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Sin3_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Cos3_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_Sin3_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;

      for (int iGap = 0; iGap < NEtaGaps; iGap++)
      {
        // for 2nd harmonics
        // EAST
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        // WEST
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;

        // for 3rd harmonics
        // EAST
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        // WEST
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Cos3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_Sin3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
      }
    }
  }

  //Read resolution
  TFile *fiRes = new TFile(resFile, "read");

  TProfile *p_res2_EP[NEtaGaps];
  TProfile *p_res3_EP[NEtaGaps];
  TProfile *p_res2_SP[NEtaGaps];
  TProfile *p_res3_SP[NEtaGaps];
  // TProfile *p_resBBC_EP[fNharmonics+1];

  TProfile *p_res2_EP_3sub_AB;
  TProfile *p_res2_EP_3sub_AC;
  TProfile *p_res2_EP_3sub_BC;

  p_res2_EP_3sub_AB = (TProfile *)fiRes->Get("p_res2_EP_3sub_AB");
  p_res2_EP_3sub_AC = (TProfile *)fiRes->Get("p_res2_EP_3sub_AC");
  p_res2_EP_3sub_BC = (TProfile *)fiRes->Get("p_res2_EP_3sub_BC");

  for (int iGap = 0; iGap < NEtaGaps; iGap++)
  {
    p_res2_EP[iGap] = (TProfile *)fiRes->Get(Form("p_res2_EP_gap%i", iGap));
    p_res3_EP[iGap] = (TProfile *)fiRes->Get(Form("p_res3_EP_gap%i", iGap));
    p_res2_SP[iGap] = (TProfile *)fiRes->Get(Form("p_res2_SP_gap%i", iGap));
    p_res3_SP[iGap] = (TProfile *)fiRes->Get(Form("p_res3_SP_gap%i", iGap));
  }

  // for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  // {
  //   p_resBBC_EP[iHarm] = (TProfile*) fiRes->Get(Form("p_resBBC_EP%i", iHarm));
  // }

  //Read resolution fit
  TFile *fiResFit = new TFile(resFitFile, "read");

  TH1D *hResBBCFULL[fNharmonics + 1];
  TH1D *hResBBCHalf[fNharmonics + 1];

  TH1D *hResZDCFULL[fNharmonics + 1];
  TH1D *hResZDCHalf[fNharmonics + 1];

  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    hResBBCFULL[iHarm] = (TH1D *)fiResFit->Get(Form("hResBBCFULL%i", iHarm));
    hResBBCHalf[iHarm] = (TH1D *)fiResFit->Get(Form("hResBBCHalf%i", iHarm));
  }

  for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  {
    hResZDCFULL[iHarm] = (TH1D *)fiResFit->Get(Form("hResZDCFULL%i", iHarm));
    hResZDCHalf[iHarm] = (TH1D *)fiResFit->Get(Form("hResZDCHalf%i", iHarm));
  }

  //Manage gain correction
  TProfile2D *BBCgain_East[fArraySize];
  TProfile2D *BBCgain_West[fArraySize];

  TFile *fiGain = new TFile(gainBBC, "read");
  fiGain->cd();

  for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  {
    BBCgain_East[iVtx] = (TProfile2D *)fiGain->Get(Form("p_gainBBC_East_Vtx%i", iVtx));
    BBCgain_West[iVtx] = (TProfile2D *)fiGain->Get(Form("p_gainBBC_West_Vtx%i", iVtx));
  }

  //Manage recentering BBC
  std::map<Double_t, Double_t> p_BBC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_BBC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_BBC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_BBC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_BBC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_BBC_Qy_West_EP[fNharmonics][fArraySize][fNCent];

  TFile *fiReCentBBC = new TFile(reCentFileBBC, "read");

  std::cout << "Reading file for BBC recentering..." << std::endl;
  for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
      prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_BBC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_BBC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;

      prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_BBC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_BBC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_BBC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_BBC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
    }
  }

  //Manage ShiftCorr BBC
  std::map<Double_t, Double_t> p_BBC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_BBC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_BBC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_BBC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_BBC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_BBC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];

  TFile *fiShiftBBC = new TFile(shiftFileBBC, "read");
  fiShiftBBC->cd();

  std::cout << "Reading file for BBC shift corrections..." << std::endl;
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      for (int iShift = 0; iShift < NShiftOrderMax * 3; iShift++)
      {
        std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
        prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_BBC_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_BBC_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;

        prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_BBC_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_BBC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_BBC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_BBC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
      }
    }
  }

  //Manage recentering ZDC
  std::map<Double_t, Double_t> p_ZDC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Qy_West_EP[fNharmonics][fArraySize][fNCent];

  TFile *fiReCentZDC = new TFile(reCentFileZDC, "read");

  std::cout << "Reading file for ZDC recentering..." << std::endl;
  for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
      prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_ZDC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_ZDC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;

      prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_ZDC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_ZDC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_ZDC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
      prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
      for (int i = 0; i < fNBins; i++)
      {
        for (int iCent = 0; iCent < fNCent; iCent++)
        {
          binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
          if (binCont == 0.)
            continue;
          p_ZDC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
        }
      }
      delete prof;
    }
  }

  //Manage ShiftCorr ZDC
  std::map<Double_t, Double_t> p_ZDC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  std::map<Double_t, Double_t> p_ZDC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];

  TFile *fiShiftZDC = new TFile(shiftFileZDC, "read");
  fiShiftZDC->cd();

  std::cout << "Reading file for ZDC shift corrections..." << std::endl;
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      for (int iShift = 0; iShift < NShiftOrderMax * 3; iShift++)
      {
        std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
        prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_ZDC_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_ZDC_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;

        prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_ZDC_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_ZDC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_ZDC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
        prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
        for (int i = 0; i < fNBins; i++)
        {
          for (int iCent = 0; iCent < fNCent; iCent++)
          {
            binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
            if (binCont == 0.)
              continue;
            p_ZDC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
          }
        }
        delete prof;
      }
    }
  }

  //Manage output TProfiles
  //TProfile *p_v2_EP[9][NEtaGaps];
  //TProfile *p_v3_EP[9][NEtaGaps];

  //TProfile *p_v2cent_EP[NEtaGaps];
  //TProfile *p_v3cent_EP[NEtaGaps];

  TProfile2D *p_v2cent_EP_ptgaps[NEtaGaps];
  TProfile2D *p_v3cent_EP_ptgaps[NEtaGaps];
  TProfile2D *p_v2cent_SP_ptgaps[NEtaGaps];
  TProfile2D *p_v3cent_SP_ptgaps[NEtaGaps];

  TProfile2D *p_v2WEST_EP_ptgaps[NEtaGaps];
  TProfile2D *p_v3WEST_EP_ptgaps[NEtaGaps];
  TProfile2D *p_v2EAST_EP_ptgaps[NEtaGaps];
  TProfile2D *p_v3EAST_EP_ptgaps[NEtaGaps];

  TProfile2D *p_vnBBC[fNharmonics + 1];
  // TProfile2D *p_vnZDC[fNharmonics+1];
  //  TProfile2D *p_vnBBCEast[fNharmonics+1];
  //  TProfile2D *p_vnBBCWest[fNharmonics+1];

  TProfile2D *p_v2run_EP[NEtaGaps][16]; // runId vs pt
  TProfile2D *p_v3run_EP[NEtaGaps][16]; // runId vs pt
  TProfile2D *p_v2run_SP[NEtaGaps][16]; // runId vs pt
  TProfile2D *p_v3run_SP[NEtaGaps][16]; // runId vs pt

  TProfile2D *p_v2run_BBC[16]; // runId vs pt
  TProfile2D *p_v2run_ZDC[16]; // runId vs pt
  TProfile2D *p_v2run_ZDC_East[16]; // runId vs pt
  TProfile2D *p_v2run_ZDC_West[16]; // runId vs pt

  //TH2D *p_ptrun[16]; // to calculate <pT> for each pT bin
  //TH2D *p1_ptrun[16]; // to calculate <pT> for each pT bin

  /*TProfile2D *p_v2cent_EP_3subTPC;
  TProfile2D *p_v2cent_EP_3subBBCEast;
  TProfile2D *p_v2cent_EP_3subBBCWest;

  TProfile2D *p_v2_VtxZsys[NEtaGaps];
  TProfile2D *p_v2_DCAsys[NEtaGaps];
  TProfile2D *p_v2_Nhitsys[NEtaGaps];

  TProfile2D *p_v3_VtxZsys[NEtaGaps];
  TProfile2D *p_v3_DCAsys[NEtaGaps];
  TProfile2D *p_v3_Nhitsys[NEtaGaps];*/

  // <vn> vs RunId & centrality - for RunQA
  // In this case, standard EP from 2-sub(TPC) method is used
  //  TProfile2D *p_v2_vs_run;
  //  TProfile2D *p_v3_vs_run;

  //  p_v2cent_EP_3subTPC = new TProfile2D(Form("p_v2cent_EP_3subTPC"),Form("p_v2cent_EP_3subTPC"),16,-0.5,15.5,50,0.15,5.15);
  //  p_v2cent_EP_3subBBCEast = new TProfile2D(Form("p_v2cent_EP_3subBBCEast"),Form("p_v2cent_EP_3subBBCEast"),16,-0.5,15.5,50,0.15,5.15);
  //  p_v2cent_EP_3subBBCWest = new TProfile2D(Form("p_v2cent_EP_3subBBCWest"),Form("p_v2cent_EP_3subBBCWest"),16,-0.5,15.5,50,0.15,5.15);

  //  p_v2_vs_run = new TProfile2D(Form("p_v2_vs_run"),Form("p_v2_vs_run"),fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  //  p_v3_vs_run = new TProfile2D(Form("p_v3_vs_run"),Form("p_v3_vs_run"),fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);

  for (int iGap = 0; iGap < NEtaGaps; iGap++)
  {
    for (int iCent = 0; iCent < 16; iCent++)
    {
      p_v2run_EP[iGap][iCent] = new TProfile2D(Form("p_v2_run_EP_gap%i_cent%i", iGap, iCent), Form("p_v2_run_EP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
      p_v3run_EP[iGap][iCent] = new TProfile2D(Form("p_v3_run_EP_gap%i_cent%i", iGap, iCent), Form("p_v3_run_EP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
      p_v2run_SP[iGap][iCent] = new TProfile2D(Form("p_v2_run_SP_gap%i_cent%i", iGap, iCent), Form("p_v2_run_SP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
      p_v3run_SP[iGap][iCent] = new TProfile2D(Form("p_v3_run_SP_gap%i_cent%i", iGap, iCent), Form("p_v3_run_SP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
      if (iGap == 0)
      {
        p_v2run_BBC[iCent] = new TProfile2D(Form("p_v2_run_BBC_cent%i", iCent), Form("p_v2_run_BBC_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
        p_v2run_ZDC[iCent] = new TProfile2D(Form("p_v2_run_ZDC_cent%i", iCent), Form("p_v2_run_ZDC_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
        p_v2run_ZDC_East[iCent] = new TProfile2D(Form("p_v2_run_ZDC_East_cent%i", iCent), Form("p_v2_run_ZDC_East_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
        p_v2run_ZDC_West[iCent] = new TProfile2D(Form("p_v2_run_ZDC_West_cent%i", iCent), Form("p_v2_run_ZDC_West_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
      }
    }
  }

  for (int iGap = 0; iGap < NEtaGaps; iGap++)
  {
    p_v2cent_EP_ptgaps[iGap] = new TProfile2D(Form("p_v2cent_EP_gap%i_ptgaps", iGap), Form("p_v2cent_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
    p_v3cent_EP_ptgaps[iGap] = new TProfile2D(Form("p_v3cent_EP_gap%i_ptgaps", iGap), Form("p_v3cent_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);

    p_v2cent_SP_ptgaps[iGap] = new TProfile2D(Form("p_v2cent_SP_gap%i_ptgaps", iGap), Form("p_v2cent_SP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
    p_v3cent_SP_ptgaps[iGap] = new TProfile2D(Form("p_v3cent_SP_gap%i_ptgaps", iGap), Form("p_v3cent_SP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);

    p_v2WEST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v2WEST_EP_gap%i_ptgaps", iGap), Form("p_v2WEST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
    p_v2EAST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v2EAST_EP_gap%i_ptgaps", iGap), Form("p_v2EAST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
    p_v3WEST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v3WEST_EP_gap%i_ptgaps", iGap), Form("p_v3WEST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
    p_v3EAST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v3EAST_EP_gap%i_ptgaps", iGap), Form("p_v3EAST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);

    /*p_v2_VtxZsys[iGap] = new TProfile2D(Form("p_v2_VtxZsys%i",iGap),Form("p_v2_VtxZsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
    p_v2_DCAsys[iGap] = new TProfile2D(Form("p_v2_DCAsys%i",iGap),Form("p_v2_DCAsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
    p_v2_Nhitsys[iGap] = new TProfile2D(Form("p_v2_Nhitsys%i",iGap),Form("p_v2_Nhitsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);

    p_v3_VtxZsys[iGap] = new TProfile2D(Form("p_v3_VtxZsys%i",iGap),Form("p_v3_VtxZsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
    p_v3_DCAsys[iGap] = new TProfile2D(Form( "p_v3_DCAsys%i",iGap),Form( "p_v3_DCAsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
    p_v3_Nhitsys[iGap] = new TProfile2D(Form("p_v3_Nhitsys%i",iGap),Form("p_v3_Nhitsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);*/
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    p_vnBBC[iHarm] = new TProfile2D(Form("p_vnBBC%i", iHarm), Form("p_vnBBC%i", iHarm), 16, -0.5, 15.5, 50, 0.15, 5.15);
    // p_vnBBCEast[iHarm] = new TProfile2D(Form("p_vnBBCEast%i",iHarm),Form("p_vnBBCEast%i",iHarm),16,-0.5,15.5,50,0.15,5.15);
    // p_vnBBCWest[iHarm] = new TProfile2D(Form("p_vnBBCWest%i",iHarm),Form("p_vnBBCWest%i",iHarm),16,-0.5,15.5,50,0.15,5.15);
  }

  //Counters
  Int_t fUsedTracks = 0;
  Int_t fQcountFull = 0;
  Int_t fQcountFullEast = 0;
  Int_t fQcountFullWest = 0;
  Int_t fQcountWest[NEtaGaps];
  Int_t fQcountEast[NEtaGaps];

  //Q-vectors
  TVector2 fQv2FullEP;
  TVector2 fQv2FullEPCorr;
  TVector2 fQv2FullSP;
  TVector2 fQv3FullEP;
  TVector2 fQv3FullSP;
  TVector2 fQv2WestEP[NEtaGaps];
  TVector2 fQv2EastEP[NEtaGaps];
  TVector2 fQv2WestSP[NEtaGaps];
  TVector2 fQv2EastSP[NEtaGaps];
  TVector2 fQv3WestEP[NEtaGaps];
  TVector2 fQv3EastEP[NEtaGaps];
  TVector2 fQv3WestSP[NEtaGaps];
  TVector2 fQv3EastSP[NEtaGaps];
  TVector2 QvTMP;
  TVector2 fQv2WestSP_Shift[NEtaGaps];
  TVector2 fQv2EastSP_Shift[NEtaGaps];
  TVector2 fQv3WestSP_Shift[NEtaGaps];
  TVector2 fQv3EastSP_Shift[NEtaGaps];
  TVector2 fUvSP;

  Double_t weight;
  Int_t VtxSign;
  Double_t psiShiftedTMP;
  Double_t PsiSin;
  Double_t PsiCos;

  Long64_t goodEventCounter = 0;

  Double_t Psi2FullEP, Psi3FullEP;
  Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
  Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
  Double_t Psi2EastSP[NEtaGaps], Psi3EastSP[NEtaGaps];
  Double_t Psi2WestSP[NEtaGaps], Psi3WestSP[NEtaGaps];
  Double_t Psi_ReCenter, dPsi;
  Int_t bin;
  Double_t res2EP, res3EP, res2EPsqr, res3EPsqr;
  Double_t flow2EP;
  Double_t flow3EP;
  Double_t res2SP, res3SP, res2SPsqr, res3SPsqr;
  Double_t flow2SP;
  Double_t flow3SP;
  Int_t shiftBin;
  Double_t psiSP_rec;
  Double_t psiEP_rec, psiEP_BBC_Full[fNharmonics + 1], shiftCos, shiftSin;
  Double_t psiEP_BBC_East[fNharmonics + 1], psiEP_BBC_West[fNharmonics + 1];

  // ZDC related EP determination
  Double_t psiEP_ZDC_Full[fNharmonics + 1], psiEP_ZDC_East[fNharmonics + 1], psiEP_ZDC_West[fNharmonics + 1];
  TVector2 qRecEastZDC[fNharmonics + 1], qRecWestZDC[fNharmonics + 1], qRecFullZDC[fNharmonics + 1];
  Double_t qxEast, qyEast, adcxEast, adcyEast;
  Double_t qxWest, qyWest, adcxWest, adcyWest;
  Double_t qxFull, qyFull, adcxFull, adcyFull;
  Double_t qxRecEast, qyRecEast;
  Double_t qxRecWest, qyRecWest;
  Double_t qxRecFull, qyRecFull;
  Double_t resEastWestZDC[fNharmonics + 1];

  // Loop over events
  for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  {

    if ((iEvent + 1) % 1000 == 0)
    {
      std::cout << "Working on event #[" << (iEvent + 1)
                << "/" << events2read << "]" << std::endl;
    }

    Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
    if (!readEvent)
    {
      std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
      break;
    }

    // Retrieve femtoDst
    StFemtoDst *dst = femtoReader->femtoDst();

    // Retrieve event information
    StFemtoEvent *event = dst->event();
    if (!event)
    {
      std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
      break;
    }

    // Event selection
    if (!isGoodEvent(event))
      continue;
    goodEventCounter++;

    TVector3 pVtx = event->primaryVertex();

    // Track analysis
    Int_t nTracks = dst->numberOfTracks();

    //Counters
    fUsedTracks = 0;
    fQcountFull = 0;
    fQcountFullEast = 0;
    fQcountFullWest = 0;
    psiShiftedTMP = 0.;
    PsiSin = 0.;
    PsiCos = 0.;
    fQv2FullEP.Set(0., 0.);
    fQv2FullSP.Set(0., 0.);
    fQv3FullEP.Set(0., 0.);
    fQv3FullSP.Set(0., 0.);
    QvTMP.Set(0., 0.);

    for (int i = 0; i < NEtaGaps; i++)
    {
      fQcountWest[i] = 0;
      fQcountEast[i] = 0;
      fQv2WestEP[i].Set(0., 0.);
      fQv2EastEP[i].Set(0., 0.);
      fQv2WestSP[i].Set(0., 0.);
      fQv2EastSP[i].Set(0., 0.);
      fQv3WestEP[i].Set(0., 0.);
      fQv3EastEP[i].Set(0., 0.);
      fQv3WestSP[i].Set(0., 0.);
      fQv3EastSP[i].Set(0., 0.);
    }

    // Vertex sign
    VtxSign = GetVzBin(pVtx.Z());
    if (VtxSign == -1)
      continue;

    //BBC EP
    Double_t adcNormEastInner = BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 1, 6) / 6.;
    Double_t adcNormEastOuter = BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 7, 16) / 10.;
    Double_t adcNormWestInner = BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 1, 6) / 6.;
    Double_t adcNormWestOuter = BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 7, 16) / 10.;

    float egain[16] = {
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 1)) / adcNormEastInner),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 2)) / adcNormEastInner),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 3)) / adcNormEastInner),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 4)) / adcNormEastInner),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 5)) / adcNormEastInner),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 6)) / adcNormEastInner),
        (float)((0.5) * BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 7)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 8)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 9)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 10)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 11)) / adcNormEastOuter),
        (float)((0.5) * BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 12)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 13)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 14)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 15)) / adcNormEastOuter),
        (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 16)) / adcNormEastOuter)};

    float wgain[16] = {
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 1)) / adcNormWestInner),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 2)) / adcNormWestInner),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 3)) / adcNormWestInner),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 4)) / adcNormWestInner),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 5)) / adcNormWestInner),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 6)) / adcNormWestInner),
        (float)((0.5) * BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 7)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 8)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 9)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 10)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 11)) / adcNormWestOuter),
        (float)((0.5) * BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 12)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 13)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 14)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 15)) / adcNormWestOuter),
        (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 16)) / adcNormWestOuter)};

    float nHitE_Gain[16], nHitW_Gain[16];
    for (int i = 0; i < 16; i++)
    {
      nHitE_Gain[i] = event->bbcAdcEast(i) / egain[i];
      nHitW_Gain[i] = event->bbcAdcWest(i) / wgain[i];
    }

    Double_t psiFull[fNharmonics + 1];
    TVector2 qGainFull[fNharmonics + 1];
    TVector2 qRecFull[fNharmonics + 1];

    Double_t psiEast[fNharmonics + 1];
    TVector2 qGainEast[fNharmonics + 1];
    TVector2 qRecEast[fNharmonics + 1];
    Double_t psiWest[fNharmonics + 1];
    TVector2 qGainWest[fNharmonics + 1];
    TVector2 qRecWest[fNharmonics + 1];

    for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
    {
      //BBC Full
      //v1*y > 0 by conventions most reasonable to the STAR event plane
      //therefore the west event plane has the right sign for the event plane
      //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)

      //The east event plane points opposite that of the west and the west finds the "correct" EP for the STAR coordinates, use West-East
      Double_t QVector_GainX = 0., QVector_GainY = 0.;
      Double_t EventPlaneGainFull = 0.;

      Float_t eXFsum_Gain = 0., wXFsum_Gain = 0., eYFsum_Gain = 0., wYFsum_Gain = 0., eWFgt_Gain = 0., wWFgt_Gain = 0.;
      TVector2 QFullGain;

      for (int iTile = 0; iTile < 16; iTile++)
      {
        eXFsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
        wXFsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];

        eYFsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
        wYFsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];

        eWFgt_Gain += nHitE_Gain[iTile];
        wWFgt_Gain += nHitW_Gain[iTile];
      }

      QVector_GainX = (eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wXFsum_Gain / wWFgt_Gain + eXFsum_Gain / eWFgt_Gain : 0.;
      QVector_GainY = (eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wYFsum_Gain / wWFgt_Gain + eYFsum_Gain / eWFgt_Gain : 0.;
      QFullGain.Set(QVector_GainX, QVector_GainY);
      if (QFullGain.Mod())
      {
        EventPlaneGainFull = QFullGain.Phi();
        if (EventPlaneGainFull < 0.0)
        {
          EventPlaneGainFull += 2. * TMath::Pi() / (iHarm + 1.);
        }
      }
      psiFull[iHarm] = EventPlaneGainFull;
      qGainFull[iHarm] = QFullGain;

      //BBC West
      //v1*y > 0 by conventions most reasonable to the STAR event plane
      //therefore the west event plane has the right sign for the event plane
      //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
      Double_t QVectorW_GainX = 0., QVectorW_GainY = 0.;
      Double_t EventPlaneGainWest = 0.;

      Float_t wXsum_Gain = 0., wYsum_Gain = 0., wWgt_Gain = 0.;
      TVector2 QWestGain;

      for (int iTile = 0; iTile < 16; iTile++)
      {
        wXsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];
        wYsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];
        wWgt_Gain += nHitW_Gain[iTile];
      }
      QVectorW_GainX = (wWgt_Gain > 0.0) ? wXsum_Gain / wWgt_Gain : 0.0;
      QVectorW_GainY = (wWgt_Gain > 0.0) ? wYsum_Gain / wWgt_Gain : 0.0;

      QWestGain.Set(QVectorW_GainX, QVectorW_GainY);
      if (QWestGain.Mod())
      {
        EventPlaneGainWest = QWestGain.Phi();
        if (EventPlaneGainWest < 0.0)
        {
          EventPlaneGainWest += 2. * TMath::Pi() / (iHarm + 1.);
        }
      }
      psiWest[iHarm] = EventPlaneGainWest;
      qGainWest[iHarm] = QWestGain;

      //BBC East
      //in BBC_GetPhi the fact that the tile numbering scheme is reversed about the y axis for the east BBC is accounted for
      //The east event plane points opposite that of the west since the spectators are on different sides -- see west comments
      Double_t QVectorE_GainX = 0., QVectorE_GainY = 0.;
      Double_t EventPlaneGainEast = 0.;

      Float_t eXsum_Gain = 0., eYsum_Gain = 0., eWgt_Gain = 0.;
      TVector2 QEastGain;

      for (int iTile = 0; iTile < 16; iTile++)
      {
        eXsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
        eYsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
        eWgt_Gain += nHitE_Gain[iTile];
      }
      QVectorE_GainX = (eWgt_Gain > 0.0) ? eXsum_Gain / eWgt_Gain : 0.0;
      QVectorE_GainY = (eWgt_Gain > 0.0) ? eYsum_Gain / eWgt_Gain : 0.0;

      QEastGain.Set(QVectorE_GainX, QVectorE_GainY);
      if (QEastGain.Mod())
      {
        EventPlaneGainEast = QEastGain.Phi();
        if (EventPlaneGainEast < 0.0)
        {
          EventPlaneGainEast += 2. * TMath::Pi() / (iHarm + 1.);
        }
      }
      psiEast[iHarm] = EventPlaneGainEast;
      qGainEast[iHarm] = QEastGain;
    }

    // Apply recentering
    for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
    {
      TVector2 recFullTMP((Double_t)p_BBC_Qx_Full_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
                          (Double_t)p_BBC_Qy_Full_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      TVector2 recEastTMP((Double_t)p_BBC_Qx_East_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
                          (Double_t)p_BBC_Qy_East_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      TVector2 recWestTMP((Double_t)p_BBC_Qx_West_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
                          (Double_t)p_BBC_Qy_West_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      qRecEast[iHarm] = qGainEast[iHarm] - recEastTMP;
      qRecWest[iHarm] = qGainWest[iHarm] - recWestTMP;
      qRecFull[iHarm] = qGainFull[iHarm] - recFullTMP;
    }

    for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
    {
      // Full BBC
      if (qRecFull[iHarm].Mod())
      {
        psiEP_rec = qRecFull[iHarm].Phi() / (iHarm + 1.);
        if (psiEP_rec < 0.)
        {
          psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
        }
      }
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
      {
        shiftCos = p_BBC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        shiftSin = p_BBC_Sin_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
                (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
                 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
      }
      psiEP_BBC_Full[iHarm] = psiEP_rec + dPsi;

      if (psiEP_BBC_Full[iHarm] < 0.0)
        psiEP_BBC_Full[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
      if (psiEP_BBC_Full[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
        psiEP_BBC_Full[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);

      // East BBC
      if (qRecEast[iHarm].Mod())
      {
        psiEP_rec = qRecEast[iHarm].Phi() / (iHarm + 1.);
        if (psiEP_rec < 0.)
        {
          psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
        }
      }
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
      {
        shiftCos = p_BBC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        shiftSin = p_BBC_Sin_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
                (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
                 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
      }
      psiEP_BBC_East[iHarm] = psiEP_rec + dPsi;

      if (psiEP_BBC_East[iHarm] < 0.0)
        psiEP_BBC_East[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
      if (psiEP_BBC_East[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
        psiEP_BBC_East[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);

      // West BBC
      if (qRecWest[iHarm].Mod())
      {
        psiEP_rec = qRecWest[iHarm].Phi() / (iHarm + 1.);
        if (psiEP_rec < 0.)
        {
          psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
        }
      }
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
      {
        shiftCos = p_BBC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        shiftSin = p_BBC_Sin_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
                (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
                 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
      }
      psiEP_BBC_West[iHarm] = psiEP_rec + dPsi;

      if (psiEP_BBC_West[iHarm] < 0.0)
        psiEP_BBC_West[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
      if (psiEP_BBC_West[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
        psiEP_BBC_West[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
    }

    //EP from BBC is ready: psiEP_BBC_Full[iHarm], where iHarm = 0 for 1st harmonics, 1 for 2nd & 2 for 3rd.

    // ZDC EP
    if (event->zdcSumAdcEast() < 1.)
      continue;
    if (event->zdcSumAdcWest() < 1.)
      continue;

    qxEast = 0.;
    qyEast = 0.;
    adcxEast = 0.;
    adcyEast = 0.;
    qxWest = 0.;
    qyWest = 0.;
    adcxWest = 0.;
    adcyWest = 0.;
    qxFull = 0.;
    qyFull = 0.;
    adcxFull = 0.;
    adcyFull = 0.;

    // Vertical ZDC-SMD
    for (int iStrip = 0; iStrip < fNZdcSmdStripsVertical; iStrip++)
    {
      qxEast += event->zdcSmdEastVertical(iStrip) * GetZDCPosition(0, 0, iStrip);
      adcxEast += event->zdcSmdEastVertical(iStrip);

      qxWest += event->zdcSmdWestVertical(iStrip) * GetZDCPosition(1, 0, iStrip);
      adcxWest += event->zdcSmdWestVertical(iStrip);
    }
    // Horizontal ZDC-SMD
    for (int iStrip = 0; iStrip < fNZdcSmdStripsHorizontal; iStrip++)
    {
      qyEast += event->zdcSmdEastHorizontal(iStrip) * GetZDCPosition(0, 1, iStrip);
      adcyEast += event->zdcSmdEastHorizontal(iStrip);

      qyWest += event->zdcSmdWestHorizontal(iStrip) * GetZDCPosition(1, 1, iStrip);
      adcyWest += event->zdcSmdWestHorizontal(iStrip);
    }

    if (adcxEast > 0)
      qxEast /= adcxEast;
    else
      qxEast = -999.;
    if (adcyEast > 0)
      qyEast /= adcyEast;
    else
      qyEast = -999.;
    if (adcxWest > 0)
      qxWest /= adcxWest;
    else
      qxWest = -999.;
    if (adcyWest > 0)
      qyWest /= adcyWest;
    else
      qyWest = -999.;

    // Full ZDC-SMD Q-vector
    if (qxEast == -999. || qxWest == -999.)
      qxFull = -999.;
    else
      qxFull = qxEast - qxWest;
    if (qyEast == -999. || qyWest == -999.)
      qyFull = -999.;
    else
      qyFull = qyEast - qyWest;

    if (qxEast != -999.)
      qxRecEast = qxEast - (Double_t)p_ZDC_Qx_East_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
    if (qyEast != -999.)
      qyRecEast = qyEast - (Double_t)p_ZDC_Qy_East_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
    if (qxWest != -999.)
      qxRecWest = qxWest - (Double_t)p_ZDC_Qx_West_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
    if (qyWest != -999.)
      qyRecWest = qyWest - (Double_t)p_ZDC_Qy_West_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
    if (qxFull != -999.)
      qxRecFull = qxFull - (Double_t)p_ZDC_Qx_Full_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
    if (qyFull != -999.)
      qyRecFull = qyFull - (Double_t)p_ZDC_Qy_Full_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];

    if (qxEast == -999.)
      qxRecEast = qxEast;
    if (qyEast == -999.)
      qyRecEast = qyEast;
    if (qxWest == -999.)
      qxRecWest = qxWest;
    if (qyWest == -999.)
      qyRecWest = qyWest;
    if (qxFull == -999.)
      qxRecFull = qxFull;
    if (qyFull == -999.)
      qyRecFull = qyFull;

    PsiSin = 0.;
    PsiCos = 0.;
    psiShiftedTMP = 0.;

    for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
    {
      qRecEastZDC[iHarm].Set(qxRecEast, qyRecEast);
      qRecWestZDC[iHarm].Set(qxRecWest, qyRecWest);
      qRecFullZDC[iHarm].Set(qxRecFull, qyRecFull);

      // Full ZDC
      if (qRecFullZDC[iHarm].Mod())
      {
        psiEP_rec = qRecFullZDC[iHarm].Phi() / (iHarm + 1.);
        if (psiEP_rec < 0.)
        {
          psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
        }
      }
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
      {
        shiftCos = p_ZDC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        shiftSin = p_ZDC_Sin_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
                (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
                 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
      }
      psiEP_ZDC_Full[iHarm] = psiEP_rec + dPsi;

      if (psiEP_ZDC_Full[iHarm] < 0.0)
        psiEP_ZDC_Full[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
      if (psiEP_ZDC_Full[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
        psiEP_ZDC_Full[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);

      // East ZDC
      if (qRecEastZDC[iHarm].Mod())
      {
        psiEP_rec = qRecEastZDC[iHarm].Phi() / (iHarm + 1.);
        if (psiEP_rec < 0.)
        {
          psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
        }
      }
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
      {
        shiftCos = p_ZDC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        shiftSin = p_ZDC_Sin_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
                (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
                 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
      }
      psiEP_ZDC_East[iHarm] = psiEP_rec + dPsi;

      if (psiEP_ZDC_East[iHarm] < 0.0)
        psiEP_ZDC_East[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
      if (psiEP_ZDC_East[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
        psiEP_ZDC_East[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);

      // West ZDC
      if (qRecWestZDC[iHarm].Mod())
      {
        psiEP_rec = qRecWestZDC[iHarm].Phi() / (iHarm + 1.);
        if (psiEP_rec < 0.)
        {
          psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
        }
      }
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
      {
        shiftCos = p_ZDC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        shiftSin = p_ZDC_Sin_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
                (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
                 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
      }
      psiEP_ZDC_West[iHarm] = psiEP_rec + dPsi;

      if (psiEP_ZDC_West[iHarm] < 0.0)
        psiEP_ZDC_West[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
      if (psiEP_ZDC_West[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
        psiEP_ZDC_West[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
    }

    // Track loop
    for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
    {

      // Retrieve i-th femto track
      StFemtoTrack *femtoTrack = dst->track(iTrk);

      if (!femtoTrack)
        continue;

      // Must be a primary track
      if (!femtoTrack->isPrimary())
        continue;

      //Track selection
      if (!isGoodTrack(femtoTrack, energy, pVtx))
        continue;

      fUsedTracks++;

      //Fill recentered Q-vectors for all TPC (FULL)
      weight = GetWeight(femtoTrack);

      QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
      QvTMP.SetY(p_Qx_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      fQv2FullEP += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
      QvTMP.Set(0., 0.);

      QvTMP.SetX(p_Qx_FULL_EP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
      QvTMP.SetY(p_Qy_FULL_EP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      fQv3FullEP += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
      QvTMP.Set(0., 0.);

      QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
      QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      fQv2FullSP += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
      QvTMP.Set(0., 0.);

      QvTMP.SetX(p_Qx_FULL_SP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
      QvTMP.SetY(p_Qy_FULL_SP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

      fQv3FullSP += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
      QvTMP.Set(0., 0.);

      fQcountFull++;
      if (femtoTrack->eta() >= 0.)
      {
        fQcountFullWest++;
      }
      else
      {
        fQcountFullEast++;
      }

      for (int iGap = 0; iGap < NEtaGaps; iGap++)
      {
        //fill Qvectors from EAST arm
        if (femtoTrack->eta() > -1. * cutEta && femtoTrack->eta() < 0.) //cutEtaGap[iGap])
        {
          QvTMP.SetX(p_Qx_EAST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_EAST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv2EastEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
          QvTMP.Set(0., 0.);

          QvTMP.SetX(p_Qx_EAST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_EAST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv3EastEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
          QvTMP.Set(0., 0.);

          QvTMP.SetX(p_Qx_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv2EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
          QvTMP.Set(0., 0.);

          QvTMP.SetX(p_Qx_EAST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_EAST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv3EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
          QvTMP.Set(0., 0.);

          fQcountEast[iGap]++;
        }
        //fill Qvectors from WEST arm
        if (femtoTrack->eta() < 1. * cutEta && femtoTrack->eta() > 0.) //cutEtaGap[iGap])
        {
          QvTMP.SetX(p_Qx_WEST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_WEST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv2WestEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
          QvTMP.Set(0., 0.);

          QvTMP.SetX(p_Qx_WEST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_WEST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv3WestEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
          QvTMP.Set(0., 0.);

          QvTMP.SetX(p_Qx_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv2WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
          QvTMP.Set(0., 0.);

          QvTMP.SetX(p_Qx_WEST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_WEST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);

          fQv3WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
          QvTMP.Set(0., 0.);

          fQcountWest[iGap]++;
        }
      }

      // Check if track has TOF signal
      if (femtoTrack->isTofTrack())
      {

      } //if( isTofTrack() )

    } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)

    //Getting PsiEP with recentering + shift correction
    Psi_ReCenter = 0.;
    Double_t mCos, mSin;
    if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
    {
      Psi2FullEP = 0.;
      Psi3FullEP = 0.;

      //FULL 2-nd order
      Psi_ReCenter = TMath::ATan2(fQv2FullEP.Y(), fQv2FullEP.X()) / 2.;
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax; iK++)
      {
        mCos = p_Cos2_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        mSin = p_Cos2_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];

        dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
      }
      Psi2FullEP = AngleShift(Psi_ReCenter + dPsi, 2.);
    }
    for (int iGap = 0; iGap < NEtaGaps; iGap++)
    {
      if (fQcountEast[iGap] > cutNcountQvGap &&
          fQcountWest[iGap] > cutNcountQvGap)
      {
        Psi2EastEP[iGap] = 0.;
        Psi3EastEP[iGap] = 0.;
        Psi2WestEP[iGap] = 0.;
        Psi3WestEP[iGap] = 0.;
        Psi2EastSP[iGap] = 0.;
        Psi3EastSP[iGap] = 0.;
        Psi2WestSP[iGap] = 0.;
        Psi3WestSP[iGap] = 0.;

        //EAST 2-nd order
        Psi_ReCenter = TMath::ATan2(fQv2EastEP[iGap].Y(), fQv2EastEP[iGap].X()) / 2.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos2_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos2_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
        }
        Psi2EastEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);

        Psi_ReCenter = TMath::ATan2(fQv2EastSP[iGap].Y(), fQv2EastSP[iGap].X()) / 2.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos2_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos2_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
        }
        Psi2EastSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
        fQv2EastSP_Shift[iGap].Set(fQv2EastSP[iGap].Mod() * TMath::Cos(2. * Psi2EastSP[iGap]), fQv2EastSP[iGap].Mod() * TMath::Sin(2. * Psi2EastSP[iGap]));

        //WEST 2-nd order
        Psi_ReCenter = TMath::ATan2(fQv2WestEP[iGap].Y(), fQv2WestEP[iGap].X()) / 2.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos2_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos2_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
        }
        Psi2WestEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);

        Psi_ReCenter = TMath::ATan2(fQv2WestSP[iGap].Y(), fQv2WestSP[iGap].X()) / 2.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos2_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos2_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
        }
        Psi2WestSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
        fQv2WestSP_Shift[iGap].Set(fQv2WestSP[iGap].Mod() * TMath::Cos(2. * Psi2WestSP[iGap]), fQv2WestSP[iGap].Mod() * TMath::Sin(2. * Psi2WestSP[iGap]));

        //EAST 3-nd order
        Psi_ReCenter = TMath::ATan2(fQv3EastEP[iGap].Y(), fQv3EastEP[iGap].X()) / 3.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos3_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos3_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
        }
        Psi3EastEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);

        Psi_ReCenter = TMath::ATan2(fQv3EastSP[iGap].Y(), fQv3EastSP[iGap].X()) / 3.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos3_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos3_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
        }
        Psi3EastSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
        fQv3EastSP_Shift[iGap].Set(fQv3EastSP[iGap].Mod() * TMath::Cos(3. * Psi3EastSP[iGap]), fQv3EastSP[iGap].Mod() * TMath::Sin(3. * Psi3EastSP[iGap]));

        //WEST 3-nd order
        Psi_ReCenter = TMath::ATan2(fQv3WestEP[iGap].Y(), fQv3WestEP[iGap].X()) / 3.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos3_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos3_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
        }
        Psi3WestEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);

        Psi_ReCenter = TMath::ATan2(fQv3WestSP[iGap].Y(), fQv3WestSP[iGap].X()) / 3.;
        dPsi = 0.;
        for (int iK = 0; iK < NShiftOrderMax; iK++)
        {
          mCos = p_Cos3_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
          mSin = p_Cos3_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];

          dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
        }
        Psi3WestSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
        fQv3WestSP_Shift[iGap].Set(fQv3WestSP[iGap].Mod() * TMath::Cos(3. * Psi3WestSP[iGap]), fQv3WestSP[iGap].Mod() * TMath::Sin(3. * Psi3WestSP[iGap]));
      }
    }

    //FLOW CALCULATIONS
    // Second Track loop
    for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
    {

      // Retrieve i-th femto track
      StFemtoTrack *femtoTrack = dst->track(iTrk);

      if (!femtoTrack)
        continue;

      // Must be a primary track
      if (!femtoTrack->isPrimary())
        continue;

      //Track selection
      if (isGoodTrackFlow(femtoTrack, energy, pVtx))
      {

        weight = GetWeight(femtoTrack);

        // Flow with BBC EP
        for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
        {
          res2EP = hResBBCFULL[iHarm]->GetBinContent(hResBBCFULL[iHarm]->FindBin(event->cent16()));
          flow2EP = TMath::Cos((iHarm + 1.) * (femtoTrack->phi() - psiEP_BBC_Full[iHarm])) / res2EP;
          p_vnBBC[iHarm]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);

          if (iHarm == 1)
          {
            p_v2run_BBC[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
          }
          /*res2EP = hResBBCHalf[iHarm]->GetBinContent(hResBBCHalf[iHarm]->FindBin(event->cent16()));
          flow2EP = TMath::Cos((iHarm+1.)*(femtoTrack->phi()-psiEP_BBC_East[iHarm])) / res2EP;
          p_vnBBCEast[iHarm]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);

          res2EP = hResBBCHalf[iHarm]->GetBinContent(hResBBCHalf[iHarm]->FindBin(event->cent16()));
          flow2EP = TMath::Cos((iHarm+1.)*(femtoTrack->phi()-psiEP_BBC_West[iHarm])) / res2EP;
          p_vnBBCWest[iHarm]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);*/
        }

        // Flow with ZDC EP
        // Full
        res2EP = hResZDCFULL[1]->GetBinContent(hResZDCFULL[1]->FindBin(event->cent16()));
        flow2EP = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_Full[0])) / res2EP;
        p_v2run_ZDC[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
        // East
        res2EP = hResZDCHalf[1]->GetBinContent(hResZDCHalf[1]->FindBin(event->cent16()));
        flow2EP = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_East[0])) / res2EP;
        p_v2run_ZDC_East[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
        // West
        res2EP = hResZDCHalf[1]->GetBinContent(hResZDCHalf[1]->FindBin(event->cent16()));
        flow2EP = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_West[0])) / res2EP;
        p_v2run_ZDC_West[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
        ///////////////////////////////////////////////////////////////////////////////////////////
        //p_ptrun[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), weight);
        //p1_ptrun[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), 1.);
        /*
        // 3-sub method
      
        //Get rid of autocorrelations
        //Since pt<2 GeV/c cut was used in EP determination, only tracks within this range should be corrected
        if (isGoodTrack(femtoTrack, energy, pVtx))
        {
          QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign]->GetBinContent(p_Qx_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent16())));
          QvTMP.SetY(p_Qy_FULL_EP[0][VtxSign]->GetBinContent(p_Qy_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent16())));

          fQv2FullEPCorr = fQv2FullEP - weight*(CalcQvector(femtoTrack,2) - QvTMP);
          QvTMP.Set(0.,0.);
          Psi_ReCenter = TMath::ATan2(fQv2FullEPCorr.Y(), fQv2FullEPCorr.X()) / 2.;
          dPsi = 0.;
          for (int iK = 0; iK < NShiftOrderMax; iK++)
          {
            bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t)event->runId(), (Double_t)event->cent16());
            mCos = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
            bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t)event->runId(), (Double_t)event->cent16());
            mSin = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);

            dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
          }
          Psi2FullEP = AngleShift(Psi_ReCenter + dPsi, 2.);
        }

        res2EPsqr = p_res2_EP_3sub_AB->GetBinContent(p_res2_EP_3sub_AB->FindBin(event->cent16())) *
                    p_res2_EP_3sub_AC->GetBinContent(p_res2_EP_3sub_AC->FindBin(event->cent16())) /
                    p_res2_EP_3sub_BC->GetBinContent(p_res2_EP_3sub_BC->FindBin(event->cent16()));
        res2EP = TMath::Sqrt(res2EPsqr);
        flow2EP = TMath::Cos(2.*(femtoTrack->phi()-Psi2FullEP)) / res2EP;
        p_v2cent_EP_3subTPC->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);

        // Calculate flow
        res2EPsqr = p_res2_EP_3sub_AB->GetBinContent(p_res2_EP_3sub_AB->FindBin(event->cent16())) *
                    p_res2_EP_3sub_BC->GetBinContent(p_res2_EP_3sub_BC->FindBin(event->cent16())) /
                    p_res2_EP_3sub_AC->GetBinContent(p_res2_EP_3sub_AC->FindBin(event->cent16()));
        res2EP = TMath::Sqrt(res2EPsqr);
        flow2EP = TMath::Cos(2.*(femtoTrack->phi()-psiEP_BBC_East[1])) / res2EP;
        p_v2cent_EP_3subBBCEast->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);

        res2EPsqr = p_res2_EP_3sub_AC->GetBinContent(p_res2_EP_3sub_AC->FindBin(event->cent16())) *
                    p_res2_EP_3sub_BC->GetBinContent(p_res2_EP_3sub_BC->FindBin(event->cent16())) /
                    p_res2_EP_3sub_AB->GetBinContent(p_res2_EP_3sub_AB->FindBin(event->cent16()));
        res2EP = TMath::Sqrt(res2EPsqr);
        flow2EP = TMath::Cos(2.*(femtoTrack->phi()-psiEP_BBC_West[1])) / res2EP;
        p_v2cent_EP_3subBBCWest->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
        ///////////////////////////////////////////////////////////////////////////////////////////
*/
        for (int iGap = 0; iGap < NEtaGaps; iGap++)
        {
          res2EPsqr = p_res2_EP[iGap]->GetBinContent(p_res2_EP[iGap]->FindBin(event->cent16()));
          res2EP = TMath::Sqrt(res2EPsqr);
          res3EPsqr = p_res3_EP[iGap]->GetBinContent(p_res3_EP[iGap]->FindBin(event->cent16()));
          res3EP = TMath::Sqrt(res3EPsqr);

          res2SPsqr = p_res2_SP[iGap]->GetBinContent(p_res2_SP[iGap]->FindBin(event->cent16()));
          if (res2SPsqr < 0)
            std::cout << "Event " << iEvent << " | Res2SP NEGATIVEE: " << res2SPsqr << std::endl;
          if (res2SPsqr == 0)
            std::cout << "Res2SP ZEROOOOO" << std::endl;
          res2SP = TMath::Sqrt(res2SPsqr);
          res3SPsqr = p_res3_SP[iGap]->GetBinContent(p_res3_SP[iGap]->FindBin(event->cent16()));
          res3SP = TMath::Sqrt(res3SPsqr);
          //EAST
          //if (track->eta() > -1.*cutEta && track->eta() < 0.)
          if (femtoTrack->eta() > -1. * cutEta && femtoTrack->eta() < -1. * cutEtaGap.at(iGap))
          {
            flow2EP = TMath::Cos(2. * (femtoTrack->phi() - Psi2WestEP[iGap])) / res2EP;
            flow3EP = TMath::Cos(3. * (femtoTrack->phi() - Psi3WestEP[iGap])) / res3EP;
            p_v2cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
            p_v3cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);
            p_v2EAST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
            p_v3EAST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);

            p_v2run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
            p_v3run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP, weight);
            /*if (pVtx.Z() < 20.)
	    {
	      p_v2_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
	      p_v3_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
	    }
            if (femtoTrack->gDCA(pVtx).Mag() < 2.)
	    {
	      p_v2_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
	      p_v3_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
	    }
            if (femtoTrack->nHitsFit() > 10.)
	    {
	      p_v2_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
	      p_v3_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
	    }

            if (iGap == 0)
            {
              p_v2_vs_run->Fill(event->runId(),event->cent16(),flow2EP);
              p_v3_vs_run->Fill(event->runId(),event->cent16(),flow3EP);
            }*/

            fUvSP.Set(1. * TMath::Cos(2. * femtoTrack->phi()), 1. * TMath::Sin(2. * femtoTrack->phi()));
            flow2SP = fUvSP * fQv2WestSP_Shift[iGap] / res2SP;
            fUvSP.Set(1. * TMath::Cos(3. * femtoTrack->phi()), 1. * TMath::Sin(3. * femtoTrack->phi()));
            flow3SP = fUvSP * fQv3WestSP_Shift[iGap] / res3SP;

            p_v2cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2SP, weight);
            p_v3cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3SP, weight);

            p_v2run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2SP, weight);
            p_v3run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3SP, weight);
          }
          //WEST
          //if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 0.)
          if (femtoTrack->eta() < 1. * cutEta && femtoTrack->eta() > 1. * cutEtaGap.at(iGap))
          {
            flow2EP = TMath::Cos(2. * (femtoTrack->phi() - Psi2EastEP[iGap])) / res2EP;
            flow3EP = TMath::Cos(3. * (femtoTrack->phi() - Psi3EastEP[iGap])) / res3EP;
            p_v2cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
            p_v3cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);
            p_v2WEST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
            p_v3WEST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);

            p_v2run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
            p_v3run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP, weight);
            /*if (pVtx.Z() < 20.)
	    {
	      p_v2_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
	      p_v3_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
	    }
            if (femtoTrack->gDCA(pVtx).Mag() < 2.)
	    {
	      p_v2_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
	      p_v3_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
	    }
            if (femtoTrack->nHitsFit() > 10.)
	    {
	      p_v2_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
	      p_v3_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
	    }

            if (iGap == 0)
            {
              p_v2_vs_run->Fill(event->runId(),event->cent16(),flow2EP);
              p_v3_vs_run->Fill(event->runId(),event->cent16(),flow3EP);
            }*/

            fUvSP.Set(1. * TMath::Cos(2. * femtoTrack->phi()), 1. * TMath::Sin(2. * femtoTrack->phi()));
            flow2SP = fUvSP * fQv2EastSP_Shift[iGap] / res2SP;
            fUvSP.Set(1. * TMath::Cos(3. * femtoTrack->phi()), 1. * TMath::Sin(3. * femtoTrack->phi()));
            flow3SP = fUvSP * fQv3EastSP_Shift[iGap] / res3SP;

            p_v2cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2SP, weight);
            p_v3cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3SP, weight);

            p_v2run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2SP, weight);
            p_v3run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3SP, weight);
          }
        }
      }
    }

  } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)

  femtoReader->Finish();

  TFile *output = new TFile(outFile, "recreate");

  for (int iGap = 0; iGap < NEtaGaps; iGap++)
  {
    // p_v2cent_EP_ptgaps[iGap]->Write();
    // p_v3cent_EP_ptgaps[iGap]->Write();
    // p_v2cent_SP_ptgaps[iGap]->Write();
    // p_v3cent_SP_ptgaps[iGap]->Write();
    p_v2EAST_EP_ptgaps[iGap]->Write();
    p_v3EAST_EP_ptgaps[iGap]->Write();
    p_v2WEST_EP_ptgaps[iGap]->Write();
    p_v3WEST_EP_ptgaps[iGap]->Write();
    for (int iCent = 0; iCent < 16; iCent++)
    {
      p_v2run_EP[iGap][iCent]->Write();
      p_v3run_EP[iGap][iCent]->Write();
      p_v2run_SP[iGap][iCent]->Write();
      p_v3run_SP[iGap][iCent]->Write();
      if (iGap == 0)
      {
        p_v2run_BBC[iCent]->Write();
      }
    }
    for (int iCent = 0; iCent < 16; iCent++)
    {
      if (iGap == 0)
      {
        p_v2run_ZDC[iCent]->Write();
        p_v2run_ZDC_East[iCent]->Write();
        p_v2run_ZDC_West[iCent]->Write();
      }
    }
    for (int iCent = 0; iCent < 16; iCent++)
    {
      if (iGap == 0)
      {
        //p_ptrun[iCent]->Write();
        //p1_ptrun[iCent]->Write();
      }
    }
    /*p_v2_VtxZsys[iGap]->Write();
    p_v2_DCAsys[iGap]->Write();
    p_v2_Nhitsys[iGap]->Write();
    p_v3_VtxZsys[iGap]->Write();
    p_v3_DCAsys[iGap]->Write();
    p_v3_Nhitsys[iGap]->Write();*/
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    p_vnBBC[iHarm]->Write();
    //p_vnBBCEast[iHarm]->Write();
    //p_vnBBCWest[iHarm]->Write();
  }
  /*p_v2cent_EP_3subTPC->Write();
  p_v2cent_EP_3subBBCEast->Write();
  p_v2cent_EP_3subBBCWest->Write();

  p_v2_vs_run->Write();
  p_v3_vs_run->Write();*/

  output->Close();

  std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
            << std::endl;
  std::cout << "Good events: " << goodEventCounter << std::endl;
  timer.Stop();
  timer.Print();
}

Bool_t isGoodEvent(StFemtoEvent *const &event)
{
  if (!event)
    return false;
  if (event == nullptr)
    return false;

  if (event->primaryVertex().Perp() > cutVtxR)
    return false;
  if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy))
    return false;

  if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz)
    return false;

  return true;
}

Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
{
  if (!track)
    return false;
  // if (!track->isPrimary()) return false;
  if (track->nHitsFit() < cutNhits)
    return false;
  if (track->nHitsPoss() <= cutNhitsPoss)
    return false;
  if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
    return false;
  if (TMath::Abs(track->eta()) >= cutEta)
    return false;

  if (track->pt() <= cutPtMin.at(_energy))
    return false;
  if (track->pt() > cutPtMax)
    return false;
  if (track->ptot() > cutPMax)
    return false;
  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
    return false;
  return true;
}

Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
{
  if (!track)
    return false;
  // if (!track->isPrimary()) return false;
  if (track->nHitsFit() < cutNhits)
    return false;
  if (track->nHitsPoss() <= cutNhitsPoss)
    return false;
  if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
    return false;
  if (TMath::Abs(track->eta()) >= cutEta)
    return false;

  if (track->pt() <= cutPtMin.at(_energy))
    return false;
  //if (track->pt() > cutPtMax) return false;
  if (track->ptot() > cutPMax)
    return false;
  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
    return false;
  return true;
}

Double_t GetWeight(StFemtoTrack *const &track)
{
  Double_t weight;
  if (track->pt() <= cutPtWeightEP)
  {
    weight = track->pt();
  }
  else
  {
    weight = cutPtWeightEP;
  }
  return weight;
}

TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
{
  TVector2 qv(0., 0.);

  qv.Set(TMath::Cos(_harm * track->phi()), TMath::Sin(_harm * track->phi()));

  return qv;
}

Double_t AngleShift(Double_t Psi, Double_t order)
{
  Double_t PsiCorr = Psi;
  if (PsiCorr > TMath::Pi() / order)
  {
    PsiCorr = Psi - 2. * TMath::Pi() / order;
  }
  if (PsiCorr < -1. * TMath::Pi() / order)
  {
    PsiCorr = Psi + 2. * TMath::Pi() / order;
  }
  return PsiCorr;
}

Int_t GetVzBin(Double_t vtxZ)
{
  for (Int_t i = 0; i < fArraySize; i++)
  {
    if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i + 1])
      return i;
  }
  return -1;
}

//--------------------------------------------------------------------------------------------------------------------------//
//this function simply connects the gain values read in to the BBC azimuthal distribution
//since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
//I prefer assigning the angle between the tiles thus "greying out" the adcs.
//Others have assigned all of the adc to one (exclusive) or the the other.
Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId)
{
  //float GetPhiInBBC(int eastWest, int bbcN) { //tileId=0 to 23
  const float Pi = TMath::Pi();
  const float Phi_div = Pi / 6;
  float bbc_phi = Phi_div;
  switch (tileId)
  {
  case 0:
    bbc_phi = 3. * Phi_div;
    break;
  case 1:
    bbc_phi = Phi_div;
    break;
  case 2:
    bbc_phi = -1. * Phi_div;
    break;
  case 3:
    bbc_phi = -3. * Phi_div;
    break;
  case 4:
    bbc_phi = -5. * Phi_div;
    break;
  case 5:
    bbc_phi = 5. * Phi_div;
    break;
  //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div;	//tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
  case 6:
    bbc_phi = 3. * Phi_div;
    break;
  case 7:
    bbc_phi = 3. * Phi_div;
    break;
  case 8:
    bbc_phi = Phi_div;
    break;
  case 9:
    bbc_phi = 0.;
    break;
  case 10:
    bbc_phi = -1. * Phi_div;
    break;
  //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div;	//tiles 13 and 15 are gained together
  case 11:
    bbc_phi = -3. * Phi_div;
    break;
  case 12:
    bbc_phi = -3. * Phi_div;
    break;
  case 13:
    bbc_phi = -5. * Phi_div;
    break;
  case 14:
    bbc_phi = Pi;
    break;
  case 15:
    bbc_phi = 5. * Phi_div;
    break;
  }

  //if we're looking at the east BBC we need to flip around x in the STAR coordinates,
  //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the
  if (0 == eastWest)
  {
    if (bbc_phi > -0.001)
    { //this is not a >= since we are talking about finite adcs -- not to important
      bbc_phi = Pi - bbc_phi;
    }
    else
    {
      bbc_phi = -Pi - bbc_phi;
    }
  }

  if (bbc_phi < 0.0)
    bbc_phi += 2. * Pi;
  if (bbc_phi > 2. * Pi)
    bbc_phi -= 2. * Pi;

  return bbc_phi;
}

Bool_t isBadRun(StFemtoEvent *const &event)
{
  Int_t run = (Int_t)event->runId();

  if (run == 12132026 ||
      run == 12134014 ||
      run == 12134051 ||
      run == 12135008 ||
      run == 12136006 ||
      run == 12136025 ||
      run == 12136047 ||
      run == 12136064)
    return true;
  else
    return false;
}

Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip)
// Get position of each slat;strip starts from 0
{

  Double_t zdcsmd_x[7] = {0.5, 2, 3.5, 5, 6.5, 8, 9.5};
  Double_t zdcsmd_y[8] = {1.25, 3.25, 5.25, 7.25, 9.25, 11.25, 13.25, 15.25};

  Double_t mZDCSMDCenterex = 4.72466;
  Double_t mZDCSMDCenterey = 5.53629;
  Double_t mZDCSMDCenterwx = 4.39604;
  Double_t mZDCSMDCenterwy = 5.19968;

  if (eastwest == 0 && verthori == 0)
    return zdcsmd_x[strip] - mZDCSMDCenterex;
  if (eastwest == 1 && verthori == 0)
    return mZDCSMDCenterwx - zdcsmd_x[strip];
  if (eastwest == 0 && verthori == 1)
    return zdcsmd_y[strip] / sqrt(2.) - mZDCSMDCenterey;
  if (eastwest == 1 && verthori == 1)
    return zdcsmd_y[strip] / sqrt(2.) - mZDCSMDCenterwy;

  return -999.;
}