StFlowMephi/FemtoDstAnalyzer_FlowPIDsys1.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 "TF1.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 isGoodTrackFlowPID(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);

Bool_t isGoodTrackFlowPIDsys_Nhits(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount);

Bool_t isGoodTrackFlowPIDsys_DCA(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount);

Bool_t isGoodTrackFlowPIDsys_Eta(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount);

Bool_t isGoodPID(StFemtoTrack *const &track);

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);

//Int_t DoPID(StFemtoTrack *const &track);
Double_t GetSigmWidth(StFemtoTrack *const &track, Int_t PID);

Double_t GetMeanNsig(StFemtoTrack *const &track, const Int_t PID);

Double_t GetWidthM2(StFemtoTrack *const &track, const Int_t PID);

Double_t GetRotAngle(StFemtoTrack *const &track);

std::pair<Double_t, Double_t> GetNewXY(StFemtoTrack *const &track);

Bool_t isMesonXY(Double_t x, Double_t y);

Double_t GetSigmNewX(StFemtoTrack *const &track, Int_t PID);

Double_t GetSigmNewY(StFemtoTrack *const &track, Int_t PID);

Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);

Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip);

//_________________
void FemtoDstAnalyzer_FlowPIDsys1(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
                                     const Char_t *outFile = "./oFlowPIDTest.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 = "",
                                     const Char_t *pidFile = "")
{

  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

//  TProfile2D *p_v2run_EP[3][NEtaGaps][16]; // runId vs pt
//  TProfile2D *p_v3run_EP[3][NEtaGaps][16]; // runId vs pt
//  TProfile2D *p_v2run_BBC[3][16];
//  TProfile2D *p_v2run_ZDC[3][16];
//  TProfile2D *p_v2run_ZDC_East[3][16];
//  TProfile2D *p_v2run_ZDC_West[3][16];

  TProfile2D *p_v2_sys_Nhits_EP[3][NEtaGaps][16][2]; // runId vs pt
  TProfile2D *p_v3_sys_Nhits_EP[3][NEtaGaps][16][2]; // runId vs pt
  TProfile2D *p_v2_sys_DCA_EP[3][NEtaGaps][16][2]; // runId vs pt
  TProfile2D *p_v3_sys_DCA_EP[3][NEtaGaps][16][2]; // runId vs pt
  TProfile2D *p_v2_sys_Eta_EP[3][NEtaGaps][16]; // runId vs pt
  TProfile2D *p_v3_sys_Eta_EP[3][NEtaGaps][16]; // runId vs pt

  for (int iPID = 0; iPID < 3; iPID++)
  {
    for (int iGap = 0; iGap < NEtaGaps; iGap++)
    {
      for (int iCent = 0; iCent < 16; iCent++)
      {
        p_v2_sys_Eta_EP[iPID][iGap][iCent] = new TProfile2D(Form("p_v2_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
                                                       Form("p_v2_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
                                                       fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
        p_v3_sys_Eta_EP[iPID][iGap][iCent] = new TProfile2D(Form("p_v3_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
                                                       Form("p_v3_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
                                                       fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
        for (int iCut = 0; iCut < 2; iCut++)
        {
          p_v2_sys_Nhits_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
              Form("p_v2_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
              Form("p_v2_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
              runId_max + 0.5, 100, 0.15, 5.15);
          p_v3_sys_Nhits_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
              Form("p_v3_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
              Form("p_v3_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
              runId_max + 0.5, 100, 0.15, 5.15);

          p_v2_sys_DCA_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
              Form("p_v2_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
              Form("p_v2_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
              runId_max + 0.5, 100, 0.15, 5.15);
          p_v3_sys_DCA_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
              Form("p_v3_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
              Form("p_v3_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
              runId_max + 0.5, 100, 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 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;

  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 Psi_ReCenter, dPsi;
  Int_t bin;
  Double_t res2EP, res3EP, res2EPsqr, res3EPsqr;
  Double_t flow2EP;
  Double_t flow3EP;
  Int_t i_part;

  Double_t flow2EPBBC, res2EPBBC;
  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 flow2EPZDC, res2EPZDC;
  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];

  std::pair<Double_t, Double_t> XY_pid;

  // 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.;

        //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.);

        //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.);

        //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.);

        //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.);
      }
    }

    //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;

      // SYSTEMATICS
      // Nhits & DCA variation
      for (int iCut = 0; iCut<2; iCut++)
      {
        // Nhits
        if (isGoodTrackFlowPIDsys_Nhits(femtoTrack, energy, pVtx,iCut))
          //&& isGoodPID(femtoTrack))
        {

          weight = GetWeight(femtoTrack);

          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);

            //Loop over particle species (pi,K,p)
            i_part = -1;
            //TPC-only

            if (!femtoTrack->isTofTrack())
            {
              //pion id
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.6 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 2)
              {
                i_part = 0;
              }
              // kaon id
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.5 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 2)
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 0.9 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 2)
              {
                i_part = 2;
              }
              //pion id
              if (femtoTrack->ptot() >= 0.6 && femtoTrack->ptot() < 0.7 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) > 2)
              {
                i_part = 0;
              }
              // kaon id
              if (femtoTrack->ptot() >= 0.5 && femtoTrack->ptot() < 0.7 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) > 3)
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.9 && femtoTrack->ptot() < 1.2 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) > 3)
              {
                i_part = 2;
              }
            }

            //TPC+TOF

            if (isGoodPID(femtoTrack))
            {
              // pion id - asymmetry cut
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
                  femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)// &&
              {
                i_part = 0;
              }
              // kaon id - asymmetry cut
              if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
                  femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)// &&
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
                  femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)// &&
              {
                i_part = 2;
              }

            }
            //}

            if (i_part == -1) continue;

            //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_v2_sys_Nhits_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow2EP,
                                                                      weight);
              p_v3_sys_Nhits_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow3EP,
                                                                      weight);
            }
            //WEST
            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_v2_sys_Nhits_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow2EP,
                                                                      weight);
              p_v3_sys_Nhits_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow3EP,
                                                                      weight);
            }
          }
        }

        // DCA
        if (isGoodTrackFlowPIDsys_DCA(femtoTrack, energy, pVtx,iCut))
          //&& isGoodPID(femtoTrack))
        {

          weight = GetWeight(femtoTrack);

          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);

            //Loop over particle species (pi,K,p)
            i_part = -1;
            //TPC-only

            if (!femtoTrack->isTofTrack())
            {
              //pion id
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.6 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 2)
              {
                i_part = 0;
              }
              // kaon id
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.5 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 2)
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 0.9 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 2)
              {
                i_part = 2;
              }
              //pion id
              if (femtoTrack->ptot() >= 0.6 && femtoTrack->ptot() < 0.7 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) > 2)
              {
                i_part = 0;
              }
              // kaon id
              if (femtoTrack->ptot() >= 0.5 && femtoTrack->ptot() < 0.7 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) > 3)
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.9 && femtoTrack->ptot() < 1.2 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) > 3)
              {
                i_part = 2;
              }
            }

            //TPC+TOF

            if (isGoodPID(femtoTrack))
            {
              // pion id - asymmetry cut
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
                  femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)// &&
              {
                i_part = 0;
              }
              // kaon id - asymmetry cut
              if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
                  femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)// &&
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
                  femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)// &&
              {
                i_part = 2;
              }

            }
            //}

            if (i_part == -1) continue;

            //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_v2_sys_DCA_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow2EP,
                                                                      weight);
              p_v3_sys_DCA_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow3EP,
                                                                      weight);
            }
            //WEST
            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_v2_sys_DCA_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow2EP,
                                                                      weight);
              p_v3_sys_DCA_EP[i_part][iGap][(Int_t) event->cent16()][iCut]->Fill(event->runId(), femtoTrack->pt(), flow3EP,
                                                                      weight);
            }
          }
        }
      }

      // Eta variation
      if (isGoodTrackFlowPIDsys_Eta(femtoTrack, energy, pVtx,0))
          //&& isGoodPID(femtoTrack))
        {

          weight = GetWeight(femtoTrack);

          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);

            //Loop over particle species (pi,K,p)
            i_part = -1;
            //TPC-only

            if (!femtoTrack->isTofTrack())
            {
              //pion id
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.6 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 2)
              {
                i_part = 0;
              }
              // kaon id
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.5 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 2)
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 0.9 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 2)
              {
                i_part = 2;
              }
              //pion id
              if (femtoTrack->ptot() >= 0.6 && femtoTrack->ptot() < 0.7 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) > 2)
              {
                i_part = 0;
              }
              // kaon id
              if (femtoTrack->ptot() >= 0.5 && femtoTrack->ptot() < 0.7 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) > 3)
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.9 && femtoTrack->ptot() < 1.2 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 2 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) > 3)
              {
                i_part = 2;
              }
            }

            //TPC+TOF

            if (isGoodPID(femtoTrack))
            {
              // pion id - asymmetry cut
              if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
                  femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)// &&
              {
                i_part = 0;
              }
              // kaon id - asymmetry cut
              if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
                  femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)// &&
              {
                i_part = 1;
              }
              // proton id
              if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
                  TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
                  femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)// &&
              {
                i_part = 2;
              }

            }
            //}

            if (i_part == -1) continue;

            //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_v2_sys_Eta_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP,
                                                                      weight);
              p_v3_sys_Eta_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP,
                                                                      weight);
            }
            //WEST
            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_v2_sys_Eta_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP,
                                                                      weight);
              p_v3_sys_Eta_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP,
                                                                      weight);
            }
          }
        }

    }

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

  femtoReader->Finish();

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

  for (int iPID = 0; iPID < 3; iPID++)
  {
    for (int iGap = 0; iGap < NEtaGaps; iGap++)
    {
      for (int iCent = 0; iCent < 16; iCent++)
      {
        p_v2_sys_Eta_EP[iPID][iGap][iCent]->Write();
        p_v3_sys_Eta_EP[iPID][iGap][iCent]->Write();
        for (int iCut=0; iCut<2; iCut++)
        {
          p_v2_sys_Nhits_EP[iPID][iGap][iCent][iCut]->Write();
          p_v3_sys_Nhits_EP[iPID][iGap][iCent][iCut]->Write();
          p_v2_sys_DCA_EP[iPID][iGap][iCent][iCut]->Write();
          p_v3_sys_DCA_EP[iPID][iGap][iCent][iCut]->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 isGoodTrackFlowPID(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() > 4.15) return false;
  if (track->ptot() > cutPMax) return false;
  if (track->gDCA(pVtx).Mag() >= cutDCA_PID) return false;
  return true;
}

Bool_t isGoodTrackFlowPIDsys_Eta(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount)
{
  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()) >= cutEtasys.at(sysCount)) return false;

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

Bool_t isGoodTrackFlowPIDsys_Nhits(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount)
{
  if (!track) return false;
  // if (!track->isPrimary()) return false;
  if (track->nHitsFit() < cutNhitssys.at(sysCount)) return false;
  if (track->nHitsPoss() <= cutNhitsPoss) return false;
  if (TMath::Abs(track->eta()) >= cutEta) return false;

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

Bool_t isGoodTrackFlowPIDsys_DCA(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount)
{
  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() > 4.15) return false;
  if (track->ptot() > cutPMax) return false;
  if (track->gDCA(pVtx).Mag() >= cutDCA_PIDsys.at(sysCount)) return false;
  return true;
}

Bool_t isGoodPID(StFemtoTrack *const &track)
{
  if (!track->isTofTrack()) return false;
  if (track->massSqr() < cutMass2Min) return false;
  //NhitsDedx cut applied in StFemtoDst?
  // ToFYLocal cut applied in StFemtoDst?
  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;
}

//Int_t DoPID(StFemtoTrack *const &track)
//{
//  if (track->isTofTrack())
//  {
//    if (track->ptot() >= 0.2 && track->ptot() < 3.4 &&
//        TMath::Abs(track->nSigmaPion()) < 3. &&
//	track->massSqr() > pionMassSqr - NpidFactor*sigMsqrPion->Eval(femtoTrack->pt()) &&
//	femtoTrack->massSqr() < pionMassSqr + NpidFactor*sigMsqrPion->Eval(femtoTrack->pt()))

//  }
//}

Double_t GetSigmWidth(StFemtoTrack *const &track, Int_t PID)
{
  Double_t pol0, pol1, pol2, pol3;
  Double_t pt = track->pt();
  if (PID != 0 && PID != 1 && PID != 2) return 0.;
  if (PID == 0) // pion
  {
    pol0 = -0.036277878;
    pol1 = 0.10230040;
    pol2 = -0.069674611;
    pol3 = 0.020434016;
  }
  if (PID == 1) // kaon
  {
    pol0 = 0.011004116;
    pol1 = -0.022878517;
    pol2 = 0.038562140;
    pol3 = -0.0060050387;
  }
  if (PID == 2) // proton
  {
    pol0 = 0.035215709;
    pol1 = -0.038627269;
    pol2 = 0.046936935;
    pol3 = -0.0063598115;
  }

  return (pol0 + pol1 * pt + pol2 * pt * pt + pol3 * pt * pt * pt);
}

Double_t GetMeanNsig(StFemtoTrack *const &track, const Int_t PID)
{
  Double_t pt, par0, par1, par2, par3;
  if (PID == 0) return 0.;
  if (PID == 1)
  {
    par0 = 5.84263e+00;
    par1 = -6.65912e+00;
    par2 = 8.32204e-01;
    par3 = 8.28440e-01;
  }
  if (PID == 2)
  {
    par0 = 1.09685e+01;
    par1 = -8.61499e+00;
    par2 = 8.92938e-01;
    par3 = 8.08397e-01;
  }
  pt = track->pt();
  if (pt != 0)
  {
    return (par0 * TMath::Power(1. / pt, par2) + par1 + par3 * pt);
  } else
  {
    return 0.;
  }
}

Double_t GetWidthM2(StFemtoTrack *const &track, const Int_t PID)
{
  Double_t pt, pol0, pol1, pol2;
  if (PID != 0) return 0.;
  pol0 = 8.74975e-04;
  pol1 = -1.62659e-03;
  pol2 = 2.89828e-02;

  pt = track->pt();
  return (pol0 + pol1 * pt + pol2 * pt * pt);
}

Double_t GetScaleFactor(StFemtoTrack *const &track)
{
  Double_t wM2pion = GetWidthM2(track, 0);
  Double_t exScale = 1. / 0.8273; // = 1. for BES energies
  if (wM2pion == 0) return 1.;
  else return (exScale / wM2pion);
}

Double_t GetRotAngle(StFemtoTrack *const &track)
{
  Double_t scaleFactor = GetScaleFactor(track);
  Double_t MeanKaonNsig = GetMeanNsig(track, 1);

  Double_t new_x = (MeanKaonNsig - 0.) / scaleFactor;
  Double_t new_y = (kaonMassSqr - pionMassSqr);

  return (-TMath::ATan2(new_y, new_x));
}

std::pair<Double_t, Double_t> GetNewXY(StFemtoTrack *const &track)
{
  std::pair<Double_t, Double_t> coord;
  Double_t scaleFactor = GetScaleFactor(track);

  Double_t new_x = (track->nSigmaPion() - 0.) / scaleFactor;
  Double_t new_y = (track->massSqr() - pionMassSqr);

  Double_t rotAngle = GetRotAngle(track);
  coord.first = new_x * TMath::Cos(rotAngle) - new_y * TMath::Sin(rotAngle);
  coord.second = new_x * TMath::Sin(rotAngle) + new_y * TMath::Cos(rotAngle);

  return coord;
}

Bool_t isMesonXY(Double_t x, Double_t y)
// Used to exclude proton peak from new XY coord distribution
{
  Double_t pol0, pol1;
  pol0 = -0.92857143;
  pol1 = 1.4285714;

  // Not a proton
  if (y > (pol0 + pol1 * x)) return true;
  // Proton
  if (y <= (pol0 + pol1 * x)) return false;

  // just in case...
  return false;
}

Double_t GetSigmNewX(StFemtoTrack *const &track, Int_t PID)
{
  if (PID != 0 && PID != 1 && PID != 2) return 0.;
  Double_t pol0, pol1, pol2, pol3;
  Double_t pt = track->pt();
  if (PID == 0) // pion
  {
    pol0 = 0.00061920122;
    pol1 = 0.00064431355;
    pol2 = 0.013672155;
    pol3 = 0.0012772833;
  }
  if (PID == 1) // kaon
  {
    pol0 = 0.0057203659;
    pol1 = 0.0013767144;
    pol2 = 0.0063761902;
    pol3 = 0.0073111853;
  }
  return (pol0 + pol1 * pt + pol2 * pt * pt + pol3 * pt * pt * pt);
}

Double_t GetSigmNewY(StFemtoTrack *const &track, Int_t PID)
{
  if (PID != 0 && PID != 1 && PID != 2) return 0.;
  Double_t pol0, pol1, pol2, pol3;
  Double_t pt = track->pt();
  if (PID == 0) // pion
  {
    pol0 = 0.0017345792;
    pol1 = -0.0059824003;
    pol2 = 0.029537517;
    pol3 = -0.0037990834;
  }
  if (PID == 1) // kaon
  {
    pol0 = 0.0041058086;
    pol1 = -0.0041708451;
    pol2 = 0.026109827;
    pol3 = -0.0018772891;
  }
  return (pol0 + pol1 * pt + pol2 * pt * pt + pol3 * pt * pt * pt);
}


//--------------------------------------------------------------------------------------------------------------------------//
//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;
}

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.;
}