FlowMC/Flow/FlowCalculator.cxx

#ifndef FLOWCALCULATOR_CXX
#define FLOWCALCULATOR_CXX

#include "FlowCalculator.h"
#include <TMath.h>
#include <TProfile.h>

#include <math.h>

// Function that normalizes angle to the range [0,2*pi/harmonic]
Double_t NormalizeAngle(Double_t phi, Double_t harm)
{
  phi = fmod(phi,2*TMath::Pi()/harm);
  if (phi < 0)
    phi += 2*TMath::Pi()/harm;
  return phi;
}

// Wrap angle to [-pi,pi]
Double_t WrapAngle(Double_t phi, Double_t harm)
{
  phi = fmod(phi + TMath::Pi()/harm, 2*TMath::Pi()/harm);
  if (phi < 0)
    phi += 2*TMath::Pi()/harm;
  return phi - TMath::Pi();
}

FlowCalculator::FlowCalculator(TTree *tree) : BaseReader(tree)
{
    centCalc = new CentralityCalculator();
    eta_gap.push_back(0.05);
    eta_gap.push_back(0.1);
    eta_gap.push_back(0.15);
    eta_gap.push_back(0.25);
    eta_rxn.push_back({1.,2.8});
    eta_rxn.push_back({1.,1.5});
    eta_rxn.push_back({1.5,2.8});
    fRotate = true;
    fDate = new TDatime();
    fRnd = new TRandom3(fDate->GetDate()+fDate->GetTime());
    gRandom->SetSeed(fDate->GetDate()+fDate->GetTime());
    fPsiRP = 0.;
}

FlowCalculator::~FlowCalculator()
{
    delete centCalc;
    eta_gap.clear();
    eta_rxn.clear();
    delete fDate;
    delete fRnd;
}

void FlowCalculator::InitCentFile(TString inName)
{
    TFile *fi = new TFile(inName.Data(),"read");
    TH1F *hist = (TH1F*) fi->Get(histRefMultName.Data());
    centCalc->SetHistogram(hist);
    centCalc->SetEfficiency(0.93);
    centCalc->Calc();
}

void FlowCalculator::InitResFile(TString inName)
{
    TFile *fi = new TFile(inName.Data(),"read");
    for (int iGap=0; iGap<4; iGap++)
    {
        h_res2TPC[iGap] = (TH1F*) fi->Get(Form("hres2TPC_gap%i",iGap));
        h_res3TPC[iGap] = (TH1F*) fi->Get(Form("hres3TPC_gap%i",iGap));
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        h_res2RXN[irxn] = (TH1F*) fi->Get(Form("hres2RXNfull_part%i",irxn));
        h_res3RXN[irxn] = (TH1F*) fi->Get(Form("hres3RXNfull_part%i",irxn));
    }

    for (int iGap=0; iGap<4; iGap++)
    {
        for (int ibin = 0; ibin<h_res2TPC[iGap]->GetNbinsX();ibin++)
        {
            fRes2TPC[iGap][ibin] = h_res2TPC[iGap]->GetBinContent(ibin+1);
        }
        for (int ibin = 0; ibin<h_res3TPC[iGap]->GetNbinsX();ibin++)
        {
            fRes3TPC[iGap][ibin] = h_res3TPC[iGap]->GetBinContent(ibin+1);
        }
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        for (int ibin = 0; ibin<h_res2RXN[irxn]->GetNbinsX();ibin++)
        {
            fRes2RXN[irxn][ibin] = h_res2RXN[irxn]->GetBinContent(ibin+1);
            std::cout << "fRes2RXN_part" << irxn << "[" << ibin << "]\t= "
                << fRes2RXN[irxn][ibin] << std::endl;
        }
        for (int ibin = 0; ibin<h_res3RXN[irxn]->GetNbinsX();ibin++)
        {
            fRes3RXN[irxn][ibin] = h_res3RXN[irxn]->GetBinContent(ibin+1);
        }
    }
}

void FlowCalculator::LoopRes(TString outName)
{
    std::cout << "FlowCalculator::LoopRes: Processing..." << std::endl;
    std::cout << "FlowCalculator::LoopRes: RP rotation: " << fRotate << std::endl;
    if (fChain == 0) return;

    centCalc->SetModel("PHSD");
    centCalc->SetEnergy(200.);

    TProfile *p_res1ZDCEvsZDCW;
    TProfile *p_res2RXNEvsRXNW[3], *p_res3RXNEvsRXNW[3];
    TProfile *p_res2TPCEvsTPCW[4], *p_res3TPCEvsTPCW[4];
    TProfile *p_res2RXNEvsTPCM[3], *p_res3RXNEvsTPCM[3];
    TProfile *p_res2RXNWvsTPCM[3], *p_res3RXNWvsTPCM[3];

    p_res1ZDCEvsZDCW = new TProfile(Form("p_res1ZDCEvsZDCW"),Form("p_res1ZDCEvsZDCW"),10,0.,100.);
    for (int iGap=0; iGap<4; iGap++)
    {
        p_res2TPCEvsTPCW[iGap] = new TProfile(Form("p_res2TPCEvsTPCW_gap%i",iGap),Form("p_res2TPCEvsTPCW_gap%i",iGap),10,0.,100.);
        p_res3TPCEvsTPCW[iGap] = new TProfile(Form("p_res3TPCEvsTPCW_gap%i",iGap),Form("p_res3TPCEvsTPCW_gap%i",iGap),10,0.,100.);
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        p_res2RXNEvsRXNW[irxn] = new TProfile(Form("p_res2RXNEvsRXNW_%i",irxn),Form("p_res2RXNEvsRXNW_%i",irxn),10,0.,100.);
        p_res3RXNEvsRXNW[irxn] = new TProfile(Form("p_res3RXNEvsRXNW_%i",irxn),Form("p_res3RXNEvsRXNW_%i",irxn),10,0.,100.);
        
        p_res2RXNEvsTPCM[irxn] = new TProfile(Form("p_res2RXNEvsTPCM_%i",irxn),Form("p_res2RXNEvsTPCM_%i",irxn),10,0.,100.);
        p_res3RXNEvsTPCM[irxn] = new TProfile(Form("p_res3RXNEvsTPCM_%i",irxn),Form("p_res3RXNEvsTPCM_%i",irxn),10,0.,100.);
        p_res2RXNWvsTPCM[irxn] = new TProfile(Form("p_res2RXNWvsTPCM_%i",irxn),Form("p_res2RXNWvsTPCM_%i",irxn),10,0.,100.);
        p_res3RXNWvsTPCM[irxn] = new TProfile(Form("p_res3RXNWvsTPCM_%i",irxn),Form("p_res3RXNWvsTPCM_%i",irxn),10,0.,100.);
    }

		// EP QA histograms (x axis - value, y axis - centrality bin 
    TH2F *hPsiRP;                                  // PsiRP distribution
    TH2F *hQv2xRXNEast[3], *hQv2yRXNEast[3];       // Qv2 from RXN East
    TH2F *hQv3xRXNEast[3], *hQv3yRXNEast[3];       // Qv3 from RXN East
    TH2F *hPsiEP2RXNEast[3], *hPsiEP3RXNEast[3];   // PsiEP from RXN East
    TH2F *hdPsiEP2RXNEast[3], *hdPsiEP3RXNEast[3]; // PsiEP - PsiRP from RXN East

    TH2F *hQv2xRXNWest[3], *hQv2yRXNWest[3];       // Qv2 from RXN West
    TH2F *hQv3xRXNWest[3], *hQv3yRXNWest[3];       // Qv3 from RXN West
    TH2F *hPsiEP2RXNWest[3], *hPsiEP3RXNWest[3];   // PsiEP from RXN West
    TH2F *hdPsiEP2RXNWest[3], *hdPsiEP3RXNWest[3]; // PsiEP - PsiRP from RXN West

    TH2F *hQv2xRXNFull[3], *hQv2yRXNFull[3];       // Qv2 from RXN Full
    TH2F *hQv3xRXNFull[3], *hQv3yRXNFull[3];       // Qv3 from RXN Full
    TH2F *hPsiEP2RXNFull[3], *hPsiEP3RXNFull[3];   // PsiEP from RXN Full
    TH2F *hdPsiEP2RXNFull[3], *hdPsiEP3RXNFull[3]; // PsiEP - PsiRP from RXN Full


    hPsiRP = new TH2F("hPsiRP","#Psi_{RP} vs centrality bins",360,0.,2*TMath::Pi(),10,0,100);
    for (int irxn=0; irxn<3; irxn++)
    {
       hQv2xRXNEast[irxn] = new TH2F(Form("hQv2xRXNEast_part%i",irxn),Form("hQv2xRXNEast_part%i",irxn),500,-50.,50.,10,0,100);
       hQv2yRXNEast[irxn] = new TH2F(Form("hQv2yRXNEast_part%i",irxn),Form("hQv2yRXNEast_part%i",irxn),500,-50.,50.,10,0,100);
       hQv3xRXNEast[irxn] = new TH2F(Form("hQv3xRXNEast_part%i",irxn),Form("hQv3xRXNEast_part%i",irxn),500,-50.,50.,10,0,100);
       hQv3yRXNEast[irxn] = new TH2F(Form("hQv3yRXNEast_part%i",irxn),Form("hQv3yRXNEast_part%i",irxn),500,-50.,50.,10,0,100);
       hPsiEP2RXNEast[irxn]  = new TH2F(Form("hPsiEP2RXNEast_part%i",irxn),Form("hPsiEP2RXNEast_part%i",irxn),360,0.,2*TMath::Pi(),10,0,100);
       hPsiEP3RXNEast[irxn]  = new TH2F(Form("hPsiEP3RXNEast_part%i",irxn),Form("hPsiEP3RXNEast_part%i",irxn),360,0.,2*TMath::Pi(),10,0,100);
       hdPsiEP2RXNEast[irxn] = new TH2F(Form("hdPsiEP2RXNEast_part%i",irxn),Form("hdPsiEP2RXNEast_part%i",irxn),360,-TMath::Pi(),TMath::Pi(),10,0,100);
       hdPsiEP3RXNEast[irxn] = new TH2F(Form("hdPsiEP3RXNEast_part%i",irxn),Form("hdPsiEP3RXNEast_part%i",irxn),360,-TMath::Pi(),TMath::Pi(),10,0,100);

       hQv2xRXNWest[irxn] = new TH2F(Form("hQv2xRXNWest_part%i",irxn),Form("hQv2xRXNWest_part%i",irxn),500,-50.,50.,10,0,100);
       hQv2yRXNWest[irxn] = new TH2F(Form("hQv2yRXNWest_part%i",irxn),Form("hQv2yRXNWest_part%i",irxn),500,-50.,50.,10,0,100);
       hQv3xRXNWest[irxn] = new TH2F(Form("hQv3xRXNWest_part%i",irxn),Form("hQv3xRXNWest_part%i",irxn),500,-50.,50.,10,0,100);
       hQv3yRXNWest[irxn] = new TH2F(Form("hQv3yRXNWest_part%i",irxn),Form("hQv3yRXNWest_part%i",irxn),500,-50.,50.,10,0,100);
       hPsiEP2RXNWest[irxn]  = new TH2F(Form("hPsiEP2RXNWest_part%i",irxn),Form("hPsiEP2RXNWest_part%i",irxn),360,0.,2*TMath::Pi(),10,0,100);
       hPsiEP3RXNWest[irxn]  = new TH2F(Form("hPsiEP3RXNWest_part%i",irxn),Form("hPsiEP3RXNWest_part%i",irxn),360,0.,2*TMath::Pi(),10,0,100);
       hdPsiEP2RXNWest[irxn] = new TH2F(Form("hdPsiEP2RXNWest_part%i",irxn),Form("hdPsiEP2RXNWest_part%i",irxn),360,-TMath::Pi(),TMath::Pi(),10,0,100);
       hdPsiEP3RXNWest[irxn] = new TH2F(Form("hdPsiEP3RXNWest_part%i",irxn),Form("hdPsiEP3RXNWest_part%i",irxn),360,-TMath::Pi(),TMath::Pi(),10,0,100);

       hQv2xRXNFull[irxn] = new TH2F(Form("hQv2xRXNFull_part%i",irxn),Form("hQv2xRXNFull_part%i",irxn),500,-50.,50.,10,0,100);
       hQv2yRXNFull[irxn] = new TH2F(Form("hQv2yRXNFull_part%i",irxn),Form("hQv2yRXNFull_part%i",irxn),500,-50.,50.,10,0,100);
       hQv3xRXNFull[irxn] = new TH2F(Form("hQv3xRXNFull_part%i",irxn),Form("hQv3xRXNFull_part%i",irxn),500,-50.,50.,10,0,100);
       hQv3yRXNFull[irxn] = new TH2F(Form("hQv3yRXNFull_part%i",irxn),Form("hQv3yRXNFull_part%i",irxn),500,-50.,50.,10,0,100);
       hPsiEP2RXNFull[irxn]  = new TH2F(Form("hPsiEP2RXNFull_part%i",irxn),Form("hPsiEP2RXNFull_part%i",irxn),360,0.,2*TMath::Pi(),10,0,100);
       hPsiEP3RXNFull[irxn]  = new TH2F(Form("hPsiEP3RXNFull_part%i",irxn),Form("hPsiEP3RXNFull_part%i",irxn),360,0.,2*TMath::Pi(),10,0,100);
       hdPsiEP2RXNFull[irxn] = new TH2F(Form("hdPsiEP2RXNFull_part%i",irxn),Form("hdPsiEP2RXNFull_part%i",irxn),360,-TMath::Pi(),TMath::Pi(),10,0,100);
       hdPsiEP3RXNFull[irxn] = new TH2F(Form("hdPsiEP3RXNFull_part%i",irxn),Form("hdPsiEP3RXNFull_part%i",irxn),360,-TMath::Pi(),TMath::Pi(),10,0,100);
    }

    Long64_t nentries = fChain->GetEntries();
    Long64_t nbytes = 0, nb = 0;
    Int_t refMult, cent;

    // EP from TPCs
    // TVector2 Qv2EastTPC[4], Qv2WestTPC[4];
    // TVector2 Qv3EastTPC[4], Qv3WestTPC[4];

    Double_t Psi2EastTPC[4], Psi2WestTPC[4];
    Double_t Psi3EastTPC[4], Psi3WestTPC[4];

    // EP from BBCs & RXN
    // TVector2 Qv2EastRXN[3], Qv2WestRXN[3];
    // TVector2 Qv3EastRXN[3], Qv3WestRXN[3];

    Double_t Psi2EastRXN[3], Psi2WestRXN[3], Psi2FullRXN[3];
    Double_t Psi3EastRXN[3], Psi3WestRXN[3], Psi3FullRXN[3];

    Double_t Psi2MiddTPC;
    Double_t Psi3MiddTPC;

    Double_t Psi1EastZDC, Psi1WestZDC;

    for (Long64_t jentry=0; jentry<nentries;jentry++)
    {
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        if (jentry%1000 == 0) 
            std::cout << "Event [" << jentry << "/" 
                << nentries << "]" << std::endl;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
        // if (Cut(ientry) < 0) continue;

        // Centrality determination
        refMult = GetMultPHENIX();
        cent    = centCalc->GetCentrality(refMult);

        if (fRotate)
        {
          fPsiRP = fRnd->Rndm()*2*TMath::Pi();
        } 
        else 
        {
          fPsiRP = 0.;
        }

        hPsiRP->Fill(fPsiRP,cent);

        // EP determination
        GetAllQv();
  
        for (int iGap=0; iGap<4; iGap++)
        {
            // Qv2EastTPC[iGap] = GetQvTPC(-1,2,iGap);
            // Qv2WestTPC[iGap] = GetQvTPC( 1,2,iGap);
            // Qv3EastTPC[iGap] = GetQvTPC(-1,3,iGap);
            // Qv3WestTPC[iGap] = GetQvTPC( 1,3,iGap);

            Psi2EastTPC[iGap] = GetPsiEP(fQv2EastTPC[iGap],2);
            Psi2WestTPC[iGap] = GetPsiEP(fQv2WestTPC[iGap],2);
            Psi3EastTPC[iGap] = GetPsiEP(fQv3EastTPC[iGap],3);
            Psi3WestTPC[iGap] = GetPsiEP(fQv3WestTPC[iGap],3);

            p_res2TPCEvsTPCW[iGap]->Fill(cent,TMath::Cos( 2*(Psi2EastTPC[iGap] - Psi2WestTPC[iGap]) ));
            p_res3TPCEvsTPCW[iGap]->Fill(cent,TMath::Cos( 3*(Psi3EastTPC[iGap] - Psi3WestTPC[iGap]) ));
        }

        Psi2MiddTPC = GetPsiEP(fQv2MiddTPC,2);
        Psi3MiddTPC = GetPsiEP(fQv3MiddTPC,3);

        Psi1EastZDC = GetPsiEP(fQv1EastZDC,1);
        Psi1WestZDC = GetPsiEP(fQv1WestZDC,1);

        p_res1ZDCEvsZDCW->Fill(cent,TMath::Cos( 1*(Psi1EastZDC - Psi1WestZDC) ));

        for (int irxn=0; irxn<3; irxn++)
        {
            if (multWestRXN[irxn] == 0 || multEastRXN[irxn] == 0) continue;
            // Qv2EastRXN[irxn] = GetQvRXN(-1,2,irxn);
            // Qv2WestRXN[irxn] = GetQvRXN( 1,2,irxn);
            // Qv3EastRXN[irxn] = GetQvRXN(-1,3,irxn);
            // Qv3WestRXN[irxn] = GetQvRXN( 1,3,irxn);

            Psi2EastRXN[irxn] = GetPsiEP(fQv2EastRXN[irxn],2);
            Psi2WestRXN[irxn] = GetPsiEP(fQv2WestRXN[irxn],2);
            Psi3EastRXN[irxn] = GetPsiEP(fQv3EastRXN[irxn],3);
            Psi3WestRXN[irxn] = GetPsiEP(fQv3WestRXN[irxn],3);

            p_res2RXNEvsRXNW[irxn]->Fill(cent,TMath::Cos(2* (Psi2EastRXN[irxn] - Psi2WestRXN[irxn]) ));
            p_res3RXNEvsRXNW[irxn]->Fill(cent,TMath::Cos(3* (Psi3EastRXN[irxn] - Psi3WestRXN[irxn]) ));

            p_res2RXNEvsTPCM[irxn]->Fill(cent,TMath::Cos(2* (Psi2EastRXN[irxn] - Psi2MiddTPC) ));
            p_res3RXNEvsTPCM[irxn]->Fill(cent,TMath::Cos(3* (Psi3EastRXN[irxn] - Psi3MiddTPC) ));

            p_res2RXNWvsTPCM[irxn]->Fill(cent,TMath::Cos(2* (Psi2WestRXN[irxn] - Psi2MiddTPC) ));
            p_res3RXNWvsTPCM[irxn]->Fill(cent,TMath::Cos(3* (Psi3WestRXN[irxn] - Psi3MiddTPC) ));

            hQv2xRXNEast[irxn]->Fill(fQv2EastRXN[irxn].X(),cent);
            hQv2yRXNEast[irxn]->Fill(fQv2EastRXN[irxn].Y(),cent);
            hQv3xRXNEast[irxn]->Fill(fQv3EastRXN[irxn].X(),cent);
            hQv3yRXNEast[irxn]->Fill(fQv3EastRXN[irxn].Y(),cent);

            hQv2xRXNWest[irxn]->Fill(fQv2WestRXN[irxn].X(),cent);
            hQv2yRXNWest[irxn]->Fill(fQv2WestRXN[irxn].Y(),cent);
            hQv3xRXNWest[irxn]->Fill(fQv3WestRXN[irxn].X(),cent);
            hQv3yRXNWest[irxn]->Fill(fQv3WestRXN[irxn].Y(),cent);

            hPsiEP2RXNEast[irxn]->Fill(NormalizeAngle(Psi2EastRXN[irxn],2),cent);
            hPsiEP2RXNWest[irxn]->Fill(NormalizeAngle(Psi2WestRXN[irxn],2),cent);
            hPsiEP3RXNEast[irxn]->Fill(NormalizeAngle(Psi3EastRXN[irxn],3),cent);
            hPsiEP3RXNWest[irxn]->Fill(NormalizeAngle(Psi3WestRXN[irxn],3),cent);

            hdPsiEP2RXNEast[irxn]->Fill(WrapAngle(Psi2EastRXN[irxn]-fPsiRP,2),cent);
            hdPsiEP2RXNWest[irxn]->Fill(WrapAngle(Psi2WestRXN[irxn]-fPsiRP,2),cent);
            hdPsiEP3RXNEast[irxn]->Fill(WrapAngle(Psi3EastRXN[irxn]-fPsiRP,3),cent);
            hdPsiEP3RXNWest[irxn]->Fill(WrapAngle(Psi3WestRXN[irxn]-fPsiRP,3),cent);

            Psi2FullRXN[irxn] = GetPsiEP(fQv2FullRXN[irxn],2);
            Psi3FullRXN[irxn] = GetPsiEP(fQv3FullRXN[irxn],3);

            hQv2xRXNFull[irxn]->Fill(fQv2FullRXN[irxn].X(),cent);
            hQv2yRXNFull[irxn]->Fill(fQv2FullRXN[irxn].Y(),cent);
            hQv3xRXNFull[irxn]->Fill(fQv3FullRXN[irxn].X(),cent);
            hQv3yRXNFull[irxn]->Fill(fQv3FullRXN[irxn].Y(),cent);

            hPsiEP2RXNFull[irxn]->Fill(NormalizeAngle(Psi2FullRXN[irxn],2),cent);
            hPsiEP3RXNFull[irxn]->Fill(NormalizeAngle(Psi3FullRXN[irxn],3),cent);

            hdPsiEP2RXNFull[irxn]->Fill(WrapAngle(Psi2FullRXN[irxn]-fPsiRP,2),cent);
            hdPsiEP3RXNFull[irxn]->Fill(WrapAngle(Psi3FullRXN[irxn]-fPsiRP,3),cent);
        }
    }
    TFile *fo = new TFile(outName.Data(),"recreate");
    fo->cd();
    for (int iGap=0; iGap<4; iGap++)
    {
        p_res2TPCEvsTPCW[iGap]->Write();
        p_res3TPCEvsTPCW[iGap]->Write();
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        p_res2RXNEvsRXNW[irxn]->Write();
        p_res3RXNEvsRXNW[irxn]->Write();
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        p_res2RXNEvsTPCM[irxn]->Write();
        p_res3RXNEvsTPCM[irxn]->Write();
        p_res2RXNWvsTPCM[irxn]->Write();
        p_res3RXNWvsTPCM[irxn]->Write();
    }
    p_res1ZDCEvsZDCW->Write();

		fo->mkdir("EventPlane");
    fo->cd("EventPlane");

    hPsiRP->Write();
    for (int irxn=0; irxn<3; irxn++)
    {
      hQv2xRXNEast[irxn]-> Write();
      hQv2yRXNEast[irxn]-> Write();
      hQv3xRXNEast[irxn]-> Write();
      hQv3yRXNEast[irxn]-> Write();
      hQv2xRXNWest[irxn]-> Write();
      hQv2yRXNWest[irxn]-> Write();
      hQv3xRXNWest[irxn]-> Write();
      hQv3yRXNWest[irxn]-> Write();
      hQv2xRXNFull[irxn]-> Write();
      hQv2yRXNFull[irxn]-> Write();
      hQv3xRXNFull[irxn]-> Write();
      hQv3yRXNFull[irxn]-> Write();

      hPsiEP2RXNEast[irxn]-> Write();
      hPsiEP2RXNWest[irxn]-> Write();
      hPsiEP2RXNFull[irxn]-> Write();
      hPsiEP3RXNEast[irxn]-> Write();
      hPsiEP3RXNWest[irxn]-> Write();
      hPsiEP3RXNFull[irxn]-> Write();

      hdPsiEP2RXNEast[irxn]-> Write();
      hdPsiEP2RXNWest[irxn]-> Write();
      hdPsiEP2RXNFull[irxn]-> Write();
      hdPsiEP3RXNEast[irxn]-> Write();
      hdPsiEP3RXNWest[irxn]-> Write();
      hdPsiEP3RXNFull[irxn]-> Write();
    }

    fo->Close();
}

void FlowCalculator::LoopFlow(TString outName)
{
    fRotate = true;
    std::cout << "FlowCalculator::LoopFlow: Processing..." << std::endl;
    std::cout << "FlowCalculator::LoopFlow: RP rotation: " << fRotate << std::endl;
    if (fChain == 0) return;

    TProfile *p_v2TPC[4][10][4], *p_v3TPC[4][10][4];
    TProfile *p_v2RXN[3][10][4], *p_v3RXN[3][10][4];

    TH2F *hPhi[4], *hdPhi[4];                      // Phi angle of the particles, (Phi - PsiRP) of the particles
    for (int iPID=0; iPID<4; iPID++)
    {
        hPhi[iPID]   = new TH2F(Form("hPhi_pid%i",iPID),Form("#varphi vs centrality bins pid %i",iPID),360,0.,2*TMath::Pi(),10,0,100);
        hdPhi[iPID]  = new TH2F(Form("hdPhi_pid%i",iPID),Form("(#varphi - #Psi_{RP}) vs centrality bins pid %i",iPID),360,0.,2*TMath::Pi(),10,0,100);
    }
    // [4][10][4] -- EtaGaps, CentralityBins, PID
    // PID: 0 - hadrons, 1 - pions, 2 - kaons, 3 - protons
    for (int iGap=0; iGap<4; iGap++)
    {   
        for (int iCent=0; iCent<10; iCent++)
        {
            for (int iPID=0; iPID<4; iPID++)
            {
                p_v2TPC[iGap][iCent][iPID] = new TProfile(Form("p_v2TPC_gap%i_cent%i_pid%i",iGap,iCent,iPID),
                    Form("p_v2TPC_gap%i_cent%i_pid%i",iGap,iCent,iPID),
                    500,0.1,5.1);
                p_v3TPC[iGap][iCent][iPID] = new TProfile(Form("p_v3TPC_gap%i_cent%i_pid%i",iGap,iCent,iPID),
                    Form("p_v3TPC_gap%i_cent%i_pid%i",iGap,iCent,iPID),
                    500,0.1,5.1);
            }
        }
    }
    for (int irxn=0; irxn<3; irxn++)
    {   
        for (int iCent=0; iCent<10; iCent++)
        {
            for (int iPID=0; iPID<4; iPID++)
            {
                p_v2RXN[irxn][iCent][iPID] = new TProfile(Form("p_v2RXN_part%i_cent%i_pid%i",irxn,iCent,iPID),
                    Form("p_v2RXN_part%i_cent%i_pid%i",irxn,iCent,iPID),
                    500,0.1,5.1);
                p_v3RXN[irxn][iCent][iPID] = new TProfile(Form("p_v3RXN_part%i_cent%i_pid%i",irxn,iCent,iPID),
                    Form("p_v3RXN_part%i_cent%i_pid%i",irxn,iCent,iPID),
                    500,0.1,5.1);
            }
        }
    }

    Long64_t nentries = fChain->GetEntries();
    std::cout << "CHECK!!!!" << std::endl;
    std::cout << nentries << std::endl;
    Long64_t nbytes = 0, nb = 0;
    Int_t refMult, cent;

    // EP from TPCs
    // TVector2 Qv2EastTPC[4], Qv2WestTPC[4];
    // TVector2 Qv3EastTPC[4], Qv3WestTPC[4];

    Double_t Psi2EastTPC[4], Psi2WestTPC[4];
    Double_t Psi3EastTPC[4], Psi3WestTPC[4];

    // EP from BBCs & RXN
    // TVector2 Qv2EastRXN[3], Qv2WestRXN[3];
    // TVector2 Qv3EastRXN[3], Qv3WestRXN[3];

    Double_t Psi2EastRXN[3], Psi2WestRXN[3];
    Double_t Psi3EastRXN[3], Psi3WestRXN[3];
    Double_t Psi2FullRXN[3], Psi3FullRXN[3];

    Double_t weight, vn, res2TPC[4], res3TPC[4], res2RXN[3], res3RXN[3];
    Int_t pid;


    for (Long64_t jentry=0; jentry<nentries;jentry++)
    {
        Long64_t ientry = LoadTree(jentry);
        if (ientry < 0) break;
        if (jentry%1000 == 0) 
            std::cout << "Event [" << jentry << "/" 
                << nentries << "]" << std::endl;
        nb = fChain->GetEntry(jentry);   nbytes += nb;
        // if (Cut(ientry) < 0) continue;

        // Centrality determination
        refMult = GetMultPHENIX();
        cent    = centCalc->GetCentrality(refMult);

        if (fRotate)
        {
          fPsiRP = gRandom->Rndm()*2*TMath::Pi(); //fRnd->Rndm()*2*TMath::Pi();
        } 
        else 
        {
          fPsiRP = 0.;
        }

        // EP determination
        GetAllQv();

        for (int iGap=0; iGap<4; iGap++)
        {
            // Qv2EastTPC[iGap] = GetQvTPC(-1,2,iGap);
            // Qv2WestTPC[iGap] = GetQvTPC( 1,2,iGap);
            // Qv3EastTPC[iGap] = GetQvTPC(-1,3,iGap);
            // Qv3WestTPC[iGap] = GetQvTPC( 1,3,iGap);

            Psi2EastTPC[iGap] = GetPsiEP(fQv2EastTPC[iGap],2);
            Psi2WestTPC[iGap] = GetPsiEP(fQv2WestTPC[iGap],2);
            Psi3EastTPC[iGap] = GetPsiEP(fQv3EastTPC[iGap],3);
            Psi3WestTPC[iGap] = GetPsiEP(fQv3WestTPC[iGap],3);
        }

        for (int irxn=0; irxn<3; irxn++)
        {
            // Qv2EastRXN[irxn] = GetQvRXN(-1,2,irxn);
            // Qv2WestRXN[irxn] = GetQvRXN( 1,2,irxn);
            // Qv3EastRXN[irxn] = GetQvRXN(-1,3,irxn);
            // Qv3WestRXN[irxn] = GetQvRXN( 1,3,irxn);

            Psi2EastRXN[irxn] = GetPsiEP(fQv2EastRXN[irxn],2);
            Psi2WestRXN[irxn] = GetPsiEP(fQv2WestRXN[irxn],2);
            Psi2FullRXN[irxn] = GetPsiEP(fQv2FullRXN[irxn],2);
            Psi3EastRXN[irxn] = GetPsiEP(fQv3EastRXN[irxn],3);
            Psi3WestRXN[irxn] = GetPsiEP(fQv3WestRXN[irxn],3);
            Psi3FullRXN[irxn] = GetPsiEP(fQv3FullRXN[irxn],3);
        }

        // Resolution
        for (int iGap=0; iGap<4; iGap++)
        {
            res2TPC[iGap] = fRes2TPC[iGap][(int)cent/10]; //h_res2TPC[iGap]->GetBinContent(h_res2TPC[iGap]->FindBin(cent));
            res3TPC[iGap] = fRes3TPC[iGap][(int)cent/10]; //h_res3TPC[iGap]->GetBinContent(h_res3TPC[iGap]->FindBin(cent));
        }
        for (int irxn=0; irxn<3; irxn++)
        {
            res2RXN[irxn] = fRes2RXN[irxn][(int)cent/10]; //h_res2RXN[irxn]->GetBinContent(h_res2RXN[irxn]->FindBin(cent));
            res3RXN[irxn] = fRes3RXN[irxn][(int)cent/10]; //h_res3RXN[irxn]->GetBinContent(h_res3RXN[irxn]->FindBin(cent));
        }

        for (int iTr=0; iTr<nh; iTr++)
        {
            //if (jentry == 154) std::cout << "\tTrack " << iTr << " out of " << nh << std::endl;
            // if (charge[iTr] == 0) continue;
            vMom->SetPxPyPzE(momx[iTr],momy[iTr],momz[iTr],ene[iTr]);
            double phi1, phi2;
            if (jentry == 0 && iTr == 0)
            {
                phi1 = vMom->Phi()*180./TMath::Pi();
                std::cout << "PsiRP = " << fPsiRP*180./TMath::Pi() << ", track.phi_before_rotation = " << phi1;
            }
            if (fRotate)
            {
              vMom->RotateZ(fPsiRP);
            }
            if (jentry == 0 && iTr == 0)
            {
                phi2 = vMom->Phi()*180./TMath::Pi();
                std::cout << " track.phi_after_rotation = " << phi2 << std::endl;
            }
            // if (vMom->Perp() < 0.1) continue;
            if (TMath::Abs(vMom->PseudoRapidity()) > 1.) continue;

            weight = (vMom->Perp() > 2.) ? 2. : vMom->Perp();

            // PID
            pid = -1;
            if (TMath::Abs(pdg[iTr]) == 211) pid = 1;
            if (TMath::Abs(pdg[iTr]) == 321) pid = 2;
            if (TMath::Abs(pdg[iTr]) == 2212) pid = 3;

            /*for (int iGap=0; iGap<4; iGap++)
            {
                // East
                if (TMath::Abs(vMom->PseudoRapidity()) >= eta_gap.at(iGap) &&
                    TMath::Abs(vMom->PseudoRapidity()) < 0.8 &&
                    vMom->PseudoRapidity() < 0)
                {
                    if (res2TPC[iGap] != 0)
                    {
                        vn = TMath::Cos(2 * (vMom->Phi() - Psi2WestTPC[iGap]) ) / res2TPC[iGap];
                        p_v2TPC[iGap][(Int_t)(cent/10)][0]->Fill(vMom->Perp(),vn);
                        if (pid != -1)
                        {
                            p_v2TPC[iGap][(Int_t)(cent/10)][pid]->Fill(vMom->Perp(),vn);
                        }
                    }
                    if (res3TPC[iGap] != 0)
                    {
                        vn = TMath::Cos(3 * (vMom->Phi() - Psi3WestTPC[iGap]) ) / res3TPC[iGap];
                        p_v3TPC[iGap][(Int_t)(cent/10)][0]->Fill(vMom->Perp(),vn);
                        if (pid != -1)
                        {
                            p_v3TPC[iGap][(Int_t)(cent/10)][pid]->Fill(vMom->Perp(),vn);
                        }
                    }
                }
                // West
                if (TMath::Abs(vMom->PseudoRapidity()) >= eta_gap.at(iGap) &&
                    TMath::Abs(vMom->PseudoRapidity()) < 0.8 &&
                    vMom->PseudoRapidity() > 0)
                {
                    if (res2TPC[iGap] != 0)
                    {
                        vn = TMath::Cos(2 * (vMom->Phi() - Psi2EastTPC[iGap]) ) / res2TPC[iGap];
                        p_v2TPC[iGap][(Int_t)(cent/10)][0]->Fill(vMom->Perp(),vn);
                        if (pid != -1)
                        {
                            p_v2TPC[iGap][(Int_t)(cent/10)][pid]->Fill(vMom->Perp(),vn);
                        }
                    }
                    if (res3TPC[iGap] != 0)
                    {
                        vn = TMath::Cos(3 * (vMom->Phi() - Psi3EastTPC[iGap]) ) / res3TPC[iGap];
                        p_v3TPC[iGap][(Int_t)(cent/10)][0]->Fill(vMom->Perp(),vn);
                        if (pid != -1)
                        {
                            p_v3TPC[iGap][(Int_t)(cent/10)][pid]->Fill(vMom->Perp(),vn);
                        }
                    }
                }
            }*/
            // std::cout << "Check!" << std::endl;
            hPhi[0]->Fill(NormalizeAngle(vMom->Phi(),1),cent);
            hdPhi[0]->Fill(NormalizeAngle(vMom->Phi()-fPsiRP,1),cent);
            if (pid != -1)
            {
              hPhi[pid]->Fill(NormalizeAngle(vMom->Phi(),1),cent);  
              hdPhi[pid]->Fill(NormalizeAngle(vMom->Phi()-fPsiRP,1),cent);
            }
            for (int irxn=0; irxn<3; irxn++)
            {
                if (TMath::Abs(vMom->PseudoRapidity()) < 0.35)
                {
                    if (multWestRXN[irxn] == 0 || multEastRXN[irxn] == 0) continue;
                    if (res2RXN[irxn] != 0 && fQv2FullRXN[irxn].X() != -999. &&
                        fQv2FullRXN[irxn].Y() != -999.)
                    {
                        vn = TMath::Cos(2 * (vMom->Phi() - Psi2FullRXN[irxn]) ) / res2RXN[irxn];
                        p_v2RXN[irxn][(Int_t)(cent/10)][0]->Fill(vMom->Perp(),vn);
                        if (pid != -1)
                        {
                            p_v2RXN[irxn][(Int_t)(cent/10)][pid]->Fill(vMom->Perp(),vn);
                        }
                    }
                    if (res3RXN[irxn] != 0  && fQv3FullRXN[irxn].X() != -999. &&
                        fQv3FullRXN[irxn].Y() != -999.)
                    {
                        vn = TMath::Cos(3 * (vMom->Phi() - Psi3FullRXN[irxn]) ) / res3RXN[irxn];
                        p_v3RXN[irxn][(Int_t)(cent/10)][0]->Fill(vMom->Perp(),vn);
                        if (pid != -1)
                        {
                            p_v3RXN[irxn][(Int_t)(cent/10)][pid]->Fill(vMom->Perp(),vn);
                        }
                    }
                }
            }

        } // end track loop
    } // end event loop
    TFile *fo = new TFile(outName.Data(),"recreate");
    fo->cd();
    for (int iGap=0; iGap<4; iGap++)
    {   
        for (int iCent=0; iCent<10; iCent++)
        {
            for (int iPID=0; iPID<4; iPID++)
            {
                p_v2TPC[iGap][iCent][iPID]->Write();
                p_v3TPC[iGap][iCent][iPID]->Write();
            }
        }
    }
    for (int irxn=0; irxn<3; irxn++)
    {   
        for (int iCent=0; iCent<10; iCent++)
        {
            for (int iPID=0; iPID<4; iPID++)
            {
                p_v2RXN[irxn][iCent][iPID]->Write();
                p_v3RXN[irxn][iCent][iPID]->Write();
            }
        }
    }
    fo->mkdir("EventPlane");
    fo->cd("EventPlane");
    for (int iPID=0; iPID<4; iPID++)
    {
        hPhi[iPID]->Write();
        hdPhi[iPID]->Write();
    }
}

TVector2 FlowCalculator::GetQvTPC(Int_t side, Int_t harm, Int_t gap)
{
    TVector2 qv(-999.,-999.);
    if (side == 0) return qv;
    if (harm<2 || harm>3) return qv;

    Double_t qvx = 0., qvy = 0.;

    for (int iTr=0; iTr<nh; iTr++)
    {
        if (charge[iTr] == 0) continue;
        vMom->SetPxPyPzE(momx[iTr],momy[iTr],momz[iTr],ene[iTr]);

        if (vMom->Perp() < 0.1) continue;
        if (vMom->Perp() > 2.) continue;
        if (side < 0 && vMom->PseudoRapidity() > 0) continue;
        if (side > 0 && vMom->PseudoRapidity() < 0) continue;
        if (TMath::Abs(vMom->PseudoRapidity()) < eta_gap.at(gap)) continue;
        if (TMath::Abs(vMom->PseudoRapidity()) > 0.8) continue;

        qvx += vMom->Perp() * TMath::Cos( harm * vMom->Phi() );
        qvy += vMom->Perp() * TMath::Sin( harm * vMom->Phi() );
    }
    qv.Set(qvx,qvy);
    return qv;
}

TVector2 FlowCalculator::GetQvBBC(Int_t side, Int_t harm)
{
    TVector2 qv(-999.,-999.);
    if (side == 0) return qv;
    if (harm<1 || harm>2) return qv;

    Double_t qvx = 0., qvy = 0.;

    for (int iTr=0; iTr<nh; iTr++)
    {
        if (charge[iTr] == 0) continue;
        vMom->SetPxPyPzE(momx[iTr],momy[iTr],momz[iTr],ene[iTr]);

        if (vMom->Perp() < 0.1) continue;
        if (side < 0 && vMom->PseudoRapidity() > 0) continue;
        if (side > 0 && vMom->PseudoRapidity() < 0) continue;
        if (TMath::Abs(vMom->PseudoRapidity()) < 3.) continue;
        if (TMath::Abs(vMom->PseudoRapidity()) > 4.) continue;

        qvx += vMom->Perp() * TMath::Cos( harm * vMom->Phi() );
        qvy += vMom->Perp() * TMath::Sin( harm * vMom->Phi() );
    }
    qv.Set(qvx,qvy);
    return qv;
}

TVector2 FlowCalculator::GetQvRXN(Int_t side, Int_t harm, Int_t part)
{
    TVector2 qv(-999.,-999.);
    if (side == 0) return qv;
    if (harm<1 || harm>3) return qv;

    Double_t eta_min, eta_max;

    if (part == 0)
    {
        eta_min = 1.;
        eta_max = 2.8;
    }
    else
    {
        if (part == 1)
        {
            eta_min = 1.;
            eta_max = 1.5;
        }
        else
        {
            if (part == 2)
            {
                eta_min = 1.5;
                eta_max = 2.8;
            }
            else
            {
                return qv;
            }
        }
        
    }
    
    Double_t qvx = 0., qvy = 0.;

    for (int iTr=0; iTr<nh; iTr++)
    {
        if (charge[iTr] == 0) continue;
        vMom->SetPxPyPzE(momx[iTr],momy[iTr],momz[iTr],ene[iTr]);

        if (vMom->Perp() < 0.1) continue;
        if (side < 0 && vMom->PseudoRapidity() > 0) continue;
        if (side > 0 && vMom->PseudoRapidity() < 0) continue;
        if (TMath::Abs(vMom->PseudoRapidity()) < eta_min) continue;
        if (TMath::Abs(vMom->PseudoRapidity()) > eta_max) continue;

        qvx += vMom->Perp() * TMath::Cos( harm * vMom->Phi() );
        qvy += vMom->Perp() * TMath::Sin( harm * vMom->Phi() );
    }
    qv.Set(qvx,qvy);
    return qv;
}

void FlowCalculator::GetAllQv()
{
    Double_t weight;
    Double_t qv2xEastTPC[4], qv3xEastTPC[4];
    Double_t qv2xEastRXN[3], qv3xEastRXN[3];
    Double_t qv2yEastTPC[4], qv3yEastTPC[4];
    Double_t qv2yEastRXN[3], qv3yEastRXN[3];
    Double_t qv2xWestTPC[4], qv3xWestTPC[4];
    Double_t qv2xWestRXN[3], qv3xWestRXN[3];
    Double_t qv2yWestTPC[4], qv3yWestTPC[4];
    Double_t qv2yWestRXN[3], qv3yWestRXN[3];
    Double_t qv2xFullRXN[3], qv3xFullRXN[3];
    Double_t qv2yFullRXN[3], qv3yFullRXN[3];
    Double_t qv2xMiddTPC,    qv3xMiddTPC;
    Double_t qv2yMiddTPC,    qv3yMiddTPC;
    Double_t qv1xEastZDC,    qv1xWestZDC, qv1xFullZDC;
    Double_t qv1yEastZDC,    qv1yWestZDC, qv1yFullZDC;
    for (int iGap=0; iGap<4; iGap++)
    {
        qv2xEastTPC[iGap] = 0.; qv2yEastTPC[iGap] = 0.;
        qv3xEastTPC[iGap] = 0.; qv3yEastTPC[iGap] = 0.;
        qv2xWestTPC[iGap] = 0.; qv2yWestTPC[iGap] = 0.;
        qv3xWestTPC[iGap] = 0.; qv3yWestTPC[iGap] = 0.;
    }
    qv2xMiddTPC = 0.; qv2yMiddTPC = 0.;
    qv1xEastZDC = 0.; qv1yEastZDC = 0.;
    qv1xWestZDC = 0.; qv1yWestZDC = 0.;
    qv1xFullZDC = 0.; qv1yFullZDC = 0.;
    multEastZDC = 0;  multWestZDC = 0;
    for (int irxn=0; irxn<3; irxn++)
    {
        qv2xEastRXN[irxn] = 0.; qv2yEastRXN[irxn] = 0.;
        qv3xEastRXN[irxn] = 0.; qv3yEastRXN[irxn] = 0.;
        qv2xFullRXN[irxn] = 0.; qv2yFullRXN[irxn] = 0.;
        qv3xFullRXN[irxn] = 0.; qv3yFullRXN[irxn] = 0.;
        qv2xWestRXN[irxn] = 0.; qv2yWestRXN[irxn] = 0.;
        qv3xWestRXN[irxn] = 0.; qv3yWestRXN[irxn] = 0.;

        multEastRXN[irxn] = 0;  multWestRXN[irxn] = 0;
    }

    for (int iTr=0; iTr<nh; iTr++)
    {
        // if (charge[iTr] == 0) continue;
        vMom->SetPxPyPzE(momx[iTr],momy[iTr],momz[iTr],ene[iTr]);
        if (fRotate)
        {
          vMom->RotateZ(fPsiRP);
        }
        // if (vMom->Perp() < 0.1) continue;

        weight = 1.; //vMom->Perp();

        if (TMath::Abs(vMom->PseudoRapidity()) < 0.35)
        {
            qv2xMiddTPC += weight * TMath::Cos( 2. * vMom->Phi() );
            qv2yMiddTPC += weight * TMath::Sin( 2. * vMom->Phi() );
            qv3xMiddTPC += weight * TMath::Cos( 3. * vMom->Phi() );
            qv3yMiddTPC += weight * TMath::Sin( 3. * vMom->Phi() );
        }

        if (TMath::Abs(vMom->PseudoRapidity()) > 2. &&
          TMath::Abs(vMom->PseudoRapidity()) < 5.)
        {
          qv1xFullZDC += weight * TMath::Cos( 1. * vMom->Phi() );
          qv1yFullZDC += weight * TMath::Sin( 1. * vMom->Phi() );
          if (vMom->PseudoRapidity() < 0)
          {
            qv1xEastZDC += weight * TMath::Cos( 1. * vMom->Phi() );
            qv1yEastZDC += weight * TMath::Sin( 1. * vMom->Phi() );
            multEastZDC++;
          }
          if (vMom->PseudoRapidity() > 0)
          {
            qv1xWestZDC += weight * TMath::Cos( 1. * vMom->Phi() );
            qv1yWestZDC += weight * TMath::Sin( 1. * vMom->Phi() );
            multWestZDC++;
          }
        }

        for (int iGap=0; iGap<4; iGap++)
        {
            if (TMath::Abs(vMom->PseudoRapidity()) >= eta_gap.at(iGap) &&
                TMath::Abs(vMom->PseudoRapidity()) < 0.8 &&
                vMom->Perp() < 2. && vMom->PseudoRapidity() < 0)
            {
                qv2xEastTPC[iGap] += weight * TMath::Cos( 2. * vMom->Phi() );
                qv2yEastTPC[iGap] += weight * TMath::Sin( 2. * vMom->Phi() );
                qv3xEastTPC[iGap] += weight * TMath::Cos( 3. * vMom->Phi() );
                qv3yEastTPC[iGap] += weight * TMath::Sin( 3. * vMom->Phi() );
            }
            if (TMath::Abs(vMom->PseudoRapidity()) >= eta_gap.at(iGap) &&
                TMath::Abs(vMom->PseudoRapidity()) < 0.8 &&
                vMom->Perp() < 2. && vMom->PseudoRapidity() > 0)
            {
                qv2xWestTPC[iGap] += weight * TMath::Cos( 2. * vMom->Phi() );
                qv2yWestTPC[iGap] += weight * TMath::Sin( 2. * vMom->Phi() );
                qv3xWestTPC[iGap] += weight * TMath::Cos( 3. * vMom->Phi() );
                qv3yWestTPC[iGap] += weight * TMath::Sin( 3. * vMom->Phi() );
            }
        }
        for (int irxn=0; irxn<3; irxn++)
        {
            if (TMath::Abs(vMom->PseudoRapidity()) >= eta_rxn.at(irxn).first &&
                TMath::Abs(vMom->PseudoRapidity()) < eta_rxn.at(irxn).second &&
                momz[iTr] < 0)
            {
                qv2xEastRXN[irxn] += weight * TMath::Cos( 2. * vMom->Phi() );
                qv2yEastRXN[irxn] += weight * TMath::Sin( 2. * vMom->Phi() );
                qv3xEastRXN[irxn] += weight * TMath::Cos( 3. * vMom->Phi() );
                qv3yEastRXN[irxn] += weight * TMath::Sin( 3. * vMom->Phi() );
                multEastRXN[irxn]++;
            }
            if (TMath::Abs(vMom->PseudoRapidity()) >= eta_rxn.at(irxn).first &&
                TMath::Abs(vMom->PseudoRapidity()) < eta_rxn.at(irxn).second &&
                momz[iTr] > 0)
            {
                qv2xWestRXN[irxn] += weight * TMath::Cos( 2. * vMom->Phi() );
                qv2yWestRXN[irxn] += weight * TMath::Sin( 2. * vMom->Phi() );
                qv3xWestRXN[irxn] += weight * TMath::Cos( 3. * vMom->Phi() );
                qv3yWestRXN[irxn] += weight * TMath::Sin( 3. * vMom->Phi() );
                multWestRXN[irxn]++;
            }
        }
    }

    fQv2MiddTPC.Set(qv2xMiddTPC,qv2yMiddTPC);
    fQv3MiddTPC.Set(qv3xMiddTPC,qv3yMiddTPC);
    fQv1EastZDC.Set(qv1xEastZDC,qv1yEastZDC);
    fQv1WestZDC.Set(qv1xWestZDC,qv1yWestZDC);
    if (multEastZDC >=2 && multWestZDC >= 2)
    {
      qv1xFullZDC = qv1xEastZDC + qv1xWestZDC;
      qv1yFullZDC = qv1yEastZDC + qv1yWestZDC;
    }
    else
    {
      qv1xFullZDC = -999.;
      qv1yFullZDC = -999.;
    }
    fQv1FullZDC.Set(qv1xFullZDC,qv1yFullZDC);

    for (int iGap=0; iGap<4; iGap++)
    {
        fQv2EastTPC[iGap].Set(qv2xEastTPC[iGap],qv2yEastTPC[iGap]);
        fQv3EastTPC[iGap].Set(qv3xEastTPC[iGap],qv3yEastTPC[iGap]);
        fQv2WestTPC[iGap].Set(qv2xWestTPC[iGap],qv2yWestTPC[iGap]);
        fQv3WestTPC[iGap].Set(qv3xWestTPC[iGap],qv3yWestTPC[iGap]);
    }
    for (int irxn=0; irxn<3; irxn++)
    {
        if (multEastRXN[irxn] >= 2 && multWestRXN[irxn] >= 2)
        {
            qv2xFullRXN[irxn] = qv2xEastRXN[irxn] + qv2xWestRXN[irxn];
            qv2yFullRXN[irxn] = qv2yEastRXN[irxn] + qv2yWestRXN[irxn];
            qv3xFullRXN[irxn] = qv3xEastRXN[irxn] + qv3xWestRXN[irxn];
            qv3yFullRXN[irxn] = qv3yEastRXN[irxn] + qv3yWestRXN[irxn];
        }
        else
        {
            qv2xFullRXN[irxn] = -999.;
            qv2yFullRXN[irxn] = -999.;
            qv3xFullRXN[irxn] = -999.;
            qv3yFullRXN[irxn] = -999.;
        }
        
        fQv2EastRXN[irxn].Set(qv2xEastRXN[irxn],qv2yEastRXN[irxn]);
        fQv3EastRXN[irxn].Set(qv3xEastRXN[irxn],qv3yEastRXN[irxn]);
        fQv2WestRXN[irxn].Set(qv2xWestRXN[irxn],qv2yWestRXN[irxn]);
        fQv3WestRXN[irxn].Set(qv3xWestRXN[irxn],qv3yWestRXN[irxn]);
        fQv2FullRXN[irxn].Set(qv2xFullRXN[irxn],qv2yFullRXN[irxn]);
        fQv3FullRXN[irxn].Set(qv3xFullRXN[irxn],qv3yFullRXN[irxn]);
    }
}

Double_t FlowCalculator::GetPsiEP(const TVector2 &qv, Double_t harm)
{
    return (Double_t) (1/harm) * TMath::ATan2(qv.Y(),qv.X());
    //return (Double_t) TMath::ATan2(qv.Y(),qv.X());
}

#endif