StFlowMephi/FemtoDstAnalyzer_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"

// 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);
Double_t GetWeight(StFemtoTrack *const &track);
TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
Int_t GetVzBin(Double_t vtxZ);

//_________________
void FemtoDstAnalyzer_ShiftCorr(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
                                const Char_t *outFile = "./oShiftTest.root",
                                const Char_t *reCentFile = "")
{

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

/*  TProfile2D *p_Qx_FULL_EP[fNharmonics][fArraySize];
  TProfile2D *p_Qy_FULL_EP[fNharmonics][fArraySize];  
  TProfile2D *p_Qx_FULL_SP[fNharmonics][fArraySize];
  TProfile2D *p_Qy_FULL_SP[fNharmonics][fArraySize];
  TProfile2D *p_Qx_EAST_EP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_EAST_EP[fNharmonics][fArraySize][NEtaGaps];  
  TProfile2D *p_Qx_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qx_WEST_EP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_WEST_EP[fNharmonics][fArraySize][NEtaGaps];  
  TProfile2D *p_Qx_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
*/
  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;

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

  //Manage output
  TProfile2D *p_Cos2_FULL_EP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Sin2_FULL_EP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Cos2_FULL_SP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Sin2_FULL_SP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Cos2_EAST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin2_EAST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Cos2_EAST_SP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin2_EAST_SP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Cos2_WEST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin2_WEST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Cos2_WEST_SP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin2_WEST_SP[fArraySize][NShiftOrderMax][NEtaGaps];

  TProfile2D *p_Cos3_FULL_EP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Sin3_FULL_EP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Cos3_FULL_SP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Sin3_FULL_SP[fArraySize][NShiftOrderMax];
  TProfile2D *p_Cos3_EAST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin3_EAST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Cos3_EAST_SP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin3_EAST_SP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Cos3_WEST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin3_WEST_EP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Cos3_WEST_SP[fArraySize][NShiftOrderMax][NEtaGaps];
  TProfile2D *p_Sin3_WEST_SP[fArraySize][NShiftOrderMax][NEtaGaps];

  for (int iVtx = 0; iVtx < 10; iVtx++)
  {
    for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
    {
      p_Cos2_FULL_EP[iVtx][iShift] = new TProfile2D(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift), Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
      p_Sin2_FULL_EP[iVtx][iShift] = new TProfile2D(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift), Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
      p_Cos2_FULL_SP[iVtx][iShift] = new TProfile2D(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift), Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
      p_Sin2_FULL_SP[iVtx][iShift] = new TProfile2D(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift), Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);

        p_Cos3_FULL_EP[iVtx][iShift] = new TProfile2D(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift), Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin3_FULL_EP[iVtx][iShift] = new TProfile2D(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift), Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Cos3_FULL_SP[iVtx][iShift] = new TProfile2D(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift), Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin3_FULL_SP[iVtx][iShift] = new TProfile2D(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift), Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
      for (int iGap = 0; iGap < NEtaGaps; iGap++)
      {
        p_Cos2_EAST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        std::cout << "Init [" << iVtx << "][" << iShift << "][" << iGap << "] p_Cos2_EAST_EP[0][0][0] Name: " << p_Cos2_EAST_EP[0][0][0]->GetName() << std::endl;
        std::cout << "\t p_Cos2_EAST_EP[0][0][0] address: " << &p_Cos2_EAST_EP[0][0][0] << std::endl;
        p_Sin2_EAST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Cos2_EAST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin2_EAST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);

        p_Cos2_WEST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin2_WEST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Cos2_WEST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin2_WEST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);

        p_Cos3_EAST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin3_EAST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Cos3_EAST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin3_EAST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);

        p_Cos3_WEST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin3_WEST_EP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Cos3_WEST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
        p_Sin3_WEST_SP[iVtx][iShift][iGap] = new TProfile2D(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap), fNBins, runId_min - 0.5, runId_max + 0.5, 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;

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

  Long64_t goodEventCounter = 0;
  Long64_t events2Fill = 0;
  Long64_t events2FillArms[NEtaGaps];
  for (int i=0; i<NEtaGaps; i++)
  {
    events2FillArms[i] = 0;
  }

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

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

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

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

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


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

    TVector3 pVtx = event->primaryVertex();

    // Track analysis
    Int_t nTracks = dst->numberOfTracks();
    
    //Counters
    fUsedTracks = 0;
    fQcountFull = 0;
    fQcountFullEast = 0;
    fQcountFullWest = 0;
    psiShiftedTMP = 0.;
    PsiSin = 0.;
    PsiCos = 0.;
    fQv2FullEP.Set(0.,0.);
    fQv2FullSP.Set(0.,0.);
    fQv3FullEP.Set(0.,0.);
    fQv3FullSP.Set(0.,0.);
    QvTMP.Set(0.,0.);

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

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

    // 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_Qy_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);

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

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

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

      QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
      QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
      
      fQv2FullSP += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
      QvTMP.Set(0.,0.);

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

      fQv3FullSP += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
      QvTMP.Set(0.,0.);
      
      fQcountFull++;
      if (femtoTrack->eta() >= 0.)
      {
        fQcountFullWest++;
      }
      else
      {
        fQcountFullEast++;
      }

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

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

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

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

          QvTMP.SetX(p_Qx_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          
          fQv2EastSP[iGap] += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
          QvTMP.Set(0.,0.);

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

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

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

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

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

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

          QvTMP.SetX(p_Qx_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          QvTMP.SetY(p_Qy_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
          
          fQv2WestSP[iGap] += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
          QvTMP.Set(0.,0.);

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

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

          fQcountWest[iGap]++;
        }
      }

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

      } //if( isTofTrack() )

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


    //Filling Shift coefficents for FULL
    if (fQcountFull > cutNcountQvFull &&
      fQcountFullEast > cutNcountQvFullEast &&
      fQcountFullWest > cutNcountQvFullWest)
    {
      events2Fill++;
      for (int iK=0; iK<NShiftOrderMax; iK++)
      {
        psiShiftedTMP = TMath::ATan2(fQv2FullEP.Y(),fQv2FullEP.X())/2.;
        PsiSin = TMath::Sin(shiftOrderCoeff2.at(iK)*psiShiftedTMP);
        PsiCos = TMath::Cos(shiftOrderCoeff2.at(iK)*psiShiftedTMP);

        p_Cos2_FULL_EP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
        p_Sin2_FULL_EP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

        psiShiftedTMP = TMath::ATan2(fQv3FullEP.Y(),fQv3FullEP.X())/3.;
        PsiSin = TMath::Sin(shiftOrderCoeff3.at(iK)*psiShiftedTMP);
        PsiCos = TMath::Cos(shiftOrderCoeff3.at(iK)*psiShiftedTMP);

        p_Cos3_FULL_EP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
        p_Sin3_FULL_EP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

        psiShiftedTMP = TMath::ATan2(fQv2FullSP.Y(),fQv2FullSP.X())/2.;
        PsiSin = TMath::Sin(shiftOrderCoeff2.at(iK)*psiShiftedTMP);
        PsiCos = TMath::Cos(shiftOrderCoeff2.at(iK)*psiShiftedTMP);

        p_Cos2_FULL_SP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
        p_Sin2_FULL_SP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

        psiShiftedTMP = TMath::ATan2(fQv3FullSP.Y(),fQv3FullSP.X())/3.;
        PsiSin = TMath::Sin(shiftOrderCoeff3.at(iK)*psiShiftedTMP);
        PsiCos = TMath::Cos(shiftOrderCoeff3.at(iK)*psiShiftedTMP);

        p_Cos3_FULL_SP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
        p_Sin3_FULL_SP[VtxSign][iK]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);
      }
    }

    for (int iGap=0; iGap<NEtaGaps; iGap++)
    {
      if (fQcountEast[iGap] > cutNcountQvGap &&
          fQcountWest[iGap] > cutNcountQvGap)
      {
        events2FillArms[iGap]++;
        for (int iK=0; iK<NShiftOrderMax; iK++)
        {
          //EAST arm
          psiShiftedTMP = TMath::ATan2(fQv2EastEP[iGap].Y(),fQv2EastEP[iGap].X())/2.;
          PsiSin = TMath::Sin(shiftOrderCoeff2.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff2.at(iK)*psiShiftedTMP);

          p_Cos2_EAST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin2_EAST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          psiShiftedTMP = TMath::ATan2(fQv3EastEP[iGap].Y(),fQv3EastEP[iGap].X())/3.;
          PsiSin = TMath::Sin(shiftOrderCoeff3.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff3.at(iK)*psiShiftedTMP);

          p_Cos3_EAST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin3_EAST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          psiShiftedTMP = TMath::ATan2(fQv2EastSP[iGap].Y(),fQv2EastSP[iGap].X())/2.;
          PsiSin = TMath::Sin(shiftOrderCoeff2.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff2.at(iK)*psiShiftedTMP);

          p_Cos2_EAST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin2_EAST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          psiShiftedTMP = TMath::ATan2(fQv3EastSP[iGap].Y(),fQv3EastSP[iGap].X())/3.;
          PsiSin = TMath::Sin(shiftOrderCoeff3.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff3.at(iK)*psiShiftedTMP);

          p_Cos3_EAST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin3_EAST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          //WEST arm
          psiShiftedTMP = TMath::ATan2(fQv2WestEP[iGap].Y(),fQv2WestEP[iGap].X())/2.;
          PsiSin = TMath::Sin(shiftOrderCoeff2.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff2.at(iK)*psiShiftedTMP);

          p_Cos2_WEST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin2_WEST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          psiShiftedTMP = TMath::ATan2(fQv3WestEP[iGap].Y(),fQv3WestEP[iGap].X())/3.;
          PsiSin = TMath::Sin(shiftOrderCoeff3.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff3.at(iK)*psiShiftedTMP);

          p_Cos3_WEST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin3_WEST_EP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          psiShiftedTMP = TMath::ATan2(fQv2WestSP[iGap].Y(),fQv2WestSP[iGap].X())/2.;
          PsiSin = TMath::Sin(shiftOrderCoeff2.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff2.at(iK)*psiShiftedTMP);

          p_Cos2_WEST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin2_WEST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);

          psiShiftedTMP = TMath::ATan2(fQv3WestSP[iGap].Y(),fQv3WestSP[iGap].X())/3.;
          PsiSin = TMath::Sin(shiftOrderCoeff3.at(iK)*psiShiftedTMP);
          PsiCos = TMath::Cos(shiftOrderCoeff3.at(iK)*psiShiftedTMP);

          p_Cos3_WEST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiCos);
          p_Sin3_WEST_SP[VtxSign][iK][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),PsiSin);
        }
      }
    }

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

  femtoReader->Finish();

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

  for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  {
    for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
    {
      p_Cos2_FULL_EP[iVtx][iShift]->Write();
      p_Sin2_FULL_EP[iVtx][iShift]->Write();
      p_Cos2_FULL_SP[iVtx][iShift]->Write();
      p_Sin2_FULL_SP[iVtx][iShift]->Write();

      p_Cos3_FULL_EP[iVtx][iShift]->Write();
      p_Sin3_FULL_EP[iVtx][iShift]->Write();
      p_Cos3_FULL_SP[iVtx][iShift]->Write();
      p_Sin3_FULL_SP[iVtx][iShift]->Write();
      for (int iGap = 0; iGap < NEtaGaps; iGap++)
      {
        p_Cos2_EAST_EP[iVtx][iShift][iGap]->Write();
        p_Sin2_EAST_EP[iVtx][iShift][iGap]->Write();
        p_Cos2_EAST_SP[iVtx][iShift][iGap]->Write();
        p_Sin2_EAST_SP[iVtx][iShift][iGap]->Write();

        p_Cos2_WEST_EP[iVtx][iShift][iGap]->Write();
        p_Sin2_WEST_EP[iVtx][iShift][iGap]->Write();
        p_Cos2_WEST_SP[iVtx][iShift][iGap]->Write();
        p_Sin2_WEST_SP[iVtx][iShift][iGap]->Write();
        
        p_Cos3_EAST_EP[iVtx][iShift][iGap]->Write();
        p_Sin3_EAST_EP[iVtx][iShift][iGap]->Write();
        p_Cos3_EAST_SP[iVtx][iShift][iGap]->Write();
        p_Sin3_EAST_SP[iVtx][iShift][iGap]->Write();

        p_Cos3_WEST_EP[iVtx][iShift][iGap]->Write();
        p_Sin3_WEST_EP[iVtx][iShift][iGap]->Write();
        p_Cos3_WEST_SP[iVtx][iShift][iGap]->Write();
        p_Sin3_WEST_SP[iVtx][iShift][iGap]->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;
  std::cout << "Good events for TProfile2D fill: " << events2Fill << std::endl;
  std::cout << "Good events for TProfile2D fill arms: " << events2FillArms[0] << 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;
}