//////////////////////////////////////////////////////////////////////////// // // $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 #include #include #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(""); // 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(" (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(""); 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(""); 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(""); 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(""); 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(""); // 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 // ////////////////////////////////////////////////////////////////////////////