StFlowMephi/FemtoDstAnalyzer_BBC_ShiftCorr.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);
Int_t GetVzBin(Double_t vtxZ);
Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);

//_________________
void FemtoDstAnalyzer_BBC_ShiftCorr(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
                                    const Char_t *outFile = "./oShiftCorrTest.root",
                                    const Char_t *gainBBC = "",
                                    const Char_t *recBBC = "")
{

  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

  //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
  TProfile2D *p_BBC_Qx_Full_EP[fNharmonics+1][fArraySize];
  TProfile2D *p_BBC_Qy_Full_EP[fNharmonics+1][fArraySize];
  TProfile2D *p_BBC_Qx_East_EP[fNharmonics+1][fArraySize];
  TProfile2D *p_BBC_Qy_East_EP[fNharmonics+1][fArraySize];
  TProfile2D *p_BBC_Qx_West_EP[fNharmonics+1][fArraySize];
  TProfile2D *p_BBC_Qy_West_EP[fNharmonics+1][fArraySize];

  TFile *fiReCentBBC = new TFile(recBBC, "read");
  fiReCentBBC->cd();

  for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      p_BBC_Qx_Full_EP[iHarm][iVtx] = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
      p_BBC_Qy_Full_EP[iHarm][iVtx] = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
      p_BBC_Qx_East_EP[iHarm][iVtx] = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
      p_BBC_Qy_East_EP[iHarm][iVtx] = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
      p_BBC_Qx_West_EP[iHarm][iVtx] = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
      p_BBC_Qy_West_EP[iHarm][iVtx] = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
    }
  }

  //Manage output
  TProfile2D *p_BBC_Cos_Full_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  TProfile2D *p_BBC_Sin_Full_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  TProfile2D *p_BBC_Cos_East_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  TProfile2D *p_BBC_Sin_East_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  TProfile2D *p_BBC_Cos_West_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  TProfile2D *p_BBC_Sin_West_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];

  for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
      {
        p_BBC_Cos_Full_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_BBC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_BBC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_BBC_Sin_Full_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_BBC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_BBC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_BBC_Cos_East_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_BBC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_BBC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_BBC_Sin_East_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_BBC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_BBC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_BBC_Cos_West_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_BBC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_BBC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_BBC_Sin_West_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_BBC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_BBC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
      }
    }
  }

  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;

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

    VtxSign = GetVzBin(event->primaryVertex().Z());
    if (VtxSign == -1)
      continue;

    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];
    Double_t psiWest[fNharmonics+1];
    Double_t psiEast[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+1; iHarm++)
    {
      TVector2 recFullTMP((Double_t)p_BBC_Qx_Full_EP[iHarm][VtxSign]->GetBinContent(p_BBC_Qx_Full_EP[iHarm][VtxSign]->FindBin(event->runId(), event->cent16())),
                          (Double_t)p_BBC_Qy_Full_EP[iHarm][VtxSign]->GetBinContent(p_BBC_Qy_Full_EP[iHarm][VtxSign]->FindBin(event->runId(), event->cent16())));

      TVector2 recEastTMP((Double_t)p_BBC_Qx_East_EP[iHarm][VtxSign]->GetBinContent(p_BBC_Qx_East_EP[iHarm][VtxSign]->FindBin(event->runId(), event->cent16())),
                          (Double_t)p_BBC_Qy_East_EP[iHarm][VtxSign]->GetBinContent(p_BBC_Qy_East_EP[iHarm][VtxSign]->FindBin(event->runId(), event->cent16())));
      
      TVector2 recWestTMP((Double_t)p_BBC_Qx_West_EP[iHarm][VtxSign]->GetBinContent(p_BBC_Qx_West_EP[iHarm][VtxSign]->FindBin(event->runId(), event->cent16())),
                          (Double_t)p_BBC_Qy_West_EP[iHarm][VtxSign]->GetBinContent(p_BBC_Qy_West_EP[iHarm][VtxSign]->FindBin(event->runId(), event->cent16())));

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

    Double_t PsiSin = 0., PsiCos = 0., psiShiftedTMP = 0.;

    for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
    {
      for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
      {
        //Full BBC
        if (qRecFull[iHarm].Mod())
        {
          psiShiftedTMP = qRecFull[iHarm].Phi() / (iHarm + 1.);
          if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
        }
        PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
        PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
        
        p_BBC_Cos_Full_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
        p_BBC_Sin_Full_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);

        //East BBC
        if (qRecEast[iHarm].Mod())
        {
          psiShiftedTMP = qRecEast[iHarm].Phi() / (iHarm + 1.);
          if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
        }
        PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
        PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
        
        p_BBC_Cos_East_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
        p_BBC_Sin_East_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);

        //West BBC
        if (qRecWest[iHarm].Mod())
        {
          psiShiftedTMP = qRecWest[iHarm].Phi() / (iHarm + 1.);
          if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
        }
        PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
        PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
        
        p_BBC_Cos_West_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
        p_BBC_Sin_West_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
      }
    }
  } // end iEvent loop

  femtoReader->Finish();

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

  for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  {
    for (int iVtx = 0; iVtx < fArraySize; iVtx++)
    {
      for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
      {
        p_BBC_Cos_Full_EP[iHarm][iVtx][iShift]->Write();
        p_BBC_Sin_Full_EP[iHarm][iVtx][iShift]->Write();
        p_BBC_Cos_East_EP[iHarm][iVtx][iShift]->Write();
        p_BBC_Sin_East_EP[iHarm][iVtx][iShift]->Write();
        p_BBC_Cos_West_EP[iHarm][iVtx][iShift]->Write();
        p_BBC_Sin_West_EP[iHarm][iVtx][iShift]->Write();
      }
    }
  }
  output->Close();

  std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
            << 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;
}

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

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