StFlowMephi/FemtoDstAnalyzer_Resolution.C

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

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

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

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

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

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

#include "Constants.h"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

  //Manage output TProfiles
  TProfile *p_res2_EP[NEtaGaps];
  TProfile *p_res3_EP[NEtaGaps];
  TProfile *p_res2_SP[NEtaGaps];
  TProfile *p_res3_SP[NEtaGaps];
  TProfile *p_res2_SP_Norm[NEtaGaps];
  TProfile *p_res3_SP_Norm[NEtaGaps];

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

  p_res2_vs_run = new TProfile2D(Form("p_res2_vs_run"),Form("p_res2_vs_run"),fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  p_res3_vs_run = new TProfile2D(Form("p_res3_vs_run"),Form("p_res3_vs_run"),fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);

  // 3-sub method
  //
  // A - TPC-full, B - BBC-East, C - BBC-West
  TProfile *p_res2_EP_3sub_AB = new TProfile(Form("p_res2_EP_3sub_AB"), Form("p_res2_EP_3sub_AB"), 16, -0.5, 15.5);
  TProfile *p_res2_EP_3sub_AC = new TProfile(Form("p_res2_EP_3sub_AC"), Form("p_res2_EP_3sub_AC"), 16, -0.5, 15.5);
  TProfile *p_res2_EP_3sub_BC = new TProfile(Form("p_res2_EP_3sub_BC"), Form("p_res2_EP_3sub_BC"), 16, -0.5, 15.5);

  TProfile *p_res3_EP_3sub_AB = new TProfile(Form("p_res3_EP_3sub_AB"), Form("p_res3_EP_3sub_AB"), 16, -0.5, 15.5);
  TProfile *p_res3_EP_3sub_AC = new TProfile(Form("p_res3_EP_3sub_AC"), Form("p_res3_EP_3sub_AC"), 16, -0.5, 15.5);
  TProfile *p_res3_EP_3sub_BC = new TProfile(Form("p_res3_EP_3sub_BC"), Form("p_res3_EP_3sub_BC"), 16, -0.5, 15.5);

  // A - BBC-full, B - TPC-East, C - TPC-West
  TProfile *p_res2_EP_3sub1_AB = new TProfile(Form("p_res2_EP_3sub1_AB"), Form("p_res2_EP_3sub1_AB"), 16, -0.5, 15.5);
  TProfile *p_res2_EP_3sub1_AC = new TProfile(Form("p_res2_EP_3sub1_AC"), Form("p_res2_EP_3sub1_AC"), 16, -0.5, 15.5);
  TProfile *p_res2_EP_3sub1_BC = new TProfile(Form("p_res2_EP_3sub1_BC"), Form("p_res2_EP_3sub1_BC"), 16, -0.5, 15.5);

  TProfile *p_res3_EP_3sub1_AB = new TProfile(Form("p_res3_EP_3sub1_AB"), Form("p_res3_EP_3sub1_AB"), 16, -0.5, 15.5);
  TProfile *p_res3_EP_3sub1_AC = new TProfile(Form("p_res3_EP_3sub1_AC"), Form("p_res3_EP_3sub1_AC"), 16, -0.5, 15.5);
  TProfile *p_res3_EP_3sub1_BC = new TProfile(Form("p_res3_EP_3sub1_BC"), Form("p_res3_EP_3sub1_BC"), 16, -0.5, 15.5);

  TProfile *p_resBBC_EP[fNharmonics+1];
  TProfile *p_resBBC_EP_firstHarm[fNharmonics+1];
  TProfile *p_resTPCBBCEast = new TProfile("p_resTPCBBCEast","p_resTPCBBCEast",16,-0.5,15.5);
  TProfile *p_resTPCBBCWest = new TProfile("p_resTPCBBCWest","p_resTPCBBCWest",16,-0.5,15.5);
  TProfile *p_resTPCWestBBCEast = new TProfile("p_resTPCWestBBCEast","p_resTPCWestBBCEast",16,-0.5,15.5);
  TProfile *p_resTPCEastBBCWest = new TProfile("p_resTPCEastBBCWest","p_resTPCEastBBCWest",16,-0.5,15.5);

  TProfile *p_resZDC_EP[fNharmonics+1];     // ZDC E,W
  TProfile *p_resZDC1_3sub[fNharmonics+1];  // ZDC E,W & TPC
  TProfile *p_resZDC2_3sub[fNharmonics+1];  // TPC E,W & ZDC

  for (int iGap = 0; iGap < NEtaGaps; iGap++)
  {
    p_res2_EP[iGap] = new TProfile(Form("p_res2_EP_gap%i", iGap), Form("p_res2_EP_gap%i", iGap), 16, -0.5, 15.5);
    p_res3_EP[iGap] = new TProfile(Form("p_res3_EP_gap%i", iGap), Form("p_res3_EP_gap%i", iGap), 16, -0.5, 15.5);
    p_res2_SP[iGap] = new TProfile(Form("p_res2_SP_gap%i", iGap), Form("p_res2_SP_gap%i", iGap), 16, -0.5, 15.5);
    p_res3_SP[iGap] = new TProfile(Form("p_res3_SP_gap%i", iGap), Form("p_res3_SP_gap%i", iGap), 16, -0.5, 15.5);
    p_res2_SP_Norm[iGap] = new TProfile(Form("p_res2_SP_Norm_gap%i", iGap), Form("p_res2_SP_Norm_gap%i", iGap), 16, -0.5, 15.5);
    p_res3_SP_Norm[iGap] = new TProfile(Form("p_res3_SP_Norm_gap%i", iGap), Form("p_res3_SP_Norm_gap%i", iGap), 16, -0.5, 15.5);
  }
  for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  {
    p_resBBC_EP[iHarm] = new TProfile(Form("p_resBBC_EP%i", iHarm), Form("p_resBBC_EP%i", iHarm), 16, -0.5, 15.5);
    p_resBBC_EP_firstHarm[iHarm] = new TProfile(Form("p_resBBC_EP_firstHarm%i", iHarm), Form("p_resBBC_EP_firstHarm%i", iHarm), 16, -0.5, 15.5);

    p_resZDC_EP[iHarm] = new TProfile(Form("p_resZDC_EP%i", iHarm), Form("p_resZDC_EP%i", iHarm), 16, -0.5, 15.5);
    p_resZDC1_3sub[iHarm] = new TProfile(Form("p_resZDC1_3sub%i", iHarm), Form("p_resZDC1_3sub%i", iHarm), 16, -0.5, 15.5);
    p_resZDC2_3sub[iHarm] = new TProfile(Form("p_resZDC2_3sub%i", iHarm), Form("p_resZDC2_3sub%i", iHarm), 16, -0.5, 15.5);
  }

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

  TVector2 fQv2WestSP_Shift[NEtaGaps];
  TVector2 fQv2EastSP_Shift[NEtaGaps];
  TVector2 fQv3WestSP_Shift[NEtaGaps];
  TVector2 fQv3EastSP_Shift[NEtaGaps];

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

  Long64_t goodEventCounter = 0;

  Double_t Psi2FullEP, Psi3FullEP;
  Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
  Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
  Double_t Psi2EastSP[NEtaGaps], Psi3EastSP[NEtaGaps];
  Double_t Psi2WestSP[NEtaGaps], Psi3WestSP[NEtaGaps];
  Double_t Psi_ReCenter, dPsi, psiEP_BBC_Full[fNharmonics+1];
  Int_t bin;
  Double_t res2EP[NEtaGaps], res3EP[NEtaGaps];
  Double_t res2SP[NEtaGaps], res3SP[NEtaGaps];

  Double_t res2EP_3subAB, res2EP_3subAC, res2EP_3subBC;
  Double_t res3EP_3subAB, res3EP_3subAC, res3EP_3subBC;

  Double_t res2EP_3sub1AB, res2EP_3sub1AC, res2EP_3sub1BC;
  Double_t res3EP_3sub1AB, res3EP_3sub1AC, res3EP_3sub1BC;

  // BBC related EP determination
  Double_t qxEastEP[fNharmonics + 1], qyEastEP[fNharmonics + 1], qxWestEP[fNharmonics + 1], qyWestEP[fNharmonics + 1];
  Double_t adcEast, adcWest, adcSumEast, adcSumWest, adcNormEastInner, adcNormWestInner, adcNormEastOuter, adcNormWestOuter;
  Int_t adcBin, shiftBin;
  Double_t psiEP_rec, psiEP_BBC_East[fNharmonics+1], psiEP_BBC_West[fNharmonics+1], shiftCos, shiftSin;
  Double_t resEastWestBBC[fNharmonics + 1], resTPCEastBBC[fNharmonics + 1], resTPCWestBBC[fNharmonics + 1], resTPCWestBBCEast[fNharmonics+1], resTPCEastBBCWest[fNharmonics+1];
  Double_t resEastWestBBCFirstHarm[fNharmonics + 1];

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

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

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

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

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

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

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

    TVector3 pVtx = event->primaryVertex();

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

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

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

      fQv2WestSP_Shift[i].Set(0., 0.);
      fQv2EastSP_Shift[i].Set(0., 0.);
      fQv3WestSP_Shift[i].Set(0., 0.);
      fQv3EastSP_Shift[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 psiWest[fNharmonics+1];
    Double_t psiEast[fNharmonics+1];
    Double_t psiFull[fNharmonics+1];
    TVector2 qGainEast[fNharmonics+1], qGainWest[fNharmonics+1], qGainFull[fNharmonics+1];
    TVector2 qRecEast[fNharmonics+1], qRecWest[fNharmonics+1], qRecFull[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;

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

    for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
    {
      resEastWestBBC[iHarm] = TMath::Cos( (iHarm+1.) * AngleShift(psiEP_BBC_East[iHarm]-psiEP_BBC_West[iHarm],(iHarm+1.)) );
      // resEastWestBBC[iHarm] = TMath::Cos( (iHarm+1.) * (qRecEast[iHarm].Phi()-qRecWest[iHarm].Phi()) );
      // resEastWestBBC[iHarm] = TMath::Cos( (iHarm+1.) * (qGainEast[iHarm].Phi()-qGainWest[iHarm].Phi()) );
      if (iHarm==0)
      {
        resEastWestBBCFirstHarm[0] = TMath::Cos( 1. * (psiEP_BBC_East[0]-psiEP_BBC_West[0]) );
        resEastWestBBCFirstHarm[1] = TMath::Cos( 2. * (psiEP_BBC_East[0]-psiEP_BBC_West[0]) );
        resEastWestBBCFirstHarm[2] = TMath::Cos( 3. * (psiEP_BBC_East[0]-psiEP_BBC_West[0]) );
      }
    }

    // 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.);
    }
    for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
    {
      resEastWestZDC[iHarm] = TMath::Cos( (iHarm+1.) * AngleShift(psiEP_ZDC_East[iHarm]-psiEP_ZDC_West[iHarm],(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() < 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() > 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.);

      //FULL 3-nd order
      Psi_ReCenter = TMath::ATan2(fQv3FullEP.Y(), fQv3FullEP.X()) / 3.;
      dPsi = 0.;
      for (int iK = 0; iK < NShiftOrderMax; iK++)
      {
        mCos = p_Cos3_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
        mSin = p_Cos3_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(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));
      }
      Psi3FullEP = AngleShift(Psi_ReCenter + dPsi, 3.);
    }
    for (int iGap = 0; iGap < NEtaGaps; iGap++)
    {
      if (fQcountEast[iGap] > cutNcountQvGap &&
          fQcountWest[iGap] > cutNcountQvGap)
      {
        Psi2EastEP[iGap] = 0.;
        Psi3EastEP[iGap] = 0.;
        Psi2WestEP[iGap] = 0.;
        Psi3WestEP[iGap] = 0.;
        Psi2EastSP[iGap] = 0.;
        Psi3EastSP[iGap] = 0.;
        Psi2WestSP[iGap] = 0.;
        Psi3WestSP[iGap] = 0.;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    //RESOLUTION CALCULATIONS
    if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
    {
      res2EP_3subAB = TMath::Cos(2 * (Psi2FullEP - psiEP_BBC_East[1]));
      res2EP_3subAC = TMath::Cos(2 * (Psi2FullEP - psiEP_BBC_West[1]));
      res2EP_3subBC = resEastWestBBC[1];

      res2EP_3sub1AB = TMath::Cos(2 * (Psi2EastEP[0] - psiEP_BBC_Full[1]));
      res2EP_3sub1AC = TMath::Cos(2 * (Psi2WestEP[0] - psiEP_BBC_Full[1]));
      res2EP_3sub1BC = TMath::Cos(2 * (Psi2EastEP[0] - Psi2WestEP[0]));

      res3EP_3subAB = TMath::Cos(3 * (Psi3FullEP - psiEP_BBC_East[2]));
      res3EP_3subAC = TMath::Cos(3 * (Psi3FullEP - psiEP_BBC_West[2]));
      res3EP_3subBC = resEastWestBBC[2];

      res3EP_3sub1AB = TMath::Cos(3 * (Psi3EastEP[0] - psiEP_BBC_Full[2]));
      res3EP_3sub1AC = TMath::Cos(3 * (Psi3WestEP[0] - psiEP_BBC_Full[2]));
      res3EP_3sub1BC = TMath::Cos(3 * (Psi3EastEP[0] - Psi3WestEP[0]));

      p_res2_EP_3sub_AB->Fill(event->cent16(), res2EP_3subAB);
      p_res2_EP_3sub_AC->Fill(event->cent16(), res2EP_3subAC);
      p_res2_EP_3sub_BC->Fill(event->cent16(), res2EP_3subBC);

      p_res2_EP_3sub1_AB->Fill(event->cent16(), res2EP_3sub1AB);
      p_res2_EP_3sub1_AC->Fill(event->cent16(), res2EP_3sub1AC);
      p_res2_EP_3sub1_BC->Fill(event->cent16(), res2EP_3sub1BC);

      p_res3_EP_3sub_AB->Fill(event->cent16(), res3EP_3subAB);
      p_res3_EP_3sub_AC->Fill(event->cent16(), res3EP_3subAC);
      p_res3_EP_3sub_BC->Fill(event->cent16(), res3EP_3subBC);

      p_res3_EP_3sub1_AB->Fill(event->cent16(), res3EP_3sub1AB);
      p_res3_EP_3sub1_AC->Fill(event->cent16(), res3EP_3sub1AC);
      p_res3_EP_3sub1_BC->Fill(event->cent16(), res3EP_3sub1BC);
    }
    for (int iGap = 0; iGap < NEtaGaps; iGap++)
    {
      if (fQcountEast[iGap] > cutNcountQvGap &&
          fQcountWest[iGap] > cutNcountQvGap)
      {
        res2EP[iGap] = TMath::Cos(2 * (Psi2WestEP[iGap] - Psi2EastEP[iGap]));
        res3EP[iGap] = TMath::Cos(3 * (Psi3WestEP[iGap] - Psi3EastEP[iGap]));

        res2SP[iGap] = (Double_t) (fQv2WestSP_Shift[iGap] * fQv2EastSP_Shift[iGap]);
        res3SP[iGap] = (Double_t) (fQv3WestSP_Shift[iGap] * fQv3EastSP_Shift[iGap]);

        p_res2_EP[iGap]->Fill(event->cent16(), res2EP[iGap]);
        p_res3_EP[iGap]->Fill(event->cent16(), res3EP[iGap]);

        p_res2_SP[iGap]->Fill(event->cent16(), res2SP[iGap]);
        p_res3_SP[iGap]->Fill(event->cent16(), res3SP[iGap]);
        p_res2_SP_Norm[iGap]->Fill(event->cent16(), fQv2WestSP_Shift[iGap].Mod() * fQv2EastSP_Shift[iGap].Mod());
        p_res3_SP_Norm[iGap]->Fill(event->cent16(), fQv3WestSP_Shift[iGap].Mod() * fQv3EastSP_Shift[iGap].Mod());

        if (iGap == 0)
        {
          p_res2_vs_run->Fill(event->runId(),event->cent16(),res2EP[iGap]);
          p_res3_vs_run->Fill(event->runId(),event->cent16(),res3EP[iGap]);
        }
        
        if (iGap == 0)
        {
          resTPCEastBBC[1] = TMath::Cos( (1.+1.) * (psiEP_BBC_East[1]-Psi2EastEP[iGap]) );
          resTPCWestBBC[1] = TMath::Cos( (1.+1.) * (psiEP_BBC_West[1]-Psi2WestEP[iGap]) );
          resTPCEastBBCWest[1] = TMath::Cos( (1.+1.) * (psiEP_BBC_West[1]-Psi2EastEP[iGap]) );
          resTPCWestBBCEast[1] = TMath::Cos( (1.+1.) * (psiEP_BBC_East[1]-Psi2WestEP[iGap]) );

          p_resTPCBBCEast->Fill(event->cent16(),resTPCEastBBC[1]);
          p_resTPCBBCWest->Fill(event->cent16(),resTPCWestBBC[1]);
          p_resTPCWestBBCEast->Fill(event->cent16(),resTPCWestBBCEast[1]);
          p_resTPCEastBBCWest->Fill(event->cent16(),resTPCEastBBCWest[1]);
        }
      }
    }
    for (int iHarm=0;iHarm<fNharmonics+1;iHarm++)
    {
      p_resBBC_EP[iHarm]->Fill(event->cent16(), resEastWestBBC[iHarm]);
      p_resBBC_EP_firstHarm[iHarm]->Fill(event->cent16(), resEastWestBBCFirstHarm[iHarm]);

      p_resZDC_EP[iHarm]->Fill(event->cent16(), resEastWestZDC[iHarm]);
    }

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

  femtoReader->Finish();

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

  for (int iGap = 0; iGap < NEtaGaps; iGap++)
  {
    p_res2_EP[iGap]->Write();
    p_res3_EP[iGap]->Write();
    p_res2_SP[iGap]->Write();
    p_res3_SP[iGap]->Write();
    p_res2_SP_Norm[iGap]->Write();
    p_res3_SP_Norm[iGap]->Write();
  }
  for (int iHarm=0;iHarm<fNharmonics+1;iHarm++)
  {
    p_resBBC_EP[iHarm]->Write();
    p_resBBC_EP_firstHarm[iHarm]->Write();
  }
  for (int iHarm=0;iHarm<fNharmonics+1;iHarm++)
  {
    p_resZDC_EP[iHarm]->Write();
  }
  p_resTPCBBCEast->Write();
  p_resTPCBBCWest->Write();
  p_resTPCWestBBCEast->Write();
  p_resTPCEastBBCWest->Write();

  p_res2_EP_3sub_AB->Write();
  p_res2_EP_3sub_AC->Write();
  p_res2_EP_3sub_BC->Write();

  p_res3_EP_3sub_AB->Write();
  p_res3_EP_3sub_AC->Write();
  p_res3_EP_3sub_BC->Write();

  p_res2_EP_3sub1_AB->Write();
  p_res2_EP_3sub1_AC->Write();
  p_res2_EP_3sub1_BC->Write();

  p_res3_EP_3sub1_AB->Write();
  p_res3_EP_3sub1_AC->Write();
  p_res3_EP_3sub1_BC->Write();

  p_res2_vs_run->Write();
  p_res3_vs_run->Write();

  output->Close();

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

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

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

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

  return true;
}

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

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

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

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

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

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

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

  return qv;
}

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

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

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

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

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

  return bbc_phi;
}

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