my_star/StRoot/StFlowAnalysisMaker/StFlowAnalysisMaker.cxx

////////////////////////////////////////////////////////////////////////////
//
// $Id: StFlowAnalysisMaker.cxx,v 1.103 2011/07/25 15:54:42 posk Exp $
//
// Authors: Raimond Snellings and Art Poskanzer, LBNL, Aug 1999
//          FTPC added by Markus Oldenburg, MPI, Dec 2000
//
////////////////////////////////////////////////////////////////////////////
//
// Description:  Maker to analyze Flow using StFlowEvent
//
////////////////////////////////////////////////////////////////////////////

#include <Stiostream.h>
#include <stdlib.h>
#include <math.h>
#include "StMaker.h"
#include "StFlowAnalysisMaker.h"
#include "StFlowMaker/StFlowMaker.h"
#include "StFlowMaker/StFlowEvent.h"
#include "StFlowMaker/StFlowConstants.h"
#include "StFlowMaker/StFlowSelection.h"
#include "StFlowMaker/StFlowCutEvent.h"
#include "StEnumerations.h"
#include "PhysicalConstants.h"
#include "SystemOfUnits.h"
#include "TVector2.h"
#include "TFile.h"
#include "TString.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TOrdCollection.h"
#include "StMessMgr.h"
#include "TMath.h"
#include "TText.h"
#include "TF1.h"
#define PR(x) cout << "##### FlowAnalysis: " << (#x) << " = " << (x) << endl;

ClassImp(StFlowAnalysisMaker)

//-----------------------------------------------------------------------

StFlowAnalysisMaker::StFlowAnalysisMaker(const Char_t* name): StMaker(name),
  MakerName(name) {
  pFlowSelect = new StFlowSelection();
  SetHistoRanges();
  SetPtRange_for_vEta(0., 0.);
  SetEtaRange_for_vPt(0., 0.);
}

StFlowAnalysisMaker::StFlowAnalysisMaker(const Char_t* name,
					 const StFlowSelection& flowSelect) :
  StMaker(name), MakerName(name) {
  pFlowSelect = new StFlowSelection(flowSelect); //copy constructor
  SetHistoRanges();
  SetPtRange_for_vEta(0., 0.);
  SetEtaRange_for_vPt(0., 0.);
}

//-----------------------------------------------------------------------

StFlowAnalysisMaker::~StFlowAnalysisMaker() {
}

//-----------------------------------------------------------------------

Int_t StFlowAnalysisMaker::Make() {
  // Make histograms

  // Get another pointer to StFlowEvent
  pFlowMaker = NULL;
  pFlowMaker = (StFlowMaker*)GetMaker("Flow");
  if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
  if (pFlowEvent && pFlowSelect->Select(pFlowEvent)) {     // event selected
    if (FillFromFlowEvent()) {               // get event quantities
      PsiShift();							 // Take care of Psi shifting - //# Added by jcampbell
      FillEventHistograms();                 // fill from FlowEvent
      FillParticleHistograms();              // fill particle histograms
    } else {
      gMessMgr->Info("##### FlowAnalysis: Event psi = 0");
    }
  } else {
    gMessMgr->Info("##### FlowAnalysis: FlowEvent pointer null");
    return kStOK;
  }// selected events
    
  if (Debug()) StMaker::PrintInfo();

  return kStOK;
}

//-----------------------------------------------------------------------

Int_t StFlowAnalysisMaker::Init() {
  // Book histograms

  StFlowMaker* pFlowMaker = NULL;
  pFlowMaker = (StFlowMaker*)GetMaker("Flow");
  Bool_t reCentCalc = pFlowMaker->ReCentCalc();

  float ptMaxPart = Flow::ptMaxPart;
  if (pFlowSelect->PtMaxPart()) {
    ptMaxPart = pFlowSelect->PtMaxPart();
  }
  int nPtBinsPart = Flow::nPtBinsPart;
  if (pFlowSelect->PtBinsPart()) {
    nPtBinsPart = pFlowSelect->PtBinsPart();
  }
  xLabel = "Pseudorapidity";
  if (strlen(pFlowSelect->PidPart()) != 0) { xLabel = "Rapidity"; }

  const float triggerMin      =  -0.5;
  const float triggerMax      =  10.5;
  const float chargeMin       =  -2.5;
  const float chargeMax       =   2.5; 
  const float dcaMin          =    0.;
  const float dcaMax          =   0.3;
  const float glDcaMax        =   3.6; 
  const float chi2Min         =    0.;
  const float chi2Max         =    5.; 
  const float fitPtsMinTpc    =  -0.5;
  const float fitPtsMaxTpc    =  60.5; 
  const float maxPtsMinTpc    =  -0.5;
  const float maxPtsMaxTpc    =  60.5; 
  const float fitPtsMinFtpc   =  -0.5;
  const float fitPtsMaxFtpc   =  12.5; 
  const float maxPtsMinFtpc   =  -0.5;
  const float maxPtsMaxFtpc   =  12.5; 
  const float fitOverMaxMin   =    0.;
  const float fitOverMaxMax   =   1.2; 
  const float origMultMin     =    0.;
  const float origMultMax     = 3000.; 
  const float MultEtaMin      =    0.;
  const float MultEtaMax      = 1000.; 
  const float totalMultMin    =    0.;
  const float totalMultMax    = 2000.; 
  const float corrMultMin     =    0.;
  const float corrMultMax     = 2000.; 
  const float multOverOrigMin =    0.;
  const float multOverOrigMax =    1.; 
  const float vertexZMin      =-100.5;
  const float vertexZMax      = 100.5; 
  const float vertexXYMin     =   -1.;
  const float vertexXYMax     =    1.; 
  const float QXYMin          =  -0.5;
  const float QXYMax          =   0.5; 
  const float etaSymZMin      = -1.15; 
  const float etaSymZMax      =  1.15; 
  const float etaSymMin       =   -6.; 
  const float etaSymMax       =    6.; 
  const float phiMin          =    0.;
  const float phiMax          = twopi; 
  const float psiMin          =    0.;
  const float psiMax          = twopi; 
  const float multMin         =    0.;
  const float multMax         = 2000.;
  const float qMin            =    0.;
  const float pidMin          =  -10.;
  const float pidMax          =   10.;
  const float centMin         =  -0.5;
  const float centMax         =   9.5;
  const float pMin            =  -2.5;
  const float pMax            =   1.5;
  const float dEdxMax       = 0.00004;
  const float qMax             =  3.5;

  enum { nTriggerBins      = 11,
	 nChargeBins       = 50,
	 nDcaBins          = 60,
	 nChi2Bins         = 50,
	 nFitPtsBinsTpc    = 61,
	 nFitPtsBinsFtpc   = 13,
	 nMaxPtsBinsTpc    = 61,
	 nMaxPtsBinsFtpc   = 13,
	 nFitOverMaxBins   = 40,
	 nOrigMultBins     = 60,
	 nMultEtaBins      = 50,
	 nTotalMultBins    = 40,
	 nMultOverOrigBins = 50,
	 nMultPartBins     = 40,
	 nVertexZBins      = 51,
	 nVertexXYBins     = 50,
	 nQXYBins          = 50,
	 nEtaSymBins       = 45,
	 nPhi3DBins        = 18,
	 nPsiBins          = 36,
	 nMultBins         = 40,
	 nPidBins          = 50,
         nCentBins         = 10,
	 nDedxBins        = 200,
	 nMomenBins       = 200,
	 n_qBins          =  50
  };

  // Trigger
  mHistTrigger = new TH1F("Flow_Trigger", "Flow_Trigger",
      nTriggerBins, triggerMin, triggerMax);
  mHistTrigger->SetXTitle("Trig: 0 mb+cen, 1 mb, 2 central, 3 laser, 10 other");
  mHistTrigger->SetYTitle("Counts");

  // Charge
  // Ftpc
  mHistChargeFtpc = new TH1F("Flow_Charge_Ftpc", "Flow_Charge_Ftpc",
      nChargeBins, chargeMin, chargeMax);
  mHistChargeFtpc->SetXTitle("Charge");
  mHistChargeFtpc->SetYTitle("Counts");
    
  // Distance of closest approach
  // Tpc
  mHistDcaTpc = new TH1F("Flow_Dca_Tpc", "Flow_Dca_Tpc",
      nDcaBins, dcaMin, dcaMax);
  mHistDcaTpc->SetXTitle("Track dca to Vertex (cm)");
  mHistDcaTpc->SetYTitle("Counts");

  // Ftpc
  mHistDcaFtpc = new TH1F("Flow_Dca_Ftpc", "Flow_Dca_Ftpc",
      nDcaBins, dcaMin, dcaMax);
  mHistDcaFtpc->SetXTitle("Track dca to Vertex (cm)");
  mHistDcaFtpc->SetYTitle("Counts");
    

  // Distance of closest approach for global tracks
  // Tpc
  mHistDcaGlobalTpc = new TH1F("Flow_DcaGlobal_Tpc", "Flow_DcaGlobal_Tpc",
      nDcaBins, dcaMin, glDcaMax);
  mHistDcaGlobalTpc->SetXTitle("Global Track dca (cm)");
  mHistDcaGlobalTpc->SetYTitle("Counts");

  // Ftpc
  mHistDcaGlobalFtpc = new TH1F("Flow_DcaGlobal_Ftpc", "Flow_DcaGlobal_Ftpc",
      nDcaBins, dcaMin, glDcaMax);
  mHistDcaGlobalFtpc->SetXTitle("Global Track dca (cm)");
  mHistDcaGlobalFtpc->SetYTitle("Counts");

    
  // Chi2
  // Tpc
  mHistChi2Tpc = new TH1F("Flow_Chi2_Tpc", "Flow_Chi2_Tpc",
      nChi2Bins, chi2Min, chi2Max);
  mHistChi2Tpc->SetXTitle("Chi square per df");
  mHistChi2Tpc->SetYTitle("Counts");
  
  // Ftpc
  mHistChi2Ftpc = new TH1F("Flow_Chi2_Ftpc", "Flow_Chi2_Ftpc",
      nChi2Bins, chi2Min, chi2Max);
  mHistChi2Ftpc->SetXTitle("Chi square per df");
  mHistChi2Ftpc->SetYTitle("Counts");
  
  
  // FitPts
  // Tpc
  mHistFitPtsTpc = new TH1F("Flow_FitPts_Tpc", "Flow_FitPts_Tpc",
      nFitPtsBinsTpc, fitPtsMinTpc, fitPtsMaxTpc);
  mHistFitPtsTpc->SetXTitle("Fit Points");
  mHistFitPtsTpc->SetYTitle("Counts");

  // Ftpc
  mHistFitPtsFtpc = new TH1F("Flow_FitPts_Ftpc", "Flow_FitPts_Ftpc",
      nFitPtsBinsFtpc, fitPtsMinFtpc, fitPtsMaxFtpc);
  mHistFitPtsFtpc->SetXTitle("Fit Points");
  mHistFitPtsFtpc->SetYTitle("Counts");
    

  // MaxPts
  // Tpc
  mHistMaxPtsTpc = new TH1F("Flow_MaxPts_Tpc ", "Flow_MaxPts_Tpc ",
      nMaxPtsBinsTpc , maxPtsMinTpc , maxPtsMaxTpc );
  mHistMaxPtsTpc ->SetXTitle("Max Points");
  mHistMaxPtsTpc ->SetYTitle("Counts");
  
  // Ftpc
  mHistMaxPtsFtpc = new TH1F("Flow_MaxPts_Ftpc", "Flow_MaxPts_Ftpc",
      nMaxPtsBinsFtpc, maxPtsMinFtpc, maxPtsMaxFtpc);
  mHistMaxPtsFtpc->SetXTitle("Max Points");
  mHistMaxPtsFtpc->SetYTitle("Counts");
    

  // FitOverMax
  // Tpc
  mHistFitOverMaxTpc = new TH1F("Flow_FitOverMax_Tpc", "Flow_FitOverMax_Tpc",
      nFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
  mHistFitOverMaxTpc->SetXTitle("Fit Points / Max Points");
  mHistFitOverMaxTpc->SetYTitle("Counts");
    
  // Ftpc
  mHistFitOverMaxFtpc = new TH1F("Flow_FitOverMax_Ftpc", "Flow_FitOverMax_Ftpc",
      nFitOverMaxBins, fitOverMaxMin, fitOverMaxMax);
  mHistFitOverMaxFtpc->SetXTitle("Fit Points / Max Points");
  mHistFitOverMaxFtpc->SetYTitle("Counts");
    

  // OrigMult
  mHistOrigMult = new TH1F("Flow_OrigMult", "Flow_OrigMult",
      nOrigMultBins, origMultMin, origMultMax);
  mHistOrigMult->SetXTitle("Original Mult");
  mHistOrigMult->SetYTitle("Counts");

  // MultEta
  mHistMultEta = new TH1F("Flow_MultEta", "Flow_MultEta",
      nMultEtaBins, MultEtaMin, MultEtaMax);
  mHistMultEta->SetXTitle("Mult for Centrality");
  mHistMultEta->SetYTitle("Counts");
    
  // Mult
  mHistMult = new TH1F("Flow_Mult", "Flow_Mult",
      nTotalMultBins, totalMultMin, totalMultMax);
  mHistMult->SetXTitle("Mult");
  mHistMult->SetYTitle("Counts");
    
  // MultOverOrig
  mHistMultOverOrig = new TH1F("Flow_MultOverOrig", "Flow_MultOverOrig",
      nMultOverOrigBins, multOverOrigMin, multOverOrigMax);
  mHistMultOverOrig->SetXTitle("Mult / Orig. Mult");
  mHistMultOverOrig->SetYTitle("Counts");
    
  // Mult correlated with the event planes
  mHistMultPart = new TH1F("Flow_MultPart", "Flow_MultPart",
      nMultPartBins, corrMultMin, corrMultMax);
  mHistMultPart->SetXTitle("Mult of Correlated Particles");
  mHistMultPart->SetYTitle("Counts");
    
  // VertexZ
  mHistVertexZ = new TH1F("Flow_VertexZ", "Flow_VertexZ",
      nVertexZBins, vertexZMin, vertexZMax);
  mHistVertexZ->SetXTitle("Vertex Z (cm)");
  mHistVertexZ->SetYTitle("Counts");
    
  // VertexXY
  mHistVertexXY2D = new TH2F("Flow_VertexXY2D", "Flow_VertexXY2D",
			     nVertexXYBins, vertexXYMin, vertexXYMax,
			     nVertexXYBins, vertexXYMin, vertexXYMax);
  mHistVertexXY2D->SetXTitle("Vertex X (cm)");
  mHistVertexXY2D->SetYTitle("Vertex Y (cm)");
    
  // EtaSym vs. Vertex Z Tpc
  mHistEtaSymVerZ2DTpc = new TH2F("Flow_EtaSymVerZ2D_Tpc", "Flow_EtaSymVerZ2D_Tpc",
    nVertexZBins, vertexZMin, vertexZMax, nEtaSymBins, etaSymZMin, etaSymZMax);
  mHistEtaSymVerZ2DTpc->SetXTitle("Vertex Z (cm)");
  mHistEtaSymVerZ2DTpc->SetYTitle("Eta Symmetry");

  // EtaSym vs. Vertex Z Ftpc
  mHistEtaSymVerZ2DFtpc = new TH2F("Flow_EtaSymVerZ2D_Ftpc", "Flow_EtaSymVerZ2D_Ftpc",
    nVertexZBins, vertexZMin, vertexZMax, nEtaSymBins, etaSymZMin, etaSymZMax);
  mHistEtaSymVerZ2DFtpc->SetXTitle("Vertex Z (cm)");
  mHistEtaSymVerZ2DFtpc->SetYTitle("Eta Symmetry");

  // EtaSym Tpc
  mHistEtaSymTpc = new TH1F("Flow_EtaSym_Tpc", "Flow_EtaSym_Tpc",
      nEtaSymBins, etaSymMin, etaSymMax);
  mHistEtaSymTpc->SetXTitle("Eta Symmetry Ratio TPC");
  mHistEtaSymTpc->SetYTitle("Counts");
    
  // EtaSym Ftpc
  mHistEtaSymFtpc = new TH1F("Flow_EtaSym_Ftpc", "Flow_EtaSym_Ftpc",
      nEtaSymBins, etaSymMin, etaSymMax);
  mHistEtaSymFtpc->SetXTitle("Eta Symmetry Ratio FTPC");
  mHistEtaSymFtpc->SetYTitle("Counts");
    
  // EtaPtPhi
  mHistEtaPtPhi3D = new TH3F("Flow_EtaPtPhi3D", "Flow_EtaPtPhi3D",
      mNEtaBins, mEtaMin, mEtaMax, Flow::nPtBins, Flow::ptMin, 
			     Flow::ptMax, nPhi3DBins, phiMin, phiMax);
  mHistEtaPtPhi3D->SetXTitle("Eta");
  mHistEtaPtPhi3D->SetYTitle("Pt (GeV/c)");
  mHistEtaPtPhi3D->SetZTitle("Phi (rad)");
    
  // Yield for all particles
  mHistYieldAll2D = new TH2D("Flow_YieldAll2D", "Flow_YieldAll2D",
    mNEtaBins, mEtaMin, mEtaMax, Flow::nPtBins, Flow::ptMin, 
			     Flow::ptMax);
  mHistYieldAll2D->Sumw2();
  mHistYieldAll2D->SetXTitle("Pseudorapidty");
  mHistYieldAll2D->SetYTitle("Pt (GeV/c)");

  // Yield for particles correlated with the event plane
  mHistYieldPart2D = new TH2D("Flow_YieldPart2D", "Flow_YieldPart2D",
    mNEtaBins, mEtaMin, mEtaMax, nPtBinsPart, Flow::ptMin, 
			      ptMaxPart);
  mHistYieldPart2D->Sumw2();
  mHistYieldPart2D->SetXTitle((char*)xLabel.Data());
  mHistYieldPart2D->SetYTitle("Pt (GeV/c)");

  // Mean Eta in each bin
  mHistBinEta = new TProfile("Flow_Bin_Eta", "Flow_Bin_Eta",
    mNEtaBins, mEtaMin, mEtaMax, mEtaMin, mEtaMax, "");
  mHistBinEta->SetXTitle((char*)xLabel.Data());
  mHistBinEta->SetYTitle("<Eta>");
  
  // Mean Pt in each bin
  mHistBinPt = new TProfile("Flow_Bin_Pt", "Flow_Bin_Pt",
    nPtBinsPart, Flow::ptMin, ptMaxPart, Flow::ptMin, ptMaxPart, "");
  mHistBinPt->SetXTitle("Pt (GeV/c)");
  mHistBinPt->SetYTitle("<Pt> (GeV/c)");
      
  // PID pi+
  mHistPidPiPlus = new TH1F("Flow_PidPiPlus", "Flow_PidPiPlus",
			    nPidBins, pidMin, pidMax);
  mHistPidPiPlus->SetXTitle("(PID - Mean) / Resolution");
  mHistPidPiPlus->SetYTitle("Counts");
    
  // PID pi-
  mHistPidPiMinus = new TH1F("Flow_PidPiMinus", "Flow_PidPiMinus",
			     nPidBins, pidMin, pidMax);
  mHistPidPiMinus->SetXTitle("(PID - Mean) / Resolution");
  mHistPidPiMinus->SetYTitle("Counts");
    
  // PID proton
  mHistPidProton = new TH1F("Flow_PidProton", "Flow_PidProton",
			    nPidBins, pidMin, pidMax);
  mHistPidProton->SetXTitle("(PID - Mean) / Resolution");
  mHistPidProton->SetYTitle("Counts");

  // PID anti proton
  mHistPidAntiProton = new TH1F("Flow_PidAntiProton", "Flow_PidAntiProton",
				nPidBins, pidMin, pidMax);
  mHistPidAntiProton->SetXTitle("(PID - Mean) / Resolution");
  mHistPidAntiProton->SetYTitle("Counts");

  // PID Kplus
  mHistPidKplus = new TH1F("Flow_PidKplus", "Flow_PidKplus",
			   nPidBins, pidMin, pidMax);
  mHistPidKplus->SetXTitle("(PID - Mean) / Resolution");
  mHistPidKplus->SetYTitle("Counts");

  // PID Kminus
  mHistPidKminus = new TH1F("Flow_PidKminus", "Flow_PidKminus",
			    nPidBins, pidMin, pidMax);
  mHistPidKminus->SetXTitle("(PID - Mean) / Resolution");
  mHistPidKminus->SetYTitle("Counts");

  // PID deuteron
  mHistPidDeuteron = new TH1F("Flow_PidDeuteron", "Flow_PidDeuteron",
			      nPidBins, pidMin, pidMax);
  mHistPidDeuteron->SetXTitle("(PID - Mean) / Resolution");
  mHistPidDeuteron->SetYTitle("Counts");

  // PID anti deuteron
  mHistPidAntiDeuteron = new TH1F("Flow_PidAntiDeuteron", 
				  "Flow_PidAntiDeuteron",
				  nPidBins, pidMin, pidMax);
  mHistPidAntiDeuteron->SetXTitle("(PID - Mean) / Resolution");
  mHistPidAntiDeuteron->SetYTitle("Counts");

  // PID electron
  mHistPidElectron = new TH1F("Flow_PidElectron", "Flow_PidElectron",
			      nPidBins, pidMin, pidMax);
  mHistPidElectron->SetXTitle("(PID - Mean) / Resolution");
  mHistPidElectron->SetYTitle("Counts");

  // PID positron
  mHistPidPositron = new TH1F("Flow_PidPositron", "Flow_PidPositron",
			      nPidBins, pidMin, pidMax);
  mHistPidPositron->SetXTitle("(PID - Mean) / Resolution");
  mHistPidPositron->SetYTitle("Counts");

  // PID pi+ selected
  mHistPidPiPlusPart = new TH1F("Flow_PidPiPlusPart", 
			       "Flow_PidPiPlusPart",
			       nPidBins, pidMin, pidMax);
  mHistPidPiPlusPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidPiPlusPart->SetYTitle("Counts");
    
  // PID pi- selected
  mHistPidPiMinusPart = new TH1F("Flow_PidPiMinusPart", 
				"Flow_PidPiMinusPart",
				nPidBins, pidMin, pidMax);
  mHistPidPiMinusPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidPiMinusPart->SetYTitle("Counts");
    
  // PID proton selected
  mHistPidProtonPart = new TH1F("Flow_PidProtonPart", 
			       "Flow_PidProtonPart",
			       nPidBins, pidMin, pidMax);
  mHistPidProtonPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidProtonPart->SetYTitle("Counts");

  // PID anti proton selected
  mHistPidAntiProtonPart = new TH1F("Flow_PidAntiProtonPart", 
				   "Flow_PidAntiProtonPart",
				   nPidBins, pidMin, pidMax);
  mHistPidAntiProtonPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidAntiProtonPart->SetYTitle("Counts");

  // PID Kplus selected
  mHistPidKplusPart = new TH1F("Flow_PidKplusPart", 
			      "Flow_PidKplusPart",
			      nPidBins, pidMin, pidMax);
  mHistPidKplusPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidKplusPart->SetYTitle("Counts");

  // PID Kminus selected
  mHistPidKminusPart = new TH1F("Flow_PidKminusPart", 
			       "Flow_PidKminusPart",
			       nPidBins, pidMin, pidMax);
  mHistPidKminusPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidKminusPart->SetYTitle("Counts");

  // PID deuteron selected
  mHistPidDeuteronPart = new TH1F("Flow_PidDeuteronPart", 
				 "Flow_PidDeuteronPart",
				 nPidBins, pidMin, pidMax);
  mHistPidDeuteronPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidDeuteronPart->SetYTitle("Counts");

  // PID anti deuteron selected
  mHistPidAntiDeuteronPart = new TH1F("Flow_PidAntiDeuteronPart", 
				     "Flow_PidAntiDeuteronPart",
				     nPidBins, pidMin, pidMax);
  mHistPidAntiDeuteronPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidAntiDeuteronPart->SetYTitle("Counts");

  // PID electron selected
  mHistPidElectronPart = new TH1F("Flow_PidElectronPart", 
				 "Flow_PidElectronPart",
				 nPidBins, pidMin, pidMax);
  mHistPidElectronPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidElectronPart->SetYTitle("Counts");

  // PID positron selected
  mHistPidPositronPart = new TH1F("Flow_PidPositronPart", 
				 "Flow_PidPositronPart",
				 nPidBins, pidMin, pidMax);
  mHistPidPositronPart->SetXTitle("(PID - Mean) / Resolution");
  mHistPidPositronPart->SetYTitle("Counts");

  // PID multiplicities selected
  mHistPidMult = new TProfile("Flow_PidMult", "Flow_PidMult",
			      13, 0.5, 13.5, 0., 10000., "");
  mHistPidMult->SetXTitle("all, h+, h-, pi+, pi-, pr+, pr-, K+, K-, d+, d-, e-, e+");
  mHistPidMult->SetYTitle("Multiplicity");
    
  // Centrality
  mHistCent = new TH1F("Flow_Cent", "Flow_Cent",
		       nCentBins, centMin, centMax);
  mHistCent->SetXTitle("Centrality Bin");
  mHistCent->SetYTitle("Counts");
    
  // CTB versus ZDC
  mHistCTBvsZDC2D = new TH2F("Flow_CTBvsZDC2D", "Flow_CTBvsZDC2D",
			       125, 0, 500,
			       125, 0, 40000);
  mHistCTBvsZDC2D->SetXTitle("ZDC sum");
  mHistCTBvsZDC2D->SetYTitle("CTB sum");

  // MeanDedxPos
  mHistMeanDedxPos2D = new TH2F("Flow_MeanDedxPos2D",
				"Flow_MeanDedxPos2D",
				nMomenBins, pMin, pMax,
				nDedxBins, 0, dEdxMax);
  mHistMeanDedxPos2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxPos2D->SetYTitle("mean dEdx");
  
  // MeanDedxNeg
  mHistMeanDedxNeg2D = new TH2F("Flow_MeanDedxNeg2D",
				"Flow_MeanDedxNeg2D",
				nMomenBins, pMin, pMax,
				nDedxBins, 0, dEdxMax);
  mHistMeanDedxNeg2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxNeg2D->SetYTitle("mean dEdx");
  
  // MeanDedx PiPlus
  mHistMeanDedxPiPlus2D = new TH2F("Flow_MeanDedxPiPlus2D", 
				   "Flow_MeanDedxPiPlus2D",
				   nMomenBins, pMin, pMax,
				   nDedxBins, 0, dEdxMax);
  mHistMeanDedxPiPlus2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxPiPlus2D->SetYTitle("mean dEdx");
  
  // MeanDedxPiMinus
  mHistMeanDedxPiMinus2D = new TH2F("Flow_MeanDedxPiMinus2D", 
				    "Flow_MeanDedxPiMinus2D",
				    nMomenBins, pMin, pMax,
				    nDedxBins, 0, dEdxMax);
  mHistMeanDedxPiMinus2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxPiMinus2D->SetYTitle("mean dEdx");
  
  // MeanDedxProton
  mHistMeanDedxProton2D = new TH2F("Flow_MeanDedxProton2D", 
				   "Flow_MeanDedxProton2D",
				   nMomenBins, pMin, pMax,
				   nDedxBins, 0, dEdxMax);
  mHistMeanDedxProton2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxProton2D->SetYTitle("mean dEdx");
  
  // MeanDedxPbar
  mHistMeanDedxPbar2D = new TH2F("Flow_MeanDedxPbar2D", 
				 "Flow_MeanDedxPbar2D",
				 nMomenBins, pMin, pMax,
				 nDedxBins, 0, dEdxMax);
  mHistMeanDedxPbar2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxPbar2D->SetYTitle("mean dEdx");
  
  // MeanDedxKplus
  mHistMeanDedxKplus2D = new TH2F("Flow_MeanDedxKplus2D", 
				  "Flow_MeanDedxKplus2D",
				  nMomenBins, pMin, pMax,
				  nDedxBins, 0, dEdxMax);
  mHistMeanDedxKplus2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxKplus2D->SetYTitle("mean dEdx");
  
  // MeanDedxKminus
  mHistMeanDedxKminus2D = new TH2F("Flow_MeanDedxKminus2D", 
				   "Flow_MeanDedxKminus2D",
				   nMomenBins, pMin, pMax,
				   nDedxBins, 0, dEdxMax);
  mHistMeanDedxKminus2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxKminus2D->SetYTitle("mean dEdx");
  
  // MeanDedxDeuteron
  mHistMeanDedxDeuteron2D = new TH2F("Flow_MeanDedxDeuteron2D", 
				     "Flow_MeanDedxDeuteron2D",
				     nMomenBins, pMin, pMax,
				     nDedxBins, 0, dEdxMax);
  mHistMeanDedxDeuteron2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxDeuteron2D->SetYTitle("mean dEdx");

  // MeanDedxAntiDeuteron
  mHistMeanDedxAntiDeuteron2D = new TH2F("Flow_MeanDedxAntiDeuteron2D", 
					 "Flow_MeanDedxAntiDeuteron2D",
					 nMomenBins, pMin, pMax,
					 nDedxBins, 0, dEdxMax);
  mHistMeanDedxAntiDeuteron2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxAntiDeuteron2D->SetYTitle("mean dEdx");
  
  // MeanDedxElectron
  mHistMeanDedxElectron2D = new TH2F("Flow_MeanDedxElectron2D", 
				     "Flow_MeanDedxElectron2D",
				     nMomenBins, pMin, pMax,
				     nDedxBins, 0, dEdxMax);
  mHistMeanDedxElectron2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxElectron2D->SetYTitle("mean dEdx");
  
  // MeanDedxPositron
  mHistMeanDedxPositron2D = new TH2F("Flow_MeanDedxPositron2D", 
				     "Flow_MeanDedxPositron2D",
				     nMomenBins, pMin, pMax,
				     nDedxBins, 0, dEdxMax);
  mHistMeanDedxPositron2D->SetXTitle("log(momentum) (GeV/c)");
  mHistMeanDedxPositron2D->SetYTitle("mean dEdx");
  
  // ZDCSMD test
  mZDC_SMD_west_vert = new TH1F("Flow_ZDC_SMD_west_vert","Flow_ZDC_SMD_west_vert",7,0.5,7.5);
  mZDC_SMD_east_vert = new TH1F("Flow_ZDC_SMD_east_vert","Flow_ZDC_SMD_east_vert",7,0.5,7.5);
  mZDC_SMD_west_hori = new TH1F("Flow_ZDC_SMD_west_hori","Flow_ZDC_SMD_west_hori",8,0.5,8.5);
  mZDC_SMD_east_hori = new TH1F("Flow_ZDC_SMD_east_hori","Flow_ZDC_SMD_east_hori",8,0.5,8.5);
  mHistZDCSMDPsiWgtEast = new TH1D("Flow_ZDCSMDPsiWgtEast","Flow_ZDCSMDPsiWgtEast",
				   Flow::zdcsmd_nPsiBins,-twopi/2.,twopi/2.);
  mHistZDCSMDPsiWgtWest = new TH1D("Flow_ZDCSMDPsiWgtWest","Flow_ZDCSMDPsiWgtWest",
				   Flow::zdcsmd_nPsiBins,-twopi/2.,twopi/2.);
  mHistZDCSMDPsiWgtTest = new TH1D("Flow_ZDCSMDPsiWgtTest","Flow_ZDCSMDPsiWgtTest",
                                   Flow::zdcsmd_nPsiBins,0.,twopi);
  mHistZDCSMDPsiWgtFull = new TH1D("Flow_ZDCSMDPsiWgtFull","Flow_ZDCSMDPsiWgtFull",
                                   Flow::zdcsmd_nPsiBins,0.,twopi);
  mHistZDCSMDPsiCorTest = new TH1D("Flow_ZDCSMDPsiCorTest","Flow_ZDCSMDPsiCorTest",
                                   Flow::zdcsmd_nPsiBins,-twopi/2.,twopi/2.);
  mHistZDCSMDPsiCorFull = new TH1D("Flow_ZDCSMDPsiCorFull","Flow_ZDCSMDPsiWgtFull",
                                   Flow::zdcsmd_nPsiBins,-twopi/2.,twopi/2.);

  TString* histTitle;
  // ------------- Group II histograms ------------ //
  //# Added by jcampbell

  // One for each selection
  for (int k = 0; k < Flow::nSels; k++) {
    for (int j = 0; j < Flow::nHars; j++) {
        histTitle = new TString("Flow_PsiFull_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistPsiFull = new TH1D(histTitle->Data(), histTitle->Data(),100,0,twopi/2);
        histFull[k].histFullHar[j].mHistPsiFull->SetXTitle("Event Plane Angle (rad)");
        histFull[k].histFullHar[j].mHistPsiFull->SetYTitle("dN/d(EP)");
        delete histTitle;

        histTitle = new TString("Flow_PsiSubs_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistPsiSub2VsSub1 = new TH2D(histTitle->Data(), histTitle->Data(),100,0,twopi/2,100,0,twopi/2);
        histFull[k].histFullHar[j].mHistPsiSub2VsSub1->SetXTitle("Psi_a");
        histFull[k].histFullHar[j].mHistPsiSub2VsSub1->SetYTitle("Psi_b");
        delete histTitle;

        histTitle = new TString("Flow_RPMult_Sub1VsFull_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistRPMultSub1VsFull = new TH2D(histTitle->Data(), histTitle->Data(),1000,-0.5,999.5,500,-0.5,499.5);
        histFull[k].histFullHar[j].mHistRPMultSub1VsFull->SetXTitle("RP Mult - Full Event");
        histFull[k].histFullHar[j].mHistRPMultSub1VsFull->SetYTitle("RP Mult - Sub A");
        delete histTitle;

        histTitle = new TString("Flow_RPMult_Sub2VsFull_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistRPMultSub2VsFull = new TH2D(histTitle->Data(), histTitle->Data(),1000,-0.5,999.5,500,-0.5,499.5);
        histFull[k].histFullHar[j].mHistRPMultSub2VsFull->SetXTitle("RP Mult - Full Event");
        histFull[k].histFullHar[j].mHistRPMultSub2VsFull->SetYTitle("RP Mult - Sub B");
        delete histTitle;

        histTitle = new TString("Flow_RPMult_Sub2VsSub1_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistRPMultSub2VsSub1 = new TH2D(histTitle->Data(), histTitle->Data(),500,-0.5,499.5,500,-0.5,499.5);
        histFull[k].histFullHar[j].mHistRPMultSub2VsSub1->SetXTitle("RP Mult - Sub A");
        histFull[k].histFullHar[j].mHistRPMultSub2VsSub1->SetYTitle("RP Mult - Sub B");
        delete histTitle;

        histTitle = new TString("Flow_SubResSquared_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistResSquaredVsMult = new TProfile(histTitle->Data(), histTitle->Data(),1000,-0.5,999.5);
        histFull[k].histFullHar[j].mHistResSquaredVsMult->SetXTitle("Centrality Index");
        histFull[k].histFullHar[j].mHistResSquaredVsMult->SetYTitle("<Cos(2(Psi_a - Psi_b))>");
        delete histTitle;

        histTitle = new TString("Flow_dEdXVsRigidity_");
        *histTitle += k+1;
          histTitle->Append("_Har");
          *histTitle += j+1;
        histFull[k].histFullHar[j].mHistdEdXVsRigidity = new TH2D(histTitle->Data(), histTitle->Data(),200,-2,2,100,0,dEdxMax);
        histFull[k].histFullHar[j].mHistdEdXVsRigidity->SetXTitle("q*p");
        histFull[k].histFullHar[j].mHistdEdXVsRigidity->SetYTitle("dE/dx");
        delete histTitle;
    }
  }



  for (int i = 0; i < Flow::nSels * Flow::nSubs; i++) {

    // for sub-events
    for (int j = 0; j < Flow::nHars; j++) {
      float order = (float)(j + 1);

	  //# Added by jcampbell
     // Multiplicities 
      histTitle = new TString("Flow_Mult_Subs");
      *histTitle += i+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histSub[i].histSubHar[j].mHistMultSubs = new TH1F(histTitle->Data(),
        histTitle->Data(), nMultBins, multMin, multMax);
      histSub[i].histSubHar[j].mHistMultSubs->Sumw2();
      histSub[i].histSubHar[j].mHistMultSubs->SetXTitle
	("Mult");
      histSub[i].histSubHar[j].mHistMultSubs->SetYTitle("Counts");
      delete histTitle;

      // event planes
      histTitle = new TString("Flow_Psi_Subs");
      *histTitle += i+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histSub[i].histSubHar[j].mHistPsiSubs = new TH1F(histTitle->Data(),
	histTitle->Data(), nPsiBins, psiMin, (psiMax / order));
      histSub[i].histSubHar[j].mHistPsiSubs->Sumw2();
      histSub[i].histSubHar[j].mHistPsiSubs->SetXTitle
	("Event Plane Angle (rad)");
      histSub[i].histSubHar[j].mHistPsiSubs->SetYTitle("Counts");
      delete histTitle;
    }
  }

  for (int k = 0; k < Flow::nSels; k++) {

    // for each selection

    // cos(n*delta_Psi)
    histTitle = new TString("Flow_Cos_Sel");
    *histTitle += k+1;
    histFull[k].mHistCos = new TProfile(histTitle->Data(), histTitle->Data(),
      Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -1., 1., "");
    histFull[k].mHistCos->SetXTitle("Harmonic");
    histFull[k].mHistCos->SetYTitle("<cos(n*delta_Psi)>");
    delete histTitle;
    
    // resolution
    histTitle = new TString("Flow_Res_Sel");
    *histTitle += k+1;
    histFull[k].mHistRes = new TH1F(histTitle->Data(), histTitle->Data(),
      Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5);
    histFull[k].mHistRes->SetXTitle("Harmonic");
    histFull[k].mHistRes->SetYTitle("Resolution");
    delete histTitle;

    // vObs
    histTitle = new TString("Flow_vObs_Sel");
    *histTitle += k+1;
    histFull[k].mHist_vObs = new TProfile(histTitle->Data(), histTitle->Data(),
      Flow::nHars, 0.5, (float)(Flow::nHars) + 0.5, -1000., 1000., "");
    histFull[k].mHist_vObs->SetXTitle("Harmonic");
    histFull[k].mHist_vObs->SetYTitle("vObs (%)");
    delete histTitle;
    
    // for each harmonic
    for (int j = 0; j < Flow::nHars; j++) {
      float order  = (float)(j+1);

      // multiplicity
      histTitle = new TString("Flow_Mul_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistMult = new TH1F(histTitle->Data(),
        histTitle->Data(), nMultBins, multMin, multMax);
      histFull[k].histFullHar[j].mHistMult->SetXTitle("Multiplicity");
      histFull[k].histFullHar[j].mHistMult->SetYTitle("Counts");
      delete histTitle;

      // event plane
      histTitle = new TString("Flow_Psi_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistPsi = new TH1F(histTitle->Data(),
        histTitle->Data(), nPsiBins, psiMin, psiMax / order);
      histFull[k].histFullHar[j].mHistPsi->Sumw2();
      histFull[k].histFullHar[j].mHistPsi->SetXTitle
	("Event Plane Angle (rad)");
      histFull[k].histFullHar[j].mHistPsi->SetYTitle("Counts");
      delete histTitle;
      
      // phi lab
      histTitle = new TString("Flow_PhiLab_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistPhiLab = new TH1F(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histFullHar[j].mHistPhiLab->SetXTitle("Particle Lab Angle (rad)");
      histFull[k].histFullHar[j].mHistPhiLab->SetYTitle("Counts");
      histFull[k].histFullHar[j].mHistPhiLab->Sumw2(); // for scale
      delete histTitle;
      
      // Recenter
      histTitle = new TString("FlowReCentX_Sel");
      *histTitle += k+1;
      *histTitle += "_Har";
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistReCentX = new TProfile(histTitle->Data(),
        histTitle->Data(), 3, 0.5, 3.5);
      histFull[k].histFullHar[j].mHistReCentX->SetXTitle("FTPCE, FTPCW, TPCE, TPCW");
      histFull[k].histFullHar[j].mHistReCentX->SetYTitle("<cos n #phi>");
      delete histTitle;
    
      histTitle = new TString("FlowReCentY_Sel");
      *histTitle += k+1;
      *histTitle += "_Har";
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistReCentY = new TProfile(histTitle->Data(),
        histTitle->Data(), 3, 0.5, 3.5);
      histFull[k].histFullHar[j].mHistReCentY->SetXTitle("FTPCE, FTPCW, TPCE, TPCW");
      histFull[k].histFullHar[j].mHistReCentY->SetYTitle("<sin n #phi>");
      delete histTitle;
    
      // Recentered, for printing, not plotting
       histTitle = new TString("FlowQreCent_Sel");
      *histTitle += k+1;
      *histTitle += "_Har";
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistQreCent = new TProfile(histTitle->Data(),
        histTitle->Data(), 2, 0.5, 2.5);
      histFull[k].histFullHar[j].mHistQreCent->SetXTitle("X, Y");
      histFull[k].histFullHar[j].mHistQreCent->SetYTitle("<Q_{n}/M>");

      // QXY
      histTitle = new TString("Flow_QXY2D_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistQXY2D = new TH2D(histTitle->Data(), histTitle->Data(),
			    nQXYBins, QXYMin, QXYMax, nQXYBins, QXYMin, QXYMax);
      histFull[k].histFullHar[j].mHistQXY2D->SetXTitle("Q_X/M");
      histFull[k].histFullHar[j].mHistQXY2D->SetYTitle("Q_Y/M");
      delete histTitle;
    
      // QSubXY for k=0 subevents
      histTitle = new TString("Flow_QFTPCSubXY2D_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistQFTPCSubXY2D = new TH2D(histTitle->Data(), histTitle->Data(),
			    nQXYBins, QXYMin*2., QXYMax*2., nQXYBins, QXYMin*2., QXYMax*2.);
      histFull[k].histFullHar[j].mHistQFTPCSubXY2D->SetXTitle("QSub_X/M");
      histFull[k].histFullHar[j].mHistQFTPCSubXY2D->SetYTitle("QSub_Y/M");
      delete histTitle;
    
      // QSubXY for k=1 subevents
      histTitle = new TString("Flow_QTPCSubXY2D_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistQTPCSubXY2D = new TH2D(histTitle->Data(), histTitle->Data(),
			    nQXYBins, QXYMin*2., QXYMax*2., nQXYBins, QXYMin*2., QXYMax*2.);
      histFull[k].histFullHar[j].mHistQTPCSubXY2D->SetXTitle("QSub_X/M");
      histFull[k].histFullHar[j].mHistQTPCSubXY2D->SetYTitle("QSub_Y/M");
      delete histTitle;
 
      // event plane difference of two selections
      histTitle = new TString("Flow_Psi_Diff_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      
      if (k == 0 ) {
	Int_t my_order = 1;
	if (j == 1) {
	  my_order = 2;
	}
	histFull[k].histFullHar[j].mHistPsi_Diff = new TH1F(histTitle->Data(),
	     histTitle->Data(), nPsiBins, -psiMax/my_order/2., psiMax/my_order/2.);
      } else {
	histFull[k].histFullHar[j].mHistPsi_Diff = new TH1F(histTitle->Data(),	
		          histTitle->Data(), nPsiBins, -psiMax/2., psiMax/2.);
      }
      
      if (k == 0) {
	if (j == 0) {
	  histFull[k].histFullHar[j].mHistPsi_Diff->SetXTitle
	    ("#Psi_{1,Sel1} - #Psi_{1,Sel2}(rad)");
	} else if (j == 1) {
	  histFull[k].histFullHar[j].mHistPsi_Diff->SetXTitle
	    ("#Psi_{2,Sel1} - #Psi_{2,Sel2}(rad)");
	}
      } else if (k == 1) {
	if (j == 0) {
	  histFull[k].histFullHar[j].mHistPsi_Diff->SetXTitle
	    ("#Psi_{1,Sel1} - #Psi_{2,Sel2}(rad)");
	} else if (j == 1) {
	  histFull[k].histFullHar[j].mHistPsi_Diff->SetXTitle
	    ("#Psi_{1,Sel1} - #Psi_{2,Sel1}(rad)");
	}
      }

      histFull[k].histFullHar[j].mHistPsi_Diff->SetYTitle("Counts");
      delete histTitle;

      // correlation of sub-event planes
      histTitle = new TString("Flow_Psi_Sub_Corr_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistPsiSubCorr = new TH1F(histTitle->Data(),
        histTitle->Data(), nPsiBins, psiMin, psiMax / order);
      histFull[k].histFullHar[j].mHistPsiSubCorr->Sumw2();
      histFull[k].histFullHar[j].mHistPsiSubCorr->SetXTitle
	("Sub-Event Correlation (rad)");
      histFull[k].histFullHar[j].mHistPsiSubCorr->SetYTitle("Counts");
      delete histTitle;
      
      // correlation of sub-event planes of different order
      histTitle = new TString("Flow_Psi_Sub_Corr_Diff_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistPsiSubCorrDiff = new 
        TH1F(histTitle->Data(), histTitle->Data(), nPsiBins, psiMin,
	     psiMax / (order+1.));
      histFull[k].histFullHar[j].mHistPsiSubCorrDiff->Sumw2();
      histFull[k].histFullHar[j].mHistPsiSubCorrDiff->SetXTitle
	("Sub-Event Correlation (rad)");
      histFull[k].histFullHar[j].mHistPsiSubCorrDiff->SetYTitle("Counts");
      delete histTitle;
      
      // q
      histTitle = new TString("Flow_q_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_q = new TH1F(histTitle->Data(),
        histTitle->Data(), n_qBins, qMin, qMax);
      histFull[k].histFullHar[j].mHist_q->Sumw2();
      histFull[k].histFullHar[j].mHist_q->SetXTitle("q = |Q|/sqrt(Mult)");
      histFull[k].histFullHar[j].mHist_q->SetYTitle("Counts");
      delete histTitle;
      
      // particle-plane azimuthal correlation
      histTitle = new TString("Flow_Phi_Corr_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistPhiCorr =	new TH1F(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax / order);
      histFull[k].histFullHar[j].mHistPhiCorr->Sumw2();
      histFull[k].histFullHar[j].mHistPhiCorr->
	SetXTitle("Particle-Plane Correlation (rad)");
      histFull[k].histFullHar[j].mHistPhiCorr->SetYTitle("Counts");
      delete histTitle;

      // Yield(eta,pt)
      histTitle = new TString("Flow_Yield2D_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHistYield2D =	new TH2D(histTitle->Data(),
        histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax,
			   Flow::nPtBins, Flow::ptMin, Flow::ptMax);
      histFull[k].histFullHar[j].mHistYield2D->Sumw2();
      histFull[k].histFullHar[j].mHistYield2D->SetXTitle("Pseudorapidty");
      histFull[k].histFullHar[j].mHistYield2D->SetYTitle("Pt (GeV/c)");
      delete histTitle;

      // Flow observed
      histTitle = new TString("Flow_vObs2D_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_vObs2D =	new TProfile2D(histTitle->Data(),
        histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax, nPtBinsPart, 
		 Flow::ptMin, ptMaxPart, -1000., 1000., "");
      histFull[k].histFullHar[j].mHist_vObs2D->SetXTitle((char*)xLabel.Data());
      histFull[k].histFullHar[j].mHist_vObs2D->SetYTitle("Pt (GeV/c)");
      delete histTitle;

      // Flow observed profiles
      histTitle = new TString("Flow_vObsEta_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_vObsEta = new TProfile(histTitle->Data(),
        histTitle->Data(), mNEtaBins, mEtaMin, mEtaMax, -1000., 1000., "");
      histFull[k].histFullHar[j].mHist_vObsEta->SetXTitle((char*)xLabel.Data());
      histFull[k].histFullHar[j].mHist_vObsEta->SetYTitle("v (%)");
      delete histTitle;

      histTitle = new TString("Flow_vObsPt_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_vObsPt = new TProfile(histTitle->Data(),
        histTitle->Data(), nPtBinsPart, Flow::ptMin, ptMaxPart, -1000., 1000., "");
      histFull[k].histFullHar[j].mHist_vObsPt->SetXTitle("Pt (GeV/c)");
      histFull[k].histFullHar[j].mHist_vObsPt->SetYTitle("v (%)");
      delete histTitle;
    }

    // for two harmonics
    for (int j = 0; j < 2; j++) {

      // Phi lab
      // Tpc (FarEast)
      histTitle = new TString("Flow_Phi_FarEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFarEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFarEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFarEast->SetYTitle("Counts");
      delete histTitle;
      
      // Tpc (East)
      histTitle = new TString("Flow_Phi_East_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiEast->SetYTitle("Counts");
      delete histTitle;
      
      // Tpc (West)
      histTitle = new TString("Flow_Phi_West_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWest->SetYTitle("Counts");
      delete histTitle;
      
      // Tpc (FarWest)
      histTitle = new TString("Flow_Phi_FarWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFarWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFarWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFarWest->SetYTitle("Counts");
      delete histTitle;
      
      // Ftpc (FarEast)
      histTitle = new TString("Flow_Phi_FtpcFarEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFtpcFarEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFtpcFarEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFtpcFarEast->SetYTitle("Counts");
      delete histTitle;
      
      // Ftpc (East)
      histTitle = new TString("Flow_Phi_FtpcEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFtpcEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFtpcEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFtpcEast->SetYTitle("Counts");
      delete histTitle;
      
      // Ftpc (West)
      histTitle = new TString("Flow_Phi_FtpcWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFtpcWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFtpcWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFtpcWest->SetYTitle("Counts");
      delete histTitle;
      
      // Ftpc (FarWest)
      histTitle = new TString("Flow_Phi_FtpcFarWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFtpcFarWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFtpcFarWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFtpcFarWest->SetYTitle("Counts");
      delete histTitle;
      
      // PhiWgt
      // Tpc (FarEast)
      histTitle = new TString("Flow_Phi_Weight_FarEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtFarEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtFarEast->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtFarEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtFarEast->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Tpc (East)
      histTitle = new TString("Flow_Phi_Weight_East_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtEast->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtEast->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Tpc (West)
      histTitle = new TString("Flow_Phi_Weight_West_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtWest->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtWest->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Tpc (FarWest)
      histTitle = new TString("Flow_Phi_Weight_FarWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtFarWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtFarWest->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtFarWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtFarWest->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Ftpc (FarEast)
      histTitle = new TString("Flow_Phi_Weight_FtpcFarEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Ftpc (East)
      histTitle = new TString("Flow_Phi_Weight_FtpcEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Ftpc (West)
      histTitle = new TString("Flow_Phi_Weight_FtpcWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest->SetYTitle("PhiWgt");
      delete histTitle;
      
      // Ftpc (FarWest)
      histTitle = new TString("Flow_Phi_Weight_FtpcFarWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest->Sumw2();
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest->SetYTitle("PhiWgt");
      delete histTitle;


      // Phi lab flattened
      // Tpc (FarEast)
      histTitle = new TString("Flow_Phi_Flat_FarEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatFarEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatFarEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatFarEast->SetYTitle("Counts");
      delete histTitle;

      // Tpc (East)
      histTitle = new TString("Flow_Phi_Flat_East_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatEast->SetYTitle("Counts");
      delete histTitle;

      // Tpc (West)
      histTitle = new TString("Flow_Phi_Flat_West_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatWest->SetYTitle("Counts");
      delete histTitle;

      // Tpc (FarWest)
      histTitle = new TString("Flow_Phi_Flat_FarWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatFarWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBins, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatFarWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatFarWest->SetYTitle("Counts");
      delete histTitle;

      // Ftpc (FarEast)
      histTitle = new TString("Flow_Phi_Flat_FtpcFarEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarEast->SetYTitle("Counts");
      delete histTitle;     

      // Ftpc (East)
      histTitle = new TString("Flow_Phi_Flat_FtpcEast_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcEast = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcEast->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcEast->SetYTitle("Counts");
      delete histTitle;
      
      // Ftpc (West)
      histTitle = new TString("Flow_Phi_Flat_FtpcWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcWest->SetYTitle("Counts");
      delete histTitle;
      
      // Ftpc (FarWest)
      histTitle = new TString("Flow_Phi_Flat_FtpcFarWest_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarWest = new TH1D(histTitle->Data(),
        histTitle->Data(), Flow::nPhiBinsFtpc, phiMin, phiMax);
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarWest->SetXTitle
	("Azimuthal Angles (rad)");
      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarWest->SetYTitle("Counts");
      delete histTitle;
   
    }
  }

  // Calculate recenter paramerters?
  TFile fileReCent("flow.reCentAna.root","R");
  if (reCentCalc) {
    gMessMgr->Info("##### FlowAnalysis: Calc reCent Pars");
    mCalcReCentPars = kTRUE;
  } else {
    mCalcReCentPars = kFALSE;
  }

  gMessMgr->SetLimit("##### FlowAnalysis", 2);
  gMessMgr->Info("##### FlowAnalysis: $Id: StFlowAnalysisMaker.cxx,v 1.103 2011/07/25 15:54:42 posk Exp $");

  //Do Psi-shift
  //# Added by jcampbell
  if(mFillPsiShift) {
    mPsiShiftFile = new TFile("flow.psiShift.root","RECREATE");
	  
	  for (int k = 0; k < Flow::nSels; k++) {
		for (int j = 0; j < Flow::nHars; j++) {
		  for (int n = 0; n < Flow::nSubs + 1; n++) {
			TString sinName = "avgSin_sel";
			TString cosName = "avgCos_sel";
			sinName += k; sinName += "_har"; sinName += j; sinName += "_sub"; sinName += n;
			cosName += k; cosName += "_har"; cosName += j; cosName += "_sub"; cosName += n;
			mPsiShiftSin[k][j][n] = new TProfile(sinName.Data(),sinName.Data(),Flow::nPsiShiftOrders,0.5,(Flow::nPsiShiftOrders+0.5));
			mPsiShiftCos[k][j][n] = new TProfile(cosName.Data(),cosName.Data(),Flow::nPsiShiftOrders,0.5,(Flow::nPsiShiftOrders+0.5));
			} // nSubs 
		} // nHars
	} // nSels
	  
  } else if(mDoPsiShift) {
    mPsiShiftFile = new TFile("flow.psiShift.root","READ");
	  for (int k = 0; k < Flow::nSels; k++) {
		for (int j = 0; j < Flow::nHars; j++) {
		  for (int n = 0; n < Flow::nSubs + 1; n++) {
			TString sinName = "avgSin_sel";
			TString cosName = "avgCos_sel";
			sinName += k; sinName += "_har"; sinName += j; sinName += "_sub"; sinName += n;
			cosName += k; cosName += "_har"; cosName += j; cosName += "_sub"; cosName += n;
			mPsiShiftSin[k][j][n] = (TProfile*)mPsiShiftFile->Get(sinName.Data());
			mPsiShiftCos[k][j][n] = (TProfile*)mPsiShiftFile->Get(cosName.Data());
			} // nSubs 
		} // nHars
	} // nSels
	  
  } 
  
  
  return StMaker::Init();
}

//-----------------------------------------------------------------------

Bool_t StFlowAnalysisMaker::FillFromFlowEvent() {
  // Get event quantities from StFlowEvent

  Bool_t kRETURN = kTRUE;

  for (int k = 0; k < Flow::nSels; k++) {
    pFlowSelect->SetSelection(k);
    for (int j = 0; j < Flow::nHars; j++) {
	  if(j == 0) continue; // SUPER UGLY HACK. BEWARE!!
      pFlowSelect->SetHarmonic(j);
      for (int n = 0; n < Flow::nSubs; n++) {
	pFlowSelect->SetSubevent(n);
	int i = Flow::nSels*k + n;
	// sub-event quantities already recentered
	mQSub[i][j] = pFlowEvent->Q(pFlowSelect);
	mPsiSub[i][j] = pFlowEvent->Psi(pFlowSelect);
 	mMultSub[i][j] = pFlowEvent->Mult(pFlowSelect);
	if (mMultSub[i][j] < 3.) kRETURN = kFALSE; // to eliminate multSub < 3
	if (mQSub[i][j].Mod()==0.) kRETURN = kFALSE; // to eliminate psiSub=0
      }//subs
      
      pFlowSelect->SetSubevent(-1);
      // full event quantities already recentered
      mQ[k][j]    = pFlowEvent->Q(pFlowSelect);
      mPsi[k][j]  = pFlowEvent->Psi(pFlowSelect);
      m_q[k][j]   = pFlowEvent->q(pFlowSelect);
      mMult[k][j] = pFlowEvent->Mult(pFlowSelect);
      if (mQ[k][j].Mod()==0.) kRETURN = kFALSE; // to eliminate psi=0	
    }//full
  }

  mZDCSMD_e_PsiWgt = pFlowEvent->ZDCSMD_PsiWgtEast();
  mZDCSMD_w_PsiWgt = pFlowEvent->ZDCSMD_PsiWgtWest();
  mZDCSMD_f_PsiWgt = (pFlowEvent->UseZDCSMD()) ? pFlowEvent->ZDCSMD_PsiWgtFull():1.;
  mFlowWeight  	   = (pFlowEvent->UseZDCSMD()) ? mZDCSMD_e_PsiWgt*mZDCSMD_w_PsiWgt*mZDCSMD_f_PsiWgt:1.;

  return kRETURN;
}

//-----------------------------------------------------------------------
//# Added by jcampbell
void StFlowAnalysisMaker::PsiShift() {
  // Parent function to handle psi-shifting 
  for (int k = 0; k < Flow::nSels; k++) {
    for (int j = 0; j < Flow::nHars; j++) {
  	    if(mFillPsiShift) {
	        int i = Flow::nSels*k;
			FillPsiShift(mPsi[k][j],k,j+1,0);
			FillPsiShift(mPsiSub[i][j],k,j+1,1);
			FillPsiShift(mPsiSub[i+1][j],k,j+1,2);
  		} else if(mDoPsiShift) {
	        int i = Flow::nSels*k;
			mPsi[k][j] += DoPsiShift(mPsi[k][j],k,j+1,0);
			mPsiSub[i][j] += DoPsiShift(mPsiSub[i][j],k,j+1,1);
			mPsiSub[i+1][j] += DoPsiShift(mPsiSub[i+1][j],k,j+1,2);
  		}
    } // nHars
  } // nSels
} 

//-----------------------------------------------------------------------
//# Added by jcampbell
void StFlowAnalysisMaker::FillPsiShift(Float_t rawPsi, Int_t sel, Int_t harm, Int_t sub) {
  // Fill Psi-shift coefficients 

	for(Int_t i = 1; i <= Flow::nPsiShiftOrders; ++i)
	{
	  mPsiShiftCos[sel][harm-1][sub]->Fill(i, cos(harm*i*rawPsi));
	  mPsiShiftSin[sel][harm-1][sub]->Fill(i, sin(harm*i*rawPsi));
	}
  
}

//-----------------------------------------------------------------------
//# Added by jcampbell
Float_t StFlowAnalysisMaker::DoPsiShift(Float_t rawPsi, Int_t sel, Int_t harm, Int_t sub) {
  // Calculate the psi-shift
 
	Float_t deltaPsi = 0, temp = 0;	
	
	for(Int_t i = 1; i <= Flow::nPsiShiftOrders; ++i)
	{
        temp = 0;
		temp += mPsiShiftCos[sel][harm-1][sub]->GetBinContent(i)*sin(harm*i*rawPsi);
		temp -= mPsiShiftSin[sel][harm-1][sub]->GetBinContent(i)*cos(harm*i*rawPsi);
		temp /= i;
        deltaPsi += temp;
	}

	deltaPsi *= (2.0 / harm);
	return deltaPsi;

}

//-----------------------------------------------------------------------
//# Added by jcampbell
void StFlowAnalysisMaker::SetFillPsiShift(Bool_t fillPsiShift) {
  // Fill histograms with event quantities
  mFillPsiShift = fillPsiShift;
}

///-----------------------------------------------------------------------
//# Added by jcampbell
Float_t  StFlowAnalysisMaker::q2(Int_t sel) {
	return m_q[sel][1];
}

///-----------------------------------------------------------------------
//# Added by jcampbell
Float_t  StFlowAnalysisMaker::Psi2(Int_t sel, Int_t sub) {
  if(sub == 0){ return mPsi[sel][1];}
  else { int i = Flow::nSels*sel + (sub - 1); return mPsiSub[i][1];}
}

//-----------------------------------------------------------------------
//# Added by jcampbell
void StFlowAnalysisMaker::SetDoPsiShift(Bool_t doPsiShift) {
  // Fill histograms with event quantities
  mDoPsiShift = doPsiShift;
}

//-----------------------------------------------------------------------

void StFlowAnalysisMaker::FillEventHistograms() {
  // Fill histograms with event quantities

  // trigger
  unsigned int triggers =  StFlowCutEvent::TriggersFound();
  mHistTrigger->Fill(triggers);

  // no selections: OrigMult, MultEta, Centrality, TotalMult, PartMult, MultOverOrig, VertexZ, VertexXY
  int origMult  = pFlowEvent->OrigMult();
  mHistOrigMult ->Fill((float)origMult);
  mHistMultEta  ->Fill((float)pFlowEvent->MultEta());
  int cent      = pFlowEvent->Centrality();
  mHistCent     ->Fill((float)cent);
  int totalMult = pFlowEvent->TrackCollection()->size();
  mHistMult     ->Fill((float)totalMult);
  UInt_t partMult = pFlowEvent->MultPart(pFlowSelect);
  mHistMultPart ->Fill((float)partMult);
  if (origMult) mHistMultOverOrig->Fill((float)totalMult / (float)origMult);

  StThreeVectorF vertex = pFlowEvent->VertexPos();
  mHistVertexZ   ->Fill(vertex.z());
  mHistVertexXY2D->Fill(vertex.x(), vertex.y());

  mHistCTBvsZDC2D->Fill(pFlowEvent->ZDCe() + pFlowEvent->ZDCw(), pFlowEvent->CTB());

  //ZDCSMD test
  for(int strip=1; strip<9; strip++) {
    mZDC_SMD_west_hori->Fill(strip,pFlowEvent->ZDCSMD(1,1,strip));
    mZDC_SMD_east_hori->Fill(strip,pFlowEvent->ZDCSMD(0,1,strip));
    if(strip==8) continue;
    mZDC_SMD_west_vert->Fill(strip,pFlowEvent->ZDCSMD(1,0,strip));
    mZDC_SMD_east_vert->Fill(strip,pFlowEvent->ZDCSMD(0,0,strip));
  }
  mHistZDCSMDPsiWgtEast->Fill(pFlowEvent->ZDCSMD_PsiEst());
  mHistZDCSMDPsiWgtWest->Fill(pFlowEvent->ZDCSMD_PsiWst());
  mHistZDCSMDPsiWgtTest->Fill(mPsi[0][0]);
  mHistZDCSMDPsiWgtFull->Fill(mPsi[0][0],mFlowWeight/mZDCSMD_f_PsiWgt);
  mHistZDCSMDPsiCorTest->Fill(pFlowEvent->ZDCSMD_PsiCorr());
  mHistZDCSMDPsiCorFull->Fill(pFlowEvent->ZDCSMD_PsiCorr(),mFlowWeight/mZDCSMD_f_PsiWgt);

  // sub-event Psi_Subs
  for (int k = 0; k < Flow::nSels; k++) {
    for (int j = 0; j < Flow::nHars; j++) {
      for (int n = 0; n < Flow::nSubs; n++) {
	int i = Flow::nSels*k + n;
	if(mPsiSub[i][j]) { histSub[i].histSubHar[j].mHistPsiSubs->Fill(mPsiSub[i][j]); }
	if (mQSub[i][j].Mod() && mMultSub[i][j]) {
		//#Added by jcampbell
	  histSub[i].histSubHar[j].mHistMultSubs->Fill((float)mMultSub[i][j]);

	  if (k==0) { // FTPC
	    double QSubx = mQSub[i][j].X() / (double)mMultSub[i][j];
	    double QSuby = mQSub[i][j].Y() / (double)mMultSub[i][j];
	    histFull[n].histFullHar[j].mHistQFTPCSubXY2D->Fill(QSubx,QSuby);
	  } else if (k==1) { // TPC
	    double QSubx = mQSub[i][j].X() / (double)mMultSub[i][j];
	    double QSuby = mQSub[i][j].Y() / (double)mMultSub[i][j];
 	    histFull[n].histFullHar[j].mHistQTPCSubXY2D->Fill(QSubx,QSuby);
	  }
	}
      }
    }
  }


  for (int k = 0; k < Flow::nSels; k++) {
    pFlowSelect->SetSelection(k);
    for (int j = 0; j < Flow::nHars; j++) {
	  if(j == 0) continue; // SUPER UGLY HACK. BEWARE!!
        int i = Flow::nSels*k;
 	    histFull[k].histFullHar[j].mHistPsiFull->Fill(mPsi[k][j]);
 	    histFull[k].histFullHar[j].mHistPsiSub2VsSub1->Fill(mPsiSub[i][j],mPsiSub[i+1][j]);
 	    histFull[k].histFullHar[j].mHistRPMultSub1VsFull->Fill(mMult[k][j],mMultSub[i][j]);
 	    histFull[k].histFullHar[j].mHistRPMultSub2VsFull->Fill(mMult[k][j],mMultSub[i+1][j]);
 	    histFull[k].histFullHar[j].mHistRPMultSub2VsSub1->Fill(mMultSub[i][j],mMultSub[i+1][j]);
 	    histFull[k].histFullHar[j].mHistResSquaredVsMult->Fill(pFlowEvent->MultEta(),cos(2*(mPsiSub[i][j]-mPsiSub[i+1][j])));
    }//full
  }//sels

  // full event Psi, PsiSubCorr, PsiSubCorrDiff, cos, mult, q, reCent
  for (int k = 0; k < Flow::nSels; k++) {
    pFlowSelect->SetSelection(k);
    for (int j = 0; j < Flow::nHars; j++) {
      pFlowSelect->SetHarmonic(j);
      float order  = (float)(j+1);

      if (mPsi[k][j]) {
	histFull[k].histFullHar[j].mHistPsi->Fill(mPsi[k][j]);
	//histFull[k].histFullHar[j].mHistPsi->Fill(mPsi[k][j],mQ[k][j].Mod());
      }
      if (mQ[k][j].Mod() && mMult[k][j]) {
	double Qx = mQ[k][j].X() / (double)mMult[k][j];
	double Qy = mQ[k][j].Y() / (double)mMult[k][j];
	histFull[k].histFullHar[j].mHistQXY2D->Fill(Qx,Qy);
      }
      if (mPsi[0][j] && mPsi[1][j]) {
	if (k < 2 && j < 2) {
	  if (k == 0) { 
	    float psi1 = mPsi[0][j];
	    float psi2 = mPsi[1][j];	  
	    float diff = psi1 - psi2;	  
	    if (diff < -twopi/2./(j+1)) {
	      diff += twopi/(j+1);
	    } else if (diff > +twopi/2./(j+1)) {
	      diff -= twopi/(j+1);
	    }	  
	    histFull[k].histFullHar[j].mHistPsi_Diff->Fill(diff);
	  } else if (k == 1) {
	    float psi1 = 0.;
	    float psi2 = 0.;
	    if (j == 0) {
	      psi1 = mPsi[0][0];
	      psi2 = mPsi[1][1];
	    }  else if (j == 1) {
	      psi1 = mPsi[0][0];
	      psi2 = mPsi[0][1];
	    }
	    float diff = psi1 - psi2;
	    diff = (TMath::Abs(diff) > twopi/2.) ? ((diff > 0.) ? -(twopi - diff) :
						    -(diff + twopi)) : diff;
	    histFull[k].histFullHar[j].mHistPsi_Diff->Fill(diff);
	  }      
	}
      }
      
      if (mPsiSub[Flow::nSels*k][j] != 0. && mPsiSub[Flow::nSels*k+1][j] != 0.) {
	float psiSubCorr;
	if(pFlowEvent->UseZDCSMD()) {
	  psiSubCorr = pFlowEvent->ZDCSMD_PsiCorr();
	} 
	else if (mV1Ep1Ep2 == kFALSE || order != 1) { // normal case
	  psiSubCorr = mPsiSub[Flow::nSels*k][j] - mPsiSub[Flow::nSels*k+1][j];
	}
	else { // mV1Ep1Ep2 == kTRUE && order == 1
	  psiSubCorr = mPsiSub[Flow::nSels*k][0] + mPsiSub[Flow::nSels*k+1][0] - 2.*mPsi[k][1];
	}
	histFull[k].mHistCos->Fill(order, (float)cos(order * psiSubCorr),mFlowWeight/mZDCSMD_f_PsiWgt);
	if (psiSubCorr < 0.) psiSubCorr += twopi / order;
	if (psiSubCorr > twopi / order) psiSubCorr -= twopi / order; // for v1Ep1Ep2 which gives -twopi < psiSubCorr < 2.*twopi
	histFull[k].histFullHar[j].mHistPsiSubCorr->Fill(psiSubCorr,mFlowWeight/mZDCSMD_f_PsiWgt);
      }

      if (j < Flow::nHars - 1) { // subevents of different harmonics
	int j1 = 0, j2 = 0;
	float psiSubCorrDiff;
	if (j==0) {
	  j1 = 1, j2 = 2;	
	} else if (j==1) {
	  j1 = 1, j2 = 3;	
	} else if (j==2) {
	  j1 = 2, j2 = 4;	
	}
	psiSubCorrDiff = fmod((double)mPsiSub[Flow::nSels*k][j1-1],
	   twopi/(double)j2) - fmod((double)mPsiSub[Flow::nSels*k+1][j2-1],
				    twopi/(double)j2);
	if (psiSubCorrDiff < 0.) psiSubCorrDiff += twopi/(float)j2;
	if (psiSubCorrDiff) { histFull[k].histFullHar[j].mHistPsiSubCorrDiff->Fill(psiSubCorrDiff); }
	psiSubCorrDiff = fmod((double)mPsiSub[Flow::nSels*k][j2-1],
	   twopi/(double)j2) - fmod((double)mPsiSub[Flow::nSels*k+1][j1-1],
				    twopi/(double)j2);
	if (psiSubCorrDiff < 0.) psiSubCorrDiff += twopi/(float)j2;
	if (psiSubCorrDiff) { histFull[k].histFullHar[j].mHistPsiSubCorrDiff->Fill(psiSubCorrDiff); }
      }

      histFull[k].histFullHar[j].mHistMult->Fill((float)mMult[k][j]);
      if (m_q[k][j]) { histFull[k].histFullHar[j].mHist_q->Fill(m_q[k][j]); }

      if (mCalcReCentPars) {
	// calculate recentering parameters, fill 4 bins
	TVector2 qReCent;
	qReCent = pFlowEvent->ReCentEPPar(pFlowSelect,"FTPCE");
	if (qReCent.X()) histFull[k].histFullHar[j].mHistReCentX->Fill(1., qReCent.X());
	if (qReCent.Y()) histFull[k].histFullHar[j].mHistReCentY->Fill(1., qReCent.Y());
	qReCent = pFlowEvent->ReCentEPPar(pFlowSelect,"FTPCW");
	if (qReCent.X()) histFull[k].histFullHar[j].mHistReCentX->Fill(2., qReCent.X());
	if (qReCent.Y()) histFull[k].histFullHar[j].mHistReCentY->Fill(2., qReCent.Y());
	qReCent = pFlowEvent->ReCentEPPar(pFlowSelect,"TPCE");
	if (qReCent.X()) histFull[k].histFullHar[j].mHistReCentX->Fill(3., qReCent.X());
	if (qReCent.Y()) histFull[k].histFullHar[j].mHistReCentY->Fill(3., qReCent.Y());
	qReCent = pFlowEvent->ReCentEPPar(pFlowSelect,"TPCW");
	if (qReCent.X()) histFull[k].histFullHar[j].mHistReCentX->Fill(4., qReCent.X());
	if (qReCent.Y()) histFull[k].histFullHar[j].mHistReCentY->Fill(4., qReCent.Y());
      }
    }
  }
}

//-----------------------------------------------------------------------

void StFlowAnalysisMaker::FillParticleHistograms() {
  // Fill histograms from the particles
  float etaSymPosTpcN  = 0.;
  float etaSymNegTpcN  = 0.;
  float etaSymPosFtpcN = 0.;
  float etaSymNegFtpcN = 0.;
  float hPlusN         = 0.;
  float hMinusN        = 0.;
  float piPlusN        = 0.;
  float piMinusN       = 0.;
  float protonN        = 0.;
  float pbarN          = 0.;
  float kMinusN        = 0.;
  float kPlusN         = 0.;
  float deuteronN      = 0.;
  float dbarN          = 0.;
  float electronN      = 0.;
  float positronN      = 0.;

  // Initialize Iterator
  StFlowTrackCollection* pFlowTracks = pFlowEvent->TrackCollection();
  StFlowTrackIterator itr;

  for (itr = pFlowTracks->begin(); itr != pFlowTracks->end(); itr++) {
    StFlowTrack* pFlowTrack = *itr;

    float phi         = pFlowTrack->Phi();
    if (phi < 0.) phi += twopi;
    float eta         = pFlowTrack->Eta();
    float zFirstPoint = 0.;
    float zLastPoint  = 0.;
    if (pFlowEvent->FirstLastPoints()) {
      zFirstPoint = pFlowTrack->ZFirstPoint();
      zLastPoint  = pFlowTrack->ZLastPoint();
    }
    float pt          = pFlowTrack->Pt();
    int   charge      = pFlowTrack->Charge();
    float dca         = pFlowTrack->Dca();
    float dcaGlobal   = pFlowTrack->DcaGlobal();
    float chi2        = pFlowTrack->Chi2();
    int   fitPts      = pFlowTrack->FitPts();
    int   maxPts      = pFlowTrack->MaxPts();
    Char_t pid[10];
    strcpy(pid, pFlowTrack->Pid());
    float totalp    = pFlowTrack->P();
    float logp      = ::log(totalp);
    float dEdx      = pFlowTrack->Dedx();
    StTrackTopologyMap map = pFlowTrack->TopologyMap();

    // no selections: Charge, Dca, DcaGlobal, Chi2, FitPts, MaxPts, FitOverMax, PID

    // distinguish between Tpc and Ftpc
    if (map.trackFtpcEast() || map.trackFtpcWest()) {
      mHistChargeFtpc   ->Fill((float)charge);
      mHistDcaFtpc      ->Fill(dca);
      mHistDcaGlobalFtpc->Fill(dcaGlobal);
      mHistChi2Ftpc     ->Fill(chi2);
      mHistFitPtsFtpc   ->Fill((float)fitPts);
      mHistMaxPtsFtpc   ->Fill((float)maxPts);
      if (maxPts) mHistFitOverMaxFtpc->Fill((float)fitPts/(float)maxPts);
    }
    
    else { // Tpc track or otherwise!!!
      
      // For PID multiplicites
      if (strcmp(pid, "pi+")  == 0)  piPlusN++;
      if (strcmp(pid, "pi-")  == 0)  piMinusN++;
      if (strcmp(pid, "pr+")  == 0)  protonN++;
      if (strcmp(pid, "pr-")  == 0)  pbarN++;
      if (strcmp(pid, "k+")   == 0)  kPlusN++;
      if (strcmp(pid, "k-")   == 0)  kMinusN++;
      if (strcmp(pid, "d+")   == 0)  deuteronN++;
      if (strcmp(pid, "d-")   == 0)  dbarN++;
      if (strcmp(pid, "e-")   == 0)  electronN++;
      if (strcmp(pid, "e+")   == 0)  positronN++;
      
      mHistDcaTpc      ->Fill(dca);
      mHistDcaGlobalTpc->Fill(dcaGlobal);
      mHistChi2Tpc     ->Fill(chi2);
      mHistFitPtsTpc   ->Fill((float)fitPts);
      mHistMaxPtsTpc   ->Fill((float)maxPts);
      if (maxPts) mHistFitOverMaxTpc->Fill((float)fitPts/(float)maxPts);
      
      if (charge == 1) {
	hPlusN++;
	mHistMeanDedxPos2D->Fill(logp, dEdx);
	float piPlus = pFlowTrack->PidPiPlus();
	mHistPidPiPlus->Fill(piPlus);
	if (strcmp(pid, "pi+") == 0) {
	  mHistMeanDedxPiPlus2D->Fill(logp, dEdx);
	  mHistPidPiPlusPart->Fill(piPlus);
	}
	float kplus  = pFlowTrack->PidKaonPlus();
	mHistPidKplus->Fill(kplus);
	if (strcmp(pid, "k+") == 0) {
	  mHistMeanDedxKplus2D->Fill(logp, dEdx);
	  mHistPidKplusPart->Fill(kplus);
	}
	float proton  = pFlowTrack->PidProton();
	mHistPidProton->Fill(proton);
	if (strcmp(pid, "pr+") == 0) {
	  mHistMeanDedxProton2D->Fill(logp, dEdx);
	  mHistPidProtonPart->Fill(proton);
	}
	float deuteron  = pFlowTrack->PidDeuteron();
	mHistPidDeuteron->Fill(deuteron);
	if (strcmp(pid, "d+") == 0) {
	  mHistMeanDedxDeuteron2D->Fill(logp, dEdx);
	  mHistPidDeuteronPart->Fill(deuteron);
	}
	float positron  = pFlowTrack->PidPositron();
	mHistPidPositron->Fill(positron);
	if (strcmp(pid, "e+") == 0) {
	  mHistMeanDedxPositron2D->Fill(logp, dEdx);
	  mHistPidPositronPart->Fill(positron);
	  }
      } else if (charge == -1) {
	hMinusN++;
	mHistMeanDedxNeg2D->Fill(logp, dEdx);
	float piMinus = pFlowTrack->PidPiMinus();
	mHistPidPiMinus->Fill(piMinus);
	if (strcmp(pid, "pi-") == 0) {
	  mHistMeanDedxPiMinus2D->Fill(logp, dEdx);
	  mHistPidPiMinusPart->Fill(piMinus);
	}
	float kminus  = pFlowTrack->PidKaonMinus();
	mHistPidKminus->Fill(kminus);
	if (strcmp(pid, "k-") == 0) {
	  mHistMeanDedxKminus2D->Fill(logp, dEdx);
	  mHistPidKminusPart->Fill(kminus);
	}
	float antiproton  = pFlowTrack->PidAntiProton();
	mHistPidAntiProton->Fill(antiproton);
	if (strcmp(pid, "pr-") == 0) {
	  mHistMeanDedxPbar2D->Fill(logp, dEdx);
	  mHistPidAntiProtonPart->Fill(antiproton);
	}
	float antideuteron  = pFlowTrack->PidAntiDeuteron();
	mHistPidAntiDeuteron->Fill(antideuteron);
	if (strcmp(pid, "d-") == 0) {
	  mHistMeanDedxAntiDeuteron2D->Fill(logp, dEdx);
	  mHistPidAntiDeuteronPart->Fill(antideuteron);
	}
	float electron  = pFlowTrack->PidElectron();
	mHistPidElectron->Fill(electron);
	if (strcmp(pid, "e-") == 0) {
	  mHistMeanDedxElectron2D->Fill(logp, dEdx);
	  mHistPidElectronPart->Fill(electron);
	}
      }
    }//all tracks
    
    // Yield3D, Yield2D, BinEta, BinPt
    mHistEtaPtPhi3D->Fill(eta, pt, phi);
    mHistYieldAll2D->Fill(eta, pt);
    if (pFlowSelect->SelectPart(pFlowTrack)) {
      if (strlen(pFlowSelect->PidPart()) != 0) { 
	float rapidity = pFlowTrack->Y();
	mHistBinEta->Fill(rapidity, rapidity);
	mHistYieldPart2D->Fill(rapidity, pt);
      } else {
	mHistBinEta->Fill(eta, eta);
	mHistYieldPart2D->Fill(eta, pt);
      }
      mHistBinPt->Fill(pt, pt);
    }//particles correlated with the EP

    // For Eta symmetry
    // Tpc
    if (map.hasHitInDetector(kTpcId)) {
      (eta > 0.) ? etaSymPosTpcN++ : etaSymNegTpcN++;
    }
    // Ftpc
    else if (map.trackFtpcEast() || map.trackFtpcWest()) {
      (eta > 0.) ? etaSymPosFtpcN++ : etaSymNegFtpcN++;
    }

    TVector2 Q, reCent;
    Double_t mult;
    // Get the angle of the Event Plane
    for (int k = 0; k < Flow::nSels; k++) {
      pFlowSelect->SetSelection(k);
      for (int j = 0; j < Flow::nHars; j++) {
	bool oddHar = (j+1) % 2;
	pFlowSelect->SetHarmonic(j);
	double order  = (double)(j+1);
	double orderEP  = order;
	float psi_i = 0., psi_2 = 0.;
	if (pFlowEvent->EtaSubs()) { // particles with the other subevent
	  int i = Flow::nSels*k;
	  psi_i = (eta > 0.) ? mPsiSub[i+1][j] : mPsiSub[i][j];
	} else if (pFlowEvent->RanSubs()) { // particles with the other subevent
	  int i = Flow::nSels*k;
	  if (pFlowTrack->Select(j,k,0)) {
	    psi_i =  mPsiSub[i+1][j];
	  } else if (pFlowTrack->Select(j,k,1)) {
	    psi_i =  mPsiSub[i][j];
	  } else { // neither
	    int r = (eta > 0.) ? 1 : 0;
	    psi_i =  mPsiSub[i+r][j]; // random
	  }
	} else if (order > 3. && !oddHar) { // 4, 6, etc
	  psi_i = mPsi[k][1];  // 2nd harmomic event plane
	  if (psi_i > twopi/order) psi_i -= twopi/order; // ???
	  if (psi_i > twopi/order) psi_i -= twopi/order;
	} else {
	  psi_i = mPsi[k][j];
	} // full EP

	if (pFlowSelect->Select(pFlowTrack)) { // particles used for the EP

	  histFull[k].histFullHar[j].mHistYield2D->Fill(eta, pt);
	  histFull[k].histFullHar[j].mHistdEdXVsRigidity->Fill(charge*totalp, dEdx);

	  // Get phiWgt and reCent from files
	  double phiWgt = pFlowEvent->PhiWeight(k, j, pFlowTrack);
	  TVector2 reCent = pFlowEvent->ReCentEP(k, j, pFlowTrack);

//	  if (pFlowMaker->PhiWgtCalc() && j < 2) { // only first two harmonics for phiWgt
	  if (j < 2) { // only first two harmonics for phiWgt

	    // Get detID
	    Bool_t kTpcFarEast  = kFALSE;
	    Bool_t kTpcEast     = kFALSE;
	    Bool_t kTpcWest     = kFALSE;
	    Bool_t kTpcFarWest  = kFALSE;
	    Bool_t kFtpcFarEast = kFALSE;
	    Bool_t kFtpcEast    = kFALSE;
	    Bool_t kFtpcWest    = kFALSE;
	    Bool_t kFtpcFarWest = kFALSE;
	    if (map.hasHitInDetector(kTpcId) || (map.data(0) == 0 && map.data(1) == 0)) {
	      // Tpc track, or TopologyMap not available
	      // Set TpcEast and West
	      if (pFlowEvent->FirstLastPoints()) {
		if (zFirstPoint > 0. && zLastPoint > 0.) {
		  kTpcFarWest = kTRUE;
		} else if (zFirstPoint > 0. && zLastPoint < 0.) {
		  kTpcWest = kTRUE;
		} else if (zFirstPoint < 0. && zLastPoint > 0.) {
		  kTpcEast = kTRUE;
		} else {
		  kTpcFarEast = kTRUE;
		}
	      } else {
		Float_t vertexZ = pFlowEvent->VertexPos().z();
		if (eta > 0. && vertexZ > 0.) {
		  kTpcFarWest = kTRUE;
		} else if (eta > 0. && vertexZ < 0.) {
		  kTpcWest = kTRUE;
		} else if (eta < 0. && vertexZ > 0.) {
		  kTpcEast = kTRUE;
		} else {
		  kTpcFarEast = kTRUE;
		}
	      }
	    } else if (map.trackFtpcEast()) { // FTPC track
	      Float_t vertexZ = pFlowEvent->VertexPos().z();
	      if (vertexZ > 0.) {
		kFtpcEast = kTRUE;
	      } else { // vertexZ < 0.
		kFtpcFarEast = kTRUE;
	      }
	    } else if (map.trackFtpcWest()) {
	      Float_t vertexZ = pFlowEvent->VertexPos().z();
	      if (vertexZ > 0.) {
		kFtpcFarWest = kTRUE;
	      } else { // vertexZ > 0.
		kFtpcWest = kTRUE;
	      }
	    }	    

	    // Calculate weights for filling histograms
	    float wt = 1.;
	    if (pFlowEvent->PtWgt()) { // pt wgt
	      wt *= (pt < pFlowEvent->PtWgtSaturation()) ? pt : 
		pFlowEvent->PtWgtSaturation();  // pt weighting going constant
	    }
	    float etaAbs = fabs(eta); // etaWgt, eta > 1.
	    //if (pFlowEvent->EtaWgt() && oddHar && etaAbs > 1.) { wt *= etaAbs; }
	    if (pFlowEvent->EtaWgt() && j==0 && etaAbs > 1.) { wt *= etaAbs; }
	    
	    // Fill histograms with selections
	    if (kFtpcFarEast) {
	      histFull[k].histTwoHar[j].mHistPhiFtpcFarEast->Fill(phi,wt);
	    } else if (kFtpcEast) {
	      histFull[k].histTwoHar[j].mHistPhiFtpcEast->Fill(phi,wt);
	    } else if (kFtpcWest) {
	      histFull[k].histTwoHar[j].mHistPhiFtpcWest->Fill(phi,wt);
	    } else if (kFtpcFarWest) {
	      histFull[k].histTwoHar[j].mHistPhiFtpcFarWest->Fill(phi,wt);
	    } else if (kTpcFarEast){
	      histFull[k].histTwoHar[j].mHistPhiFarEast->Fill(phi,wt);
	    } else if (kTpcEast){
	      histFull[k].histTwoHar[j].mHistPhiEast->Fill(phi,wt);
	    } else if (kTpcWest){
	      histFull[k].histTwoHar[j].mHistPhiWest->Fill(phi,wt);
	    } else if (kTpcFarWest){
	      histFull[k].histTwoHar[j].mHistPhiFarWest->Fill(phi,wt);
	    }
	    
	    // Which flat hist?
	    if (pFlowEvent->FirstLastPoints() && !pFlowEvent->FirstLastPhiWgt()) {
	      kTpcFarEast = kFALSE;
	      kTpcEast    = kFALSE;
	      kTpcWest    = kFALSE;
	      kTpcFarWest = kFALSE;
	      Float_t vertexZ = pFlowEvent->VertexPos().z();
	      if (eta > 0. && vertexZ > 0.) {
		kTpcFarWest = kTRUE;
	      } else if (eta > 0. && vertexZ < 0.) {
		kTpcWest = kTRUE;
	      } else if (eta < 0. && vertexZ > 0.) {
		kTpcEast = kTRUE;
	      } else {
		kTpcFarEast = kTRUE;
	      }
	    }
	    
	    //if (oddHar && eta < 0.) phiWgt /= -1.; // only for flat hists
	    if (j==0 && eta < 0.) phiWgt /= -1.; // only for flat hists
	    // Fill Flat histograms
	    if (kFtpcFarEast) {
	      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarEast->Fill(phi, phiWgt);
	    } else if (kFtpcEast) {
	      histFull[k].histTwoHar[j].mHistPhiFlatFtpcEast->Fill(phi, phiWgt);
	    } else if (kFtpcWest) {
	      histFull[k].histTwoHar[j].mHistPhiFlatFtpcWest->Fill(phi, phiWgt);
	    } else if (kFtpcFarWest) {
	      histFull[k].histTwoHar[j].mHistPhiFlatFtpcFarWest->Fill(phi, phiWgt);
	    } else if (kTpcFarEast) {
	      histFull[k].histTwoHar[j].mHistPhiFlatFarEast->Fill(phi, phiWgt);
	    } else if (kTpcEast) {
	      histFull[k].histTwoHar[j].mHistPhiFlatEast->Fill(phi, phiWgt);
	    } else if (kTpcWest) {
	      histFull[k].histTwoHar[j].mHistPhiFlatWest->Fill(phi, phiWgt);
	    } else if (kTpcFarWest) {
	      histFull[k].histTwoHar[j].mHistPhiFlatFarWest->Fill(phi, phiWgt);
	    }
	    if (j==0 && eta < 0.) phiWgt *= -1.; // restore value	    
	  }	  

	  // Remove autocorrelations with full EP
	  if (!pFlowEvent->EtaSubs() && !pFlowEvent->RanSubs()) {
	    TVector2 Q_i;
	    if (order > 3. && !oddHar) { // 2nd harmonic event plane
	      orderEP = 2;
	    }
	    double Qx = phiWgt * cos(orderEP * phi) - reCent.X();
	    double Qy = phiWgt * sin(orderEP * phi) - reCent.Y();
	    Q_i.Set(Qx, Qy);
	    TVector2 mQ_i = mQ[k][j] - Q_i;
	    psi_i = mQ_i.Phi() / orderEP;
	    if (psi_i < 0.) psi_i += twopi / orderEP;
	  }
	}//particles used for the EP
	  
	// Remove autocorrelations of the 2nd order 'particles' which were used for v1{EP1,EP2}.
	if (mV1Ep1Ep2 == kTRUE && order == 1) {
	  StFlowSelection usedForPsi2 = *pFlowSelect;
	  usedForPsi2.SetHarmonic(1);
	  usedForPsi2.SetSubevent(-1);
	  if (usedForPsi2.Select(pFlowTrack)) {
	    TVector2 Q_i;
	    double phiWgt = pFlowEvent->PhiWeight(k, 1, pFlowTrack);
	    TVector2 reCent = pFlowEvent->ReCentEP(k, 1, pFlowTrack);
	    double Qx = phiWgt * cos(2 * phi) - reCent.X();
	    double Qy = phiWgt * sin(2 * phi) - reCent.Y();
	    Q_i.Set(Qx, Qy);
	    TVector2 mQ_i = mQ[k][1] - Q_i;
	    psi_2 = mQ_i.Phi() / 2.;
	    if (psi_2 < 0.) psi_2 += twopi / 2.;	  
	  }
	  else { // particle was not used for Psi2 
	    psi_2 = mPsi[k][1];
	  }
	}

	// test recentering of Q per particle
	mult = (double)(pFlowEvent->Mult(pFlowSelect));
	if (mult > 0.) {
	  histFull[k].histFullHar[j].mHistQreCent->Fill(1., mQ[k][j].X()/mult);
	  histFull[k].histFullHar[j].mHistQreCent->Fill(2., mQ[k][j].Y()/mult);
	}

      	// Calculate v for all particles selected for correlation analysis
	if (pFlowSelect->SelectPart(pFlowTrack)) { // particles correlated with the EP
	  float v;
	  if (pFlowEvent->UseZDCSMD()) {
	    v = cos(order *(phi-mQ[k][1].Phi()))/perCent;
	  }
	  else if (mV1Ep1Ep2 == kFALSE || order != 1) {
	    v = cos(order * (phi - psi_i))/perCent; // normal method
	  }
	  else { // mV1Ep1Ep2 == kTRUE && order == 1
	    v = cos(phi + psi_i - 2.*psi_2)/perCent;
	  }
	  float vFlip = v;
	  if (eta < 0 && j==0) vFlip *= -1; // for 1st harmonic only
	  if (strlen(pFlowSelect->PidPart()) != 0) { // pid, fill rapidity 
	    float rapidity = pFlowTrack->Y();
	    histFull[k].histFullHar[j].mHist_vObs2D->Fill(rapidity, pt, v,mFlowWeight);
	    
	    if (mPtRange_for_vEta[1] > mPtRange_for_vEta[0]) { // cut is used
	      if (pt < mPtRange_for_vEta[1] && pt >= mPtRange_for_vEta[0]) {
		// check cut range, fill if in range
		histFull[k].histFullHar[j].mHist_vObsEta->Fill(rapidity, v,mFlowWeight);
	      }
	    }
	    else { // cut is not used, fill in any case
	      histFull[k].histFullHar[j].mHist_vObsEta->Fill(rapidity, v, mFlowWeight);
	    }

	  } else { // no pid, fill eta
	    histFull[k].histFullHar[j].mHist_vObs2D->Fill(eta, pt, v,mFlowWeight);
	   
	    if (mPtRange_for_vEta[1] > mPtRange_for_vEta[0]) { // cut is used
	      if (pt < mPtRange_for_vEta[1] && pt >= mPtRange_for_vEta[0]) {
		// check cut range, fill if in range
		histFull[k].histFullHar[j].mHist_vObsEta->Fill(eta, v, mFlowWeight);
	      }
	    }
	    else { // cut is not used, fill in any case
	      histFull[k].histFullHar[j].mHist_vObsEta->Fill(eta, v, mFlowWeight);
	    }
	  }

	  if (mEtaRange_for_vPt[1] > mEtaRange_for_vPt[0]) { // cut is used
	    if (TMath::Abs(eta) < mEtaRange_for_vPt[1] && TMath::Abs(eta) >= mEtaRange_for_vPt[0]) {
	      // check cut range, fill if in range
	      histFull[k].histFullHar[j].mHist_vObsPt->Fill(pt, vFlip, mFlowWeight);
	    }
	  }
	  else { // cut is not used, fill in any case
	    histFull[k].histFullHar[j].mHist_vObsPt->Fill(pt, vFlip,mFlowWeight);
	  }

	  // v_
	  Bool_t etaPtNoCut = kTRUE;
	  if (mPtRange_for_vEta[1] > mPtRange_for_vEta[0] && 
	      (pt < mPtRange_for_vEta[0] || pt >= mPtRange_for_vEta[1])) {
	    etaPtNoCut = kFALSE;
	  }
	  if (mEtaRange_for_vPt[1] > mEtaRange_for_vPt[0] && 
	      (TMath::Abs(eta) < mEtaRange_for_vPt[0] || 
	       TMath::Abs(eta) >= mEtaRange_for_vPt[1])) {
	    etaPtNoCut = kFALSE;
	  }
	  if (etaPtNoCut) histFull[k].mHist_vObs->Fill(order, vFlip,mFlowWeight);
	  
	  // PhiLab and Correlation of Phi of all particles with Psi
	  if (mPsi[k][j]) {
	    histFull[k].histFullHar[j].mHistPhiLab->Fill(phi);
	    histFull[k].histFullHar[j].mHistPhiLab->Fill(phi);
	    float phi_i = phi;
	    if (eta < 0 && j==0) {
	      phi_i += pi; // backward particle and 1st harmonic
	      if (phi_i > twopi) phi_i -= twopi;
	    }
	    float dPhi = phi_i - psi_i;
	    if (dPhi < 0.) dPhi += twopi;
	    histFull[k].histFullHar[j].mHistPhiCorr->
	      Fill(fmod((double)dPhi, twopi / order));
	  }
	}//particles correlated with the EP
      }
    }  
  }

  // EtaSym
  float etaSymTpc  = (etaSymPosTpcN  - etaSymNegTpcN)  / (etaSymPosTpcN  + etaSymNegTpcN);
  float etaSymFtpc = (etaSymPosFtpcN - etaSymNegFtpcN) / (etaSymPosFtpcN + etaSymNegFtpcN);

  StThreeVectorF vertex = pFlowEvent->VertexPos();
  Float_t vertexZ = vertex.z();
  mHistEtaSymVerZ2DTpc ->Fill(vertexZ , etaSymTpc);
  mHistEtaSymVerZ2DFtpc->Fill(vertexZ , etaSymFtpc);
  
  // Tpc
  float etaSymZInterceptTpc = 0.00023;  // new values introduced for 200 GeV
  float etaSymZSlopeTpc     = -0.00394; // data based on full statistics
  etaSymTpc -= (etaSymZInterceptTpc + etaSymZSlopeTpc * vertexZ); // corrected for acceptance
  etaSymTpc *= ::sqrt((etaSymPosTpcN + etaSymNegTpcN)); // corrected for statistics
  mHistEtaSymTpc->Fill(etaSymTpc);

  // Ftpc
  float etaSymZInterceptFtpc = -0.0077; // values for the FTPC based on 200 GeV data with
  float etaSymZSlopeFtpc = 0.0020;      // all sectors and 'bad runs' (323-325) excluded
  etaSymFtpc -= (etaSymZInterceptFtpc + etaSymZSlopeFtpc * vertexZ); // corrected for acceptance
  etaSymFtpc *= ::sqrt((etaSymPosFtpcN + etaSymNegFtpcN));  // corrected for statistics
  mHistEtaSymFtpc->Fill(etaSymFtpc);

  // PID multiplicities
  float totalMult = (float)pFlowEvent->TrackCollection()->size();
  mHistPidMult->Fill(1.,  totalMult);
  mHistPidMult->Fill(2.,  hPlusN);
  mHistPidMult->Fill(3.,  hMinusN);
  mHistPidMult->Fill(4.,  piPlusN);
  mHistPidMult->Fill(5.,  piMinusN);
  mHistPidMult->Fill(6.,  protonN);
  mHistPidMult->Fill(7.,  pbarN);
  mHistPidMult->Fill(8.,  kPlusN);
  mHistPidMult->Fill(9.,  kMinusN);
  mHistPidMult->Fill(10., deuteronN);
  mHistPidMult->Fill(11., dbarN);
  mHistPidMult->Fill(12., electronN);
  mHistPidMult->Fill(13., positronN);

}

//-----------------------------------------------------------------------

static Double_t resEventPlane(double chi) {
  // Calculates the event plane resolution as a function of chi

  double con = 0.626657;                   // sqrt(pi/2)/2
  double arg = chi * chi / 4.;

  Double_t res = con * chi * exp(-arg) * (TMath::BesselI0(arg) +
					  TMath::BesselI1(arg)); 

  return res;
}

//-----------------------------------------------------------------------

static Double_t resEventPlaneK2(double chi) {
  // Calculates the event plane resolution as a function of chi
  //  for the case k=2.

  double con = 0.626657;                   // sqrt(pi/2)/2
  double arg = chi * chi / 4.;

  double besselOneHalf    = ::sqrt(arg/halfpi) * sinh(arg)/arg;
  double besselThreeHalfs = ::sqrt(arg/halfpi) * (cosh(arg)/arg - sinh(arg)/(arg*arg));
  Double_t res = con * chi * exp(-arg) * (besselOneHalf + besselThreeHalfs); 

  return res;
}

//-----------------------------------------------------------------------

static Double_t resEventPlaneK3(double chi) {
  // Calculates the event plane resolution as a function of chi
  //  for the case k=3.

  double con = 0.626657;                   // sqrt(pi/2)/2
  double arg = chi * chi / 4.;

  Double_t res = con * chi * exp(-arg) * (TMath::BesselI1(arg) +
					  TMath::BesselI(2, arg)); 

  return res;
}

//-----------------------------------------------------------------------

static Double_t resEventPlaneK4(double chi) {
  // Calculates the event plane resolution as a function of chi
  //  for the case k=4.

  double con = 0.626657;                   // sqrt(pi/2)/2
  double arg = chi * chi / 4.;

  double besselOneHalf    = ::sqrt(arg/halfpi) * sinh(arg)/arg;
  double besselThreeHalfs = ::sqrt(arg/halfpi) * (cosh(arg)/arg - sinh(arg)/(arg*arg));
  double besselFiveHalfs  = besselOneHalf - 3*besselThreeHalfs/arg;

  Double_t res = con * chi * exp(-arg) * (besselThreeHalfs + besselFiveHalfs); 

  return res;
}

//-----------------------------------------------------------------------

Double_t chi(double res) {
  // Calculates chi from the event plane resolution

  double chi = 2.0;
  double delta = 1.0;

  for (int i = 0; i < 20; i++) {
      while(resEventPlane(chi) < res) {chi += delta;}
      delta = delta / 2.;
      while(resEventPlane(chi) > res) {chi -= delta;}
      delta = delta / 2.;
  }

  return chi;
}

//-----------------------------------------------------------------------

Int_t StFlowAnalysisMaker::Finish() {
  // Calculates resolution and mean flow values
  // Outputs phiWgt and reCent values

  cout << endl << "##### Analysis Maker:" << endl;
  TString* histTitle;

  // PhiWgt histogram collection
  TOrdCollection* phiWgtHistNames = new TOrdCollection(Flow::nSels*Flow::nHars+3);

  // If mCalcReCentPars write out the recentering parameters and print the recentered values
  TOrdCollection* savedHistReCentNames  = new TOrdCollection(Flow::nSels * Flow::nHars * 2);
  if (mCalcReCentPars) {
    Float_t reCentX, reCentY;
    cout << "ReCentered Q vector per particle:" << endl;
    for (int k = 0; k < Flow::nSels; k++) {
      for (int j = 0; j < Flow::nHars; j++) {
	savedHistReCentNames->AddLast(histFull[k].histFullHar[j].mHistReCentX);
	savedHistReCentNames->AddLast(histFull[k].histFullHar[j].mHistReCentY);
	reCentX   = histFull[k].histFullHar[j].mHistQreCent->GetBinContent(1);
	reCentY   = histFull[k].histFullHar[j].mHistQreCent->GetBinContent(2);
	cout << setprecision(3) << "Sel = " << k+1 << ", Har = " << j+1 << " : reCentedQ_x = "
	     << reCentX << ",\t reCentedQ_y = " << reCentY << endl;
      }
    }
  }
  cout << endl;

  // Calculate resolution from sqrt(mHistCos)
  double cosPair[Flow::nSels][Flow::nHars];
  double cosPairErr[Flow::nSels][Flow::nHars];
  double content;
  double error;
  double totalError; for (int k = 0; k < Flow::nSels; k++) {
    // Create the 1D v histogram
    histTitle = new TString("Flow_v_Sel");
    *histTitle += k+1;
    histFull[k].mHist_v = 
      histFull[k].mHist_vObs->ProjectionX(histTitle->Data());
    histFull[k].mHist_v->SetTitle(histTitle->Data());
    histFull[k].mHist_v->SetXTitle("Harmonic");
    histFull[k].mHist_v->SetYTitle("v (%)");
    delete histTitle;
    AddHist(histFull[k].mHist_v);

    for (int j = 0; j < Flow::nHars; j++) {
      double order  = (double)(j+1);
      cosPair[k][j]    = histFull[k].mHistCos->GetBinContent(j+1);
      cosPairErr[k][j] = histFull[k].mHistCos->GetBinError(j+1);

      if (pFlowEvent->UseZDCSMD()) { // ZDCSMD used to determine RP resolution
	double ZDCSMD_deltaResSub = 0.005,ZDCSMD_mResDelta=0.;
	double ZDCSMD_resSub = (histFull[k].mHistCos->GetBinContent(1)>0.) ? 
	  ::sqrt(histFull[k].mHistCos->GetBinContent(1)) : 0.;
	double ZDCSMD_resSubErr = (histFull[k].mHistCos->GetBinContent(1)>0.) ? 
	  histFull[k].mHistCos->GetBinError(1)/(2.*ZDCSMD_resSub) : 0.;
	double ZDCSMD_chiSub = chi(ZDCSMD_resSub);
	double ZDCSMD_chiSubDelta = chi((ZDCSMD_resSub+ZDCSMD_deltaResSub));
	if (j==0) {
	  mRes[k][j] = resEventPlane(::sqrt(2.) * ZDCSMD_chiSub);
	  ZDCSMD_mResDelta = resEventPlane(::sqrt(2.) * ZDCSMD_chiSubDelta);
	}
	if (j==1) {
	  mRes[k][j] = resEventPlaneK2(::sqrt(2.) * ZDCSMD_chiSub);
	  ZDCSMD_mResDelta = resEventPlaneK2(::sqrt(2.) * ZDCSMD_chiSubDelta);
	}
	if (j==2) {
	  mRes[k][j] = resEventPlaneK3(::sqrt(2.) * ZDCSMD_chiSub);
	  ZDCSMD_mResDelta = resEventPlaneK3(::sqrt(2.) * ZDCSMD_chiSubDelta);
	}
	if (j==3) {
	  mRes[k][j] = resEventPlaneK4(::sqrt(2.) * ZDCSMD_chiSub);
	  ZDCSMD_mResDelta = resEventPlaneK4(::sqrt(2.) * ZDCSMD_chiSubDelta);
	}
	mResErr[k][j] = ZDCSMD_resSubErr * fabs ((double)mRes[k][j]
		- ZDCSMD_mResDelta) / ZDCSMD_deltaResSub;
	
      }//UseZDCSMD
      else {
	if (cosPair[k][j] > 0.) {
	  double resSub, resSubErr;
	  double res2, res2error;
	  if (mV1Ep1Ep2 == kTRUE && order == 1) { // mixed harmonics
	    // calculate resolution of second order event plane first
	    if (histFull[k].mHistCos->GetBinContent(2) > 0.) {
	      if (histFull[k].mHistCos->GetBinContent(2) > 0.92) { // resolution saturates
		res2 = 0.99;
		res2error = 0.007;
	      } else {
		double deltaRes2Sub = 0.005;  // differential for the error propagation
		double res2Sub = ::sqrt(histFull[k].mHistCos->GetBinContent(2));
		double res2SubErr = histFull[k].mHistCos->GetBinError(2) / (2. * res2Sub);
		double chiSub2 = chi(res2Sub);
		double chiSub2Delta = chi(res2Sub + deltaRes2Sub);
		res2 = resEventPlane(::sqrt(2.) * chiSub2); // full event plane res.
		double mRes2Delta = resEventPlane(::sqrt(2.) * chiSub2Delta);
		res2error = res2SubErr * fabs((double)res2 - mRes2Delta) 
		  / deltaRes2Sub;
	      }
	    } else { // neg. corr.
	      res2 = 0.;
	      res2error = 0.;
	    }
	    // now put everything together with first order event plane
	    mRes[k][j]    = ::sqrt(cosPair[k][0]*res2);
	    mResErr[k][j] = 1./(2.*mRes[k][j])*::sqrt(res2*res2*cosPairErr[k][0]*cosPairErr[k][0]
	    + cosPair[k][0]*cosPair[k][0]*res2error*res2error); // Gaussian error propagation
	    if (!pFlowEvent->EtaSubs()) {
	      // correct to full event plane resolution
	      // 1st order res for k=1 is small and linear
	      mRes[k][j] *= ::sqrt(2.);
	      mResErr[k][j] *= ::sqrt(2.);
	    }
	  } else if (pFlowEvent->EtaSubs() || pFlowEvent->RanSubs()) { // sub res only
	    resSub = ::sqrt(cosPair[k][j]);
	    resSubErr = cosPairErr[k][j] / (2. * resSub);
	    mRes[k][j]    = resSub;
	    mResErr[k][j] = resSubErr;
	  } else if (order==4. || order==6.|| order==8.) { // 2nd harmonic event plane
	    double deltaResSub = 0.005;  // differential for the error propagation
	    double resSub = ::sqrt(cosPair[k][1]);
	    double resSubErr = cosPairErr[k][1] / (2. * resSub);
	    double chiSub = chi(resSub);
	    double chiSubDelta = chi(resSub + deltaResSub);
	    double mResDelta;
	    if (order==4.) {
	      mRes[k][j] = resEventPlaneK2(::sqrt(2.) * chiSub); // full event plane res.
	      mResDelta = resEventPlaneK2(::sqrt(2.) * chiSubDelta);
	    } else if (order==6.) {
	      mRes[k][j] = resEventPlaneK3(::sqrt(2.) * chiSub); // full event plane res.
	      mResDelta = resEventPlaneK3(::sqrt(2.) * chiSubDelta);
	    } else {
	      mRes[k][j] = resEventPlaneK4(::sqrt(2.) * chiSub); // full event plane res.
	      mResDelta = resEventPlaneK4(::sqrt(2.) * chiSubDelta);
	    }
	    mResErr[k][j] = resSubErr * fabs((double)mRes[k][j] - mResDelta) / deltaResSub;
	  } else { //normal case
	    if (cosPair[k][j] > 0.92) { // resolution saturates
	      mRes[k][j]    = 0.99;
	      mResErr[k][j] = 0.007;
	    } else {
	      double deltaResSub = 0.005;  // differential for the error propagation
	      double resSub = ::sqrt(cosPair[k][j]);
	      double resSubErr = cosPairErr[k][j] / (2. * resSub);
	      double chiSub = chi(resSub);
	      double chiSubDelta = chi(resSub + deltaResSub);
	      mRes[k][j] = resEventPlane(::sqrt(2.) * chiSub); // full event plane res.
	      double mResDelta = resEventPlane(::sqrt(2.) * chiSubDelta);
	      mResErr[k][j] = resSubErr * fabs((double)mRes[k][j] - mResDelta) 
		/ deltaResSub;
	    }
	  }
	} else {      // subevent correlation must be positive
	  mRes[k][j]    = 0.;
	  mResErr[k][j] = 0.;
	}
      }//else :standard way if(!pFlowEvent->UseZDCSMD())
      histFull[k].mHistRes->SetBinContent(j+1, mRes[k][j]);
      histFull[k].mHistRes->SetBinError(j+1, mResErr[k][j]);

	// Create the v 2D histogram
      histTitle = new TString("Flow_v2D_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_v2D = 
	histFull[k].histFullHar[j].mHist_vObs2D->ProjectionXY(histTitle->Data());
      histFull[k].histFullHar[j].mHist_v2D->SetTitle(histTitle->Data());
      histFull[k].histFullHar[j].mHist_v2D->SetXTitle((char*)xLabel.Data());
      histFull[k].histFullHar[j].mHist_v2D->SetYTitle("Pt (GeV/c)");
      histFull[k].histFullHar[j].mHist_v2D->SetZTitle("v (%)");
      delete histTitle;
      AddHist(histFull[k].histFullHar[j].mHist_v2D);

      // Create the 1D v histograms
      histTitle = new TString("Flow_vEta_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_vEta = 
	histFull[k].histFullHar[j].mHist_vObsEta->ProjectionX(histTitle->Data());
      histFull[k].histFullHar[j].mHist_vEta->SetTitle(histTitle->Data());
      histFull[k].histFullHar[j].mHist_vEta->SetXTitle((char*)xLabel.Data());
      histFull[k].histFullHar[j].mHist_vEta->SetYTitle("v (%)");
      delete histTitle;
      AddHist(histFull[k].histFullHar[j].mHist_vEta);

      TString* histTitle = new TString("Flow_vPt_Sel");
      *histTitle += k+1;
      histTitle->Append("_Har");
      *histTitle += j+1;
      histFull[k].histFullHar[j].mHist_vPt = 
	histFull[k].histFullHar[j].mHist_vObsPt->ProjectionX(histTitle->Data());
      histFull[k].histFullHar[j].mHist_vPt->SetTitle(histTitle->Data());
      histFull[k].histFullHar[j].mHist_vPt->SetXTitle("Pt (GeV/c)");
      histFull[k].histFullHar[j].mHist_vPt->SetYTitle("v (%)");
      delete histTitle;
      AddHist(histFull[k].histFullHar[j].mHist_vPt);

      // Calulate v = vObs / Resolution
      if (mRes[k][j]) {
	cout << endl << "##### Resolution of the " << j+1 << "th harmonic = " << 
	  mRes[k][j] << " +/- " << mResErr[k][j] << endl;
	// The systematic error of the resolution is not folded in.
	histFull[k].histFullHar[j].mHist_v2D-> Scale(1. / mRes[k][j]);
	histFull[k].histFullHar[j].mHist_vEta->Scale(1. / mRes[k][j]);
	histFull[k].histFullHar[j].mHist_vPt ->Scale(1. / mRes[k][j]);
	content = histFull[k].mHist_v->GetBinContent(j+1);
	content /=  mRes[k][j];
	histFull[k].mHist_v->SetBinContent(j+1, content);
	// The systematic error of the resolution is folded in.
	error = histFull[k].mHist_v->GetBinError(j+1);
	error /= mRes[k][j];
	totalError = fabs(content) * ::sqrt((error/content)*(error/content) +
	       (mResErr[k][j]/mRes[k][j])*(mResErr[k][j]/mRes[k][j]));
	histFull[k].mHist_v->SetBinError(j+1, totalError);

	// Calculate non-flatness correction
	TF1* funcCosSin = new TF1("funcCosSin",
		   "[3]*(1.+[0]*2./100.*cos([2]*x)+[1]*2./100.*sin([2]*x))", 0., twopi);
	//funcCosSin->SetParNames("100*cos", "100*sin", "har");
	funcCosSin->SetParameters(0, 0, j+1); // initial values
	funcCosSin->SetParLimits(2, 1, 1); // har is fixed
	histFull[k].histFullHar[j].mHistPhiLab->Fit("funcCosSin","Q,N");
	Double_t cosParLab = funcCosSin->GetParameter(0);
	Double_t sinParLab = funcCosSin->GetParameter(1);
	histFull[k].histFullHar[j].mHistPsi->Fit("funcCosSin","Q,N");
	Double_t cosParEP = funcCosSin->GetParameter(0);
	Double_t sinParEP = funcCosSin->GetParameter(1);
	cout << "100*cosLab = " << cosParLab << ", 100*sinLab = " << sinParLab << endl;
	cout << "100*cosEP = " << cosParEP << ", 100*sinEP = " << sinParEP << endl;
	delete funcCosSin;
	float nonflat = (cosParEP*cosParLab + sinParEP*sinParLab)/mRes[k][j]/100.;
	cout << "nonflat = " << nonflat << endl;

	cout << "##### v" << j+1 << "= (" << content << " +/- " << error << ")%" << endl;
	cout << "##### v" << j+1 << "= (" << content - nonflat << " (with nonflat corr.) +/- "
	     << totalError << " (with syst.) )%" << endl;
	histFull[k].mHist_v->SetBinContent(j+1, content - nonflat);


      } else {
	cout << "##### Resolution of the " << j+1 << "th harmonic was zero."
	     << endl;
	histFull[k].histFullHar[j].mHist_v2D-> Reset();
	histFull[k].histFullHar[j].mHist_vEta->Reset();
	histFull[k].histFullHar[j].mHist_vPt ->Reset();
	histFull[k].mHist_v->SetBinContent(j+1, 0.);
	histFull[k].mHist_v->SetBinError(j+1, 0.);
      }//v

      // Calculate PhiWgt
//      if (pFlowMaker->PhiWgtCalc()) {
      if (1) {
	if (j < 2) {
	  double meanFarEast = histFull[k].histTwoHar[j].mHistPhiFarEast->Integral() 
	    / (double)Flow::nPhiBins;
	  double meanEast = histFull[k].histTwoHar[j].mHistPhiEast->Integral() 
	    / (double)Flow::nPhiBins;
	  double meanWest = histFull[k].histTwoHar[j].mHistPhiWest->Integral() 
	    / (double)Flow::nPhiBins;
	  double meanFarWest = histFull[k].histTwoHar[j].mHistPhiFarWest->Integral() 
	    / (double)Flow::nPhiBins;
	  double meanFtpcFarEast = histFull[k].histTwoHar[j].mHistPhiFtpcFarEast->Integral() 
	    / (double)Flow::nPhiBinsFtpc;
	  double meanFtpcEast = histFull[k].histTwoHar[j].mHistPhiFtpcEast->Integral() 
	    / (double)Flow::nPhiBinsFtpc;
	  double meanFtpcWest = histFull[k].histTwoHar[j].mHistPhiFtpcWest->Integral() 
	    / (double)Flow::nPhiBinsFtpc;
	  double meanFtpcFarWest = histFull[k].histTwoHar[j].mHistPhiFtpcFarWest->Integral() 
	    / (double)Flow::nPhiBinsFtpc;
	  
	  // Tpc
	  for (int i = 0; i < Flow::nPhiBins; i++) {
	    histFull[k].histTwoHar[j].mHistPhiWgtFarEast->SetBinContent(i+1,meanFarEast);
	    histFull[k].histTwoHar[j].mHistPhiWgtFarEast->SetBinError(i+1, 0.);
	    histFull[k].histTwoHar[j].mHistPhiWgtEast   ->SetBinContent(i+1, meanEast);
	    histFull[k].histTwoHar[j].mHistPhiWgtEast   ->SetBinError(i+1, 0.);
	    histFull[k].histTwoHar[j].mHistPhiWgtWest   ->SetBinContent(i+1, meanWest);
	    histFull[k].histTwoHar[j].mHistPhiWgtWest   ->SetBinError(i+1, 0.);
	    histFull[k].histTwoHar[j].mHistPhiWgtFarWest->SetBinContent(i+1,meanFarWest);
	    histFull[k].histTwoHar[j].mHistPhiWgtFarWest->SetBinError(i+1, 0.);
	  }
	  
	  // Ftpc
	  for (int i = 0; i < Flow::nPhiBinsFtpc; i++) {
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast->SetBinContent(i+1,meanFtpcFarEast);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast->SetBinError(i+1, 0.);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast   ->SetBinContent(i+1,meanFtpcEast);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast   ->SetBinError(i+1, 0.);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest   ->SetBinContent(i+1,meanFtpcWest);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest   ->SetBinError(i+1, 0.);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest->SetBinContent(i+1,meanFtpcFarWest);
	    histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest->SetBinError(i+1, 0.);
	  }
	  
	  // Tpc
	  histFull[k].histTwoHar[j].mHistPhiWgtFarEast->
	    Divide(histFull[k].histTwoHar[j].mHistPhiFarEast);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtFarEast);
	  histFull[k].histTwoHar[j].mHistPhiWgtEast->
	    Divide(histFull[k].histTwoHar[j].mHistPhiEast);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtEast);
	  histFull[k].histTwoHar[j].mHistPhiWgtWest->
	    Divide(histFull[k].histTwoHar[j].mHistPhiWest);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtWest);
	  histFull[k].histTwoHar[j].mHistPhiWgtFarWest->
	    Divide(histFull[k].histTwoHar[j].mHistPhiFarWest);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtFarWest);
	  
	  // Ftpc
	  histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast->
	    Divide(histFull[k].histTwoHar[j].mHistPhiFtpcFarEast);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarEast);
	  histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast->
	    Divide(histFull[k].histTwoHar[j].mHistPhiFtpcEast);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtFtpcEast);
	  histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest->
	    Divide(histFull[k].histTwoHar[j].mHistPhiFtpcWest);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtFtpcWest);
	  histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest->
	    Divide(histFull[k].histTwoHar[j].mHistPhiFtpcFarWest);
	  phiWgtHistNames->AddLast(histFull[k].histTwoHar[j].mHistPhiWgtFtpcFarWest);
	}
      }//phiWgt
    }
  }
  phiWgtHistNames->AddLast(mHistZDCSMDPsiWgtEast);
  phiWgtHistNames->AddLast(mHistZDCSMDPsiWgtWest);
  TFile* pPhiWgtFile = new TFile("flowPhiWgt.hist.root", "READ");
  if (pPhiWgtFile->IsOpen())
    { phiWgtHistNames->AddLast(mHistZDCSMDPsiWgtFull); }

  // Write all histograms
  TFile histFile("flow.hist.root", "RECREATE");
  //GetHistList()->ls();
  GetHistList()->Write();
  histFile.Close();
  
  // Write PhiWgt histograms preceded by documenting text
  if (1) {
//  if (pFlowMaker->PhiWgtCalc()) {
    TFile phiWgtNewFile("flowPhiWgtNew.hist.root", "RECREATE");
    TText* textInfo = 0;
    if (pFlowEvent->FirstLastPoints()) {
      char chInfo[400];
      sprintf(chInfo, "%s%d%s%d%s", " pt weight= ", pFlowEvent->PtWgt(),
	      ", eta weight= ", pFlowEvent->EtaWgt(), "\n");
      textInfo = new TText(0,0,chInfo);
      textInfo->Write("info");
    }
    phiWgtNewFile.cd();
    phiWgtHistNames->Write();
    phiWgtNewFile.Close();
    if (pFlowEvent->FirstLastPoints()) delete textInfo;
  }
  delete phiWgtHistNames;

  // Write reCent values
  if (mCalcReCentPars) {
    TFile fileReCent("flow.reCentAnaNew.root", "RECREATE");
    fileReCent.cd();
    savedHistReCentNames->Write();
    fileReCent.Close();
  }
  delete savedHistReCentNames;
  
  delete pFlowSelect;

  //# Added by jcampbell
  if(mPsiShiftFile) { 
    mPsiShiftFile->cd();
    for (int k = 0; k < Flow::nSels; k++) {
      for (int j = 0; j < Flow::nHars; j++) {
        for (int n = 0; n < Flow::nSubs+1; n++) {
	    	if(mFillPsiShift) {
			  mPsiShiftSin[k][j][n]->Write();
			  mPsiShiftCos[k][j][n]->Write();
			}
			delete mPsiShiftSin[k][j][n];
			delete mPsiShiftCos[k][j][n];
        } // nSubs
      } // nHars
	} // nSels
	mPsiShiftFile->Close();
	delete mPsiShiftFile;
  } // If mPsiShiftFile

  return StMaker::Finish();
}

//-----------------------------------------------------------------------

void StFlowAnalysisMaker::SetHistoRanges(Bool_t ftpc_included) {

    if (ftpc_included) {
	  mEtaMin = Flow::etaMin;
	  mEtaMax = Flow::etaMax;
	mNEtaBins = Flow::nEtaBins;
    }
    else {
	  mEtaMin = Flow::etaMinTpcOnly;
 	  mEtaMax = Flow::etaMaxTpcOnly;
	mNEtaBins = Flow::nEtaBinsTpcOnly;
    }

    return;
}

//------------------------------------------------------------------------

void StFlowAnalysisMaker::SetPtRange_for_vEta(Float_t lo, Float_t hi) {

  // Sets the pt range for the v(eta) histograms.

  mPtRange_for_vEta[0] = lo;
  mPtRange_for_vEta[1] = hi;

  return;
}

//------------------------------------------------------------------------

void StFlowAnalysisMaker::SetEtaRange_for_vPt(Float_t lo, Float_t hi) {
  
  // Sets the |eta| range for the v(pt) histograms.

  mEtaRange_for_vPt[0] = lo;
  mEtaRange_for_vPt[1] = hi;

  return;
}

//------------------------------------------------------------------------

void StFlowAnalysisMaker::SetV1Ep1Ep2(Bool_t v1Ep1Ep2) {
  
  // Switches the v_1{EP1,EP2} calculation on/off.

  mV1Ep1Ep2 = v1Ep1Ep2;

  return;
}


////////////////////////////////////////////////////////////////////////////
//
// $Log: StFlowAnalysisMaker.cxx,v $
// Revision 1.103  2011/07/25 15:54:42  posk
// Added correction for non-flatness of event plane.
//
// Revision 1.102  2011/03/10 18:56:20  posk
// Added histogram for laboratory azimuthal distribution of particles.
//
// Revision 1.101  2010/09/30 19:28:09  posk
// Instead of reversing the weight for negative pseudrapidity for odd harmonics,
// it is now done only for the first harmonic.
// Recentering is now done for all harmonics.
//
// Revision 1.100  2010/02/15 12:01:58  canson
// Changed mHistCTBvsZDC2D from filling with ZDC_e + ZDC_e to filling with ZDC_e + ZDC_w
//
// Revision 1.99  2009/11/24 19:29:11  posk
// Added reCenter to remove acceptance correlations as an option instead of phiWgt.
//
// Revision 1.98  2007/07/13 22:18:29  posk
// Method chi() revised by Wang Gang: "With this modification,
// it will give good results not only for realistic resolutions, but also
// for res=1 and res=0."
//
// Revision 1.97  2007/02/06 19:00:39  posk
// In Lee Yang Zeros method, introduced recentering of Q vector.
// Reactivated eta symmetry cut.
//
// Revision 1.96  2006/07/10 21:03:48  posk
// For profile histograms of v, changed the limits to -1000, 1000.
//
// Revision 1.95  2006/02/22 19:36:21  posk
// Minor updates.
//
// Revision 1.94  2005/08/26 19:00:15  posk
// plot style back to bold
//
// Revision 1.93  2005/08/05 20:13:35  posk
// Improved first guess for qDist fit.
//
// Revision 1.92  2005/02/11 23:17:14  posk
// Fixed trigger histogram.
//
// Revision 1.91  2005/02/08 22:37:53  posk
// Fixed trigger histogram for year=4.
//
// Revision 1.90  2004/12/20 19:41:24  aihong
// crashes when run without ZDCSMD. bug fixed
//
// Revision 1.89  2004/12/17 22:33:35  aihong
// add in full Psi weight for ZDC SMD and fix a few bugs, done by Gang
//
// Revision 1.88  2004/12/09 23:47:05  posk
// Minor changes in code formatting.
// Added hist for TPC primary dca to AnalysisMaker.
//
// Revision 1.87  2004/12/07 23:10:19  posk
// Only odd and even phiWgt hists. If the old phiWgt file contains more than
// two harmonics, only the first two are read. Now writes only the first two.
//
// Revision 1.86  2004/08/24 20:22:36  oldi
// Minor modifications to avoid compiler warnings.
//
// Revision 1.85  2004/08/18 00:18:59  oldi
// Several changes were necessary to comply with latest changes of MuDsts and StEvent:
//
// nHits, nFitPoints, nMaxPoints
// -----------------------------
// From now on
//  - the fit points used in StFlowMaker are the fit points within the TPC xor FTPC (vertex excluded).
//  - the max. possible points used in StFlowMAker are the max. possible points within the TPC xor FTPC (vertex excluded).
//  - the number of points (nHits; not used for analyses so far) are the total number of points on a track, i. e.
//    TPC + SVT + SSD + FTPCeast + FTPCwest [reading from HBT event gives a warning, but it seems like nobody uses it anyhow].
// - The fit/max plot (used to be (fit-1)/max) was updated accordingly.
// - The default cuts for fit points were changed (only for the FTPC, since TPC doesn't set default cuts).
// - All these changes are backward compatible, as long as you change your cuts for the fit points by 1 (the vertex used to
//   be included and is not included anymore). In other words, your results won't depend on old or new MuDst, StEvent,
//   PicoDsts as long as you use the new flow software (together with the latest MuDst and StEvent software version).
// - For backward compatibility reasons the number of fit points which is written out to the flowpicoevent.root file
//   includes the vertex. It is subtracted internally while reading back the pico files. This is completely hidden from the
//   user.
//
// zFirstPoint
// -----------
// The positions of the first point of tracks which have points in the TPC can lie outside of the TPC (the tracks can start in
// the SVT or SSD now). In this case, the first point of the track is obtained by extrapolating the track helix to the inner
// radius of the TPC.
//
// Revision 1.84  2004/05/31 20:09:22  oldi
// PicoDst format changed (Version 7) to hold ZDC SMD information.
// Trigger cut modified to comply with TriggerCollections.
// Centrality definition for 62 GeV data introduced.
// Minor bug fixes.
//
// Revision 1.83  2004/05/05 21:13:47  aihong
// Gang's code for ZDC-SMD added
//
// Revision 1.82  2004/03/11 18:00:03  posk
// Added Random Subs analysis method.
//
// Revision 1.81  2003/12/12 02:34:40  oldi
// Removal of some major bugs in the v1{EP1,EP2} method.
//
// Revision 1.80  2003/12/09 01:40:15  oldi
// Removed 'inline' of some functions to cope with new compiler.
//
// Revision 1.79  2003/11/14 20:00:40  oldi
// Implementation of v1{EP1,EP2}. This method is set to be the default for v1 now!
// Minor code clean-ups.
//
// Revision 1.78  2003/09/02 17:58:10  perev
// gcc 3.2 updates + WarnOff
//
// Revision 1.77  2003/08/26 21:10:10  posk
// Calculates v8 if nHars=8.
//
// Revision 1.76  2003/08/06 20:54:06  oldi
// Introduction of possibility to exclude pt ranges for v(eta) and eta regions
// for v(pt) histograms. Default behavior stays the same (all available tracks
// are included in v(pt) and v(eta)).
//
// Revision 1.75  2003/07/30 22:08:25  oldi
// Several code fixes for EtaSym plots introduced (esp. the acceptance correction
// is done now for 200 GeV data and for the FTPCs as well).
// PtWgtSaturation parameter introduced.
//
// Revision 1.74  2003/07/07 21:58:16  posk
// Made units of momentum GeV/c instead of GeV.
//
// Revision 1.73  2003/06/27 21:25:41  posk
// v4 and v6 are with repect to the 2nd harmonic event plane.
//
// Revision 1.72  2003/05/16 20:44:46  posk
// First commit of StFlowPhiWgtMaker
//
// Revision 1.71  2003/05/06 21:33:04  posk
// Removed some histograms.
//
// Revision 1.70  2003/05/06 18:38:05  posk
// Removed StFlowTagMaker.
//
// Revision 1.69  2003/01/16 16:02:27  posk
// Some plotting changes.
//
// Revision 1.68  2003/01/10 16:40:16  oldi
// Several changes to comply with FTPC tracks:
// - Switch to include/exclude FTPC tracks introduced.
//   The same switch changes the range of the eta histograms.
// - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
// - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
//   West, FarWest (depending on vertex.z()).
// - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
// - Cut to exclude mu-events with no primary vertex introduced.
//   (This is possible for UPC events and FTPC tracks.)
// - Global DCA cut for FTPC tracks added.
// - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
// - Charge cut for FTPC tracks added.
//
// Revision 1.67  2003/01/08 19:28:07  posk
// PhiWgt hists sorted on sign of z of first and last points.
//
// Revision 1.66  2002/10/28 19:45:52  posk
// Eliminate events with Psi=0.
//
// Revision 1.65  2002/06/15 23:03:56  posk
// Changed Fit/Max histogram to (Fit - 1)/Max.
//
// Revision 1.64  2002/05/23 18:57:08  posk
// changed label on MultHist histogram
//
// Revision 1.63  2002/05/21 18:42:15  posk
// Kirill's correction to minBias.C for bins with one count.
//
// Revision 1.62  2002/02/13 22:31:27  posk
// Pt Weight now also weights Phi Weight. Added Eta Weught, default=FALSE.
//
// Revision 1.61  2002/01/31 01:09:10  posk
// *** empty log message ***
//
// Revision 1.60  2002/01/14 23:42:21  posk
// Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
//
// Revision 1.59  2001/12/18 19:27:06  posk
// "proton" and "antiproton" replaced by "pr+" and "pr-".
//
// Revision 1.58  2001/12/11 22:03:41  posk
// Four sets of phiWgt histograms.
// StFlowMaker StFlowEvent::PhiWeight() changes.
// Cumulant histogram names changed.
//
// Revision 1.57  2001/11/13 22:47:17  posk
// Documentation updated. Fit to q function moved to macro.
//
// Revision 1.56  2001/11/10 01:09:05  posk
// Moved some constants into StFlowConstants.
//
// Revision 1.55  2001/11/09 21:14:33  posk
// Switched from CERNLIB to TMath. Using global dca instead of dca.
//
// Revision 1.54  2001/08/02 20:42:01  snelling
// removed hist title conflict
//
// Revision 1.53  2001/08/02 17:41:49  snelling
// Added trigger histogram
//
// Revision 1.52  2001/05/22 20:10:55  posk
// Changed dEdx graphs.
//
// Revision 1.51  2001/04/25 17:43:24  perev
// HPcorrs
//
// Revision 1.50  2001/04/03 17:46:06  oldi
// Bug fix that excluded FTPC tracks from the determination of the reaction plane.
//
// Revision 1.49  2000/12/12 15:01:10  posk
// Put log comments at end of file.
//
// Revision 1.48  2000/12/10 02:02:01  oldi
// A new member (StTrackTopologyMap mTopology) was added to StFlowPicoTrack.
// The evaluation of either a track originates from the FTPC or not is
// unambiguous now. The evaluation itself is easily extendible for other
// detectors (e.g. SVT+TPC). Old flowpicoevent.root files are treated as if
// they contain TPC tracks only (backward compatibility).
//
// Revision 1.47  2000/12/08 17:04:09  oldi
// Phi weights for both FTPCs included.
//
// Revision 1.43  2000/09/22 22:01:38  posk
// Doubly integrated v now contains resolution error.
//
// Revision 1.42  2000/09/16 22:23:04  snelling
// Auto magically switch to rapidity when identified particles are used
//
// Revision 1.41  2000/09/15 22:52:53  posk
// Added Pt weighting for event plane calculation.
//
// Revision 1.40  2000/09/12 01:31:00  snelling
// Added pid histograms for e- e+ and dbar
//
// Revision 1.39  2000/09/07 17:02:45  snelling
// Made the hist file standard root compatible
//
// Revision 1.38  2000/09/05 16:12:12  snelling
// Added the new particles to the pid histogram
//
// Revision 1.37  2000/08/31 18:50:29  posk
// Added plotCen.C to plot from a series of files with different centralities.
//
// Revision 1.36  2000/08/12 20:20:13  posk
// More centrality bins.
//
// Revision 1.35  2000/08/09 21:38:59  snelling
// Added monitor histograms
//
// Revision 1.34  2000/08/01 21:51:18  posk
// Added doubly integrated v.
//
// Revision 1.33  2000/07/12 17:49:37  posk
// Changed EtaSym plots.
//
// Revision 1.32  2000/06/30 14:51:18  posk
// Using MessageMgr. Added graph for Eta Symmetry vs. Vertex Z.
//
// Revision 1.31  2000/06/01 18:29:56  posk
// When resolution=0 reset histograms.
//
// Revision 1.30  2000/05/26 21:25:20  posk
// Use TProfile2D class and profile projection methods.
// Correction needed for >2 subevents.
//
// Revision 1.27  2000/04/13 22:34:13  posk
// Resolution correction is now made.
//
// Revision 1.26  2000/03/28 23:25:36  posk
// Allow multiple instances.
//
// Revision 1.25  2000/03/21 00:24:43  posk
// Added GetCVS and changed some plot names.
//
// Revision 1.24  2000/03/15 23:32:03  posk
// Added StFlowSelection.
//
// Revision 1.23  2000/03/02 22:55:32  posk
// Changed header file extensions from .hh to .h .
//
// Revision 1.22  2000/02/29 21:55:12  posk
// Removed static const int& statements.
//
// Revision 1.21  2000/02/18 23:44:52  posk
// Added PID and centrality.
//
// Revision 1.20  2000/02/10 01:47:30  snelling
// Make changes for HP compiler
//
// Revision 1.19  2000/02/04 16:26:41  posk
// Added correct calculation of event plane resolution for large flow.
//
// Revision 1.15  2000/01/14 01:35:52  snelling
// changed include path ../FlowMaker/ to FlowMaker/
//
// Revision 1.14  2000/01/14 01:13:34  snelling
// modified spt (sum pt) to mpt (mean pt) because FlowTag changed
//
// Revision 1.11  1999/12/04 00:15:39  posk
// Works with StFlowEvent which works with the new StEvent
//
// Revision 1.10  1999/11/24 18:14:05  posk
// Now reads event quantities with StFlowEvent methods
//
// Revision 1.8  1999/11/05 00:02:02  posk
// Changed the flow vector, Q, to a TVector2.
//
// Revision 1.7  1999/10/05 16:54:07  posk
// Added getPhiWeight method for making the event plane isotropic.
//
// Revision 1.6  1999/09/24 01:23:06  fisyak
// Reduced Include Path
//
// Revision 1.4  1999/09/03 01:05:59  fisyak
// replace iostream/stdlib by Stiostream.h/stdlib.h
//
// Revision 1.3  1999/08/24 18:02:37  posk
// Calculates event plane resolution.
// Added macros for plotting histograms.
//
// Revision 1.1.1.1  1999/08/09 19:50:37  posk
//
// Revision 1.0  1999/08/02 
//  
////////////////////////////////////////////////////////////////////////////