my_star/StRoot/StFlowAnalysisMaker/plot.C

///////////////////////////////////////////////////////////////////////////////
//
// $Id: plot.C,v 1.70 2011/07/25 15:54:53 posk Exp $
//
// Author:       Art Poskanzer, LBNL, Aug 1999
//               FTPC added by Markus Oldenburg, MPI, Dec 2000
// Description:  Macro to plot histograms made by StFlowAnalysisMaker,
//                 StFlowLeeYangZerosMaker, and StFlowScalarProdMaker.
//               If selN = 0 plot all selections and harmonics.
//               First time type .x plot.C() to see the menu.
//               Run Number appended to "ana" is entered in the bottom, left box.
//               Hist file is anaXX.root where XX is the number.
//               Default hist file is flow.hist.root .
//               After the first execution, just type plot(N) .
//               A negative N plots all pages starting with page N.
//               Place a symbolic link to this file in StRoot/macros/analysis .
//
///////////////////////////////////////////////////////////////////////////////

#include <math.h> 
#include "TMath.h" 
#include <iostream.h>
#include <iomanip.h>


const    Int_t nHars    = 4; // 4 
const    Int_t nSels    = 2;
const    Int_t nSubs    = 2;
Int_t    runNumber      = 0;
char     runName[6];
char     fileNumber[4]  = "x";
char     fileName[30];
TFile*   histFile;
char     tmp[10];
char     runNo[10];
TCanvas* can;

TCanvas* plot(Int_t pageNumber=0, Int_t selN=0, Int_t harN=0){
  gInterpreter->ProcessLine(".O0");

  Bool_t includeFtpc = kTRUE;
  Bool_t reCent = kFALSE;

  bool multiGraph  = kFALSE;                            // set flags
  bool singleGraph = kFALSE;
  if (selN == 0) multiGraph = kTRUE;

  TCanvas* cOld = (TCanvas*)gROOT->GetListOfCanvases(); // delete old canvas
  if (cOld) cOld->Delete();
    
  //gROOT->SetStyle("Pub");                              // set style
  gROOT->SetStyle("Bold");                              // set style
  gROOT->ForceStyle();
  gStyle->SetOptStat(10);
  gStyle->SetPalette(1);

  // names of histograms made by StFlowAnalysisMaker
  // also projections of some of these histograms
  const char* baseName1[] = { // with FTPCs
    "Flow_Res_Sel",
    "Flow_Trigger",
    "Flow_VertexZ",
    "Flow_VertexXY2D",
    "Flow_EtaSymVerZ2D_Tpc",
    "Flow_EtaSym_Tpc",
    "Flow_EtaSymVerZ2D_Ftpc",
    "Flow_EtaSym_Ftpc",
    "Flow_CTBvsZDC2D",
    "Flow_Cent",
    "Flow_Mult",
    "Flow_OrigMult",
    "Flow_MultOverOrig",
    "Flow_MultEta",
    "Flow_MultPart",
    "Flow_Charge_Ftpc",
    "Flow_DcaGlobal_Tpc",
    "Flow_DcaGlobal_Ftpc",
    "Flow_Dca_Tpc",
    "Flow_Dca_Ftpc",
    "Flow_Chi2_Tpc",
    "Flow_Chi2_Ftpc",
    "Flow_FitPts_Tpc",
    "Flow_MaxPts_Tpc",
    "Flow_FitOverMax_Tpc",
    "Flow_FitPts_Ftpc",
    "Flow_MaxPts_Ftpc",
    "Flow_FitOverMax_Ftpc",
    //"Flow_EtaPtPhi3D",
    "Flow_EtaPtPhi2D.PhiEta",
    "Flow_EtaPtPhi2D.PhiPt",
    "Flow_YieldAll2D",
    "Flow_YieldAll.Eta",
    "Flow_YieldAll.Pt",
    "Flow_YieldPart2D",
    "Flow_YieldPart.Eta",
    "Flow_YieldPart.Pt",
    //"Flow_Bin_Eta",
    //"Flow_Bin_Pt",
    "Flow_MeanDedxPos2D",
    "Flow_MeanDedxNeg2D",
//     "Flow_PidPiPlusPart",
//     "Flow_PidPiMinusPart",
//     "Flow_PidProtonPart",
//     "Flow_PidAntiProtonPart",
//     "Flow_PidKplusPart",
//     "Flow_PidKminusPart",
//     "Flow_PidDeuteronPart",
//     "Flow_PidAntiDeuteronPart",
//     "Flow_PidElectronPart",
//     "Flow_PidPositronPart",
    "Flow_PidMult",
    "Flow_Phi_FarEast_Sel",                      // first multi graph hist
    "Flow_Phi_Flat_FarEast_Sel",
    "Flow_Phi_East_Sel",
    "Flow_Phi_Flat_East_Sel",
    "Flow_Phi_West_Sel",
    "Flow_Phi_Flat_West_Sel",
    "Flow_Phi_FarWest_Sel",
    "Flow_Phi_Flat_FarWest_Sel",
    "Flow_Phi_FtpcFarEast_Sel",
    "Flow_Phi_Flat_FtpcFarEast_Sel",
    "Flow_Phi_FtpcEast_Sel",
    "Flow_Phi_Flat_FtpcEast_Sel",
    "Flow_Phi_FtpcWest_Sel",
    "Flow_Phi_Flat_FtpcWest_Sel",
    "Flow_Phi_FtpcFarWest_Sel",
    "Flow_Phi_Flat_FtpcFarWest_Sel",
    "Flow_Mul_Sel",
    "Flow_Yield2D_Sel",
    "Flow_Yield.Eta_Sel",
    "Flow_Yield.Pt_Sel",
    "Flow_Psi_Subs",
    "Flow_Psi_Sel",
    "Flow_PhiLab_Sel",
    "Flow_Psi_Diff_Sel",
    "Flow_Psi_Sub_Corr_Sel",
    "Flow_Psi_Sub_Corr_Diff_Sel",
    "Flow_Phi_Corr_Sel",
    "Flow_vObs2D_Sel",
    "Flow_vObsEta_Sel",
    "Flow_vObsPt_Sel",
    "Flow_v2D_Sel",
    "Flow_vEta_Sel",
    "Flow_vPt_Sel",
    "Flow_q_Sel",
    "FlowLYZ_vEta_Sel",
    "FlowLYZ_vPt_Sel",
    "Flow_vObs2D_ScalarProd_Sel",
    "Flow_v2D_ScalarProd_Sel",
    "Flow_vEta_ScalarProd_Sel",
    "Flow_vPt_ScalarProd_Sel"
  };

  const char* baseName2[] = { // without FTPCs
    "Flow_Res_Sel",
    "Flow_Trigger",
    "Flow_VertexZ",
    "Flow_VertexXY2D",
    "Flow_EtaSymVerZ2D_Tpc",
    "Flow_EtaSym_Tpc",
    "Flow_CTBvsZDC2D",
    "Flow_Cent",
    "Flow_Mult",
    "Flow_OrigMult",
    "Flow_MultOverOrig",
    "Flow_MultEta",
    "Flow_MultPart",
    "Flow_Dca_Tpc",
    "Flow_DcaGlobal_Tpc",
    "Flow_Chi2_Tpc",
    "Flow_FitPts_Tpc",
    "Flow_MaxPts_Tpc",
    "Flow_FitOverMax_Tpc",
    //"Flow_EtaPtPhi3D",
    "Flow_EtaPtPhi2D.PhiEta",
    "Flow_EtaPtPhi2D.PhiPt",
    "Flow_YieldAll2D",
    "Flow_YieldAll.Eta",
    "Flow_YieldAll.Pt",
    "Flow_YieldPart2D",
    "Flow_YieldPart.Eta",
    "Flow_YieldPart.Pt",
    //"Flow_Bin_Eta",
    //"Flow_Bin_Pt",
//     "Flow_MeanDedxPos2D",
//     "Flow_MeanDedxNeg2D",
//     "Flow_PidPiPlusPart",
//     "Flow_PidPiMinusPart",
//     "Flow_PidProtonPart",
//     "Flow_PidAntiProtonPart",
//     "Flow_PidKplusPart",
//     "Flow_PidKminusPart",
//     "Flow_PidDeuteronPart",
//     "Flow_PidAntiDeuteronPart",
//     "Flow_PidElectronPart",
//     "Flow_PidPositronPart",
    "Flow_PidMult",
    "Flow_Phi_FarEast_Sel",                      // first multi graph hist
    "Flow_Phi_Flat_FarEast_Sel",
    "Flow_Phi_East_Sel",
    "Flow_Phi_Flat_East_Sel",
    "Flow_Phi_West_Sel",
    "Flow_Phi_Flat_West_Sel",
    "Flow_Phi_FarWest_Sel",
    "Flow_Phi_Flat_FarWest_Sel",
    "Flow_Mul_Sel",
    "Flow_Yield2D_Sel",
    "Flow_Yield.Eta_Sel",
    "Flow_Yield.Pt_Sel",
    "Flow_Psi_Subs",
    "Flow_Psi_Sel",
    "Flow_PhiLab_Sel",
    "Flow_Psi_Diff_Sel",
    "Flow_Psi_Sub_Corr_Sel",
    "Flow_Psi_Sub_Corr_Diff_Sel",
    "Flow_Phi_Corr_Sel",
    "Flow_vObs2D_Sel",
    "Flow_vObsEta_Sel",
    "Flow_vObsPt_Sel",
    "Flow_v2D_Sel",
    "Flow_vEta_Sel",
    "Flow_vPt_Sel",
    "Flow_q_Sel",
    "FlowLYZ_vEta_Sel",
    "FlowLYZ_vPt_Sel",
    "Flow_vObs2D_ScalarProd_Sel",
    "Flow_v2D_ScalarProd_Sel",
    "Flow_vEta_ScalarProd_Sel",
    "Flow_vPt_ScalarProd_Sel"
  };

  const char* baseName3[] = { // for reCent
    "Flow_Res_Sel",
    "Flow_Trigger",
    "Flow_VertexZ",
    "Flow_VertexXY2D",
    "Flow_EtaSymVerZ2D_Tpc",
    "Flow_EtaSym_Tpc",
    "Flow_EtaSymVerZ2D_Ftpc",
    "Flow_EtaSym_Ftpc",
    "Flow_CTBvsZDC2D",
    "Flow_Cent",
    "Flow_Mult",
    "Flow_OrigMult",
    "Flow_MultOverOrig",
    "Flow_MultEta",
    "Flow_MultPart",
    "Flow_Charge_Ftpc",
    "Flow_DcaGlobal_Tpc",
    "Flow_DcaGlobal_Ftpc",
    "Flow_Dca_Tpc",
    "Flow_Dca_Ftpc",
    "Flow_Chi2_Tpc",
    "Flow_Chi2_Ftpc",
    "Flow_FitPts_Tpc",
    "Flow_MaxPts_Tpc",
    "Flow_FitOverMax_Tpc",
    "Flow_FitPts_Ftpc",
    "Flow_MaxPts_Ftpc",
    "Flow_FitOverMax_Ftpc",
    //"Flow_EtaPtPhi3D",
    "Flow_EtaPtPhi2D.PhiEta",
    "Flow_EtaPtPhi2D.PhiPt",
    "Flow_YieldAll2D",
    "Flow_YieldAll.Eta",
    "Flow_YieldAll.Pt",
    "Flow_YieldPart2D",
    "Flow_YieldPart.Eta",
    "Flow_YieldPart.Pt",
    //"Flow_Bin_Eta",
    //"Flow_Bin_Pt",
    "Flow_MeanDedxPos2D",
    "Flow_MeanDedxNeg2D",
//     "Flow_PidPiPlusPart",
//     "Flow_PidPiMinusPart",
//     "Flow_PidProtonPart",
//     "Flow_PidAntiProtonPart",
//     "Flow_PidKplusPart",
//     "Flow_PidKminusPart",
//     "Flow_PidDeuteronPart",
//     "Flow_PidAntiDeuteronPart",
//     "Flow_PidElectronPart",
//     "Flow_PidPositronPart",
    "Flow_PidMult",
    "Flow_Mul_Sel",                    // first multi graph hist
    "Flow_Yield2D_Sel",
    "Flow_Yield.Eta_Sel",
    "Flow_Yield.Pt_Sel",
    "Flow_Psi_Subs",
    "Flow_Psi_Sel",
    "Flow_PhiLab_Sel",
    "Flow_Psi_Diff_Sel",
    "Flow_Psi_Sub_Corr_Sel",
    "Flow_Psi_Sub_Corr_Diff_Sel",
    "Flow_Phi_Corr_Sel",
    "Flow_QXY2D_Sel",
    "Flow_QFTPCSubXY2D_Sel",
    "Flow_QTPCSubXY2D_Sel",
    "Flow_vObs2D_Sel",
    "Flow_vObsEta_Sel",
    "Flow_vObsPt_Sel",
    "Flow_v2D_Sel",
    "Flow_vEta_Sel",
    "Flow_vPt_Sel",
    "Flow_q_Sel",
    "FlowLYZ_vEta_Sel",
    "FlowLYZ_vPt_Sel",
    "Flow_vObs2D_ScalarProd_Sel",
    "Flow_v2D_ScalarProd_Sel",
    "Flow_vEta_ScalarProd_Sel",
    "Flow_vPt_ScalarProd_Sel"
  };

  const Int_t firstMultFtpc = 38; // these must be set by hand
  const Int_t firstMultTpc  = 27;
  Int_t nSingles;
  Int_t nName;
  if (reCent) {
    nName = sizeof(baseName3) / sizeof(char*);
    nSingles = firstMultFtpc + 1;
  } else if (includeFtpc) {
    nName = sizeof(baseName1) / sizeof(char*);
    nSingles = firstMultFtpc + 1;
  } else {
    nName = sizeof(baseName2) / sizeof(char*);
    nSingles = firstMultTpc + 1;
  }
  const Int_t nNames = nName;

  // construct arrays of base and short names
  char* baseName[nNames];
  char* shortName[nNames];
  for (int n = 0; n < nNames; n++) {
    baseName[n] = new char[35];
    shortName[n] = new char[35];
    float etaMax;
    if (reCent) {
      strcpy(baseName[n], baseName3[n]);
      etaMax  =   4.5;
    } else if (includeFtpc) {
      strcpy(baseName[n], baseName1[n]);
      etaMax  =   4.5;
    } else {
      strcpy(baseName[n], baseName2[n]);
      etaMax  =   1.5;
    }
    strcpy(shortName[n], baseName[n]);
    char* cp = strstr(shortName[n],"_Sel");
    if (cp) *cp = '\0';                                  // truncate
  }

  // input the run number
  if (runNumber == 0) {
    cout << "     run number? ";    
    fgets(runNo, sizeof(runNo), stdin);
    runNumber = atoi(runNo);
    sprintf(runName, "ana%2d", runNumber);               // add ana prefix
    cout << " run name = " << runName << endl;
  }

  // input the file number (0 opens flow.hist.root)
  if (strstr(fileNumber, "x")!=0) {
    cout << "     anaXX.root file number? [0= flow.hist.root] ";
    fgets(fileNumber, sizeof(fileNumber), stdin);
    fileNumber[strlen(fileNumber)-1] = '\0';
    if (strlen(fileNumber) == 1 && strstr(fileNumber,"0")) {
      sprintf(fileName, "flow.hist.root");
    } else {
      sprintf(fileName, "ana%s.root", fileNumber);       // insert
    }
    cout << " file name = " << fileName << endl;
    histFile = new TFile(fileName);
  }
  
  // input the page number
  while (pageNumber <= 1 || pageNumber > nNames) {
    if (pageNumber < 0) {                                 // plot all
      plotAll(nNames, selN, harN, -pageNumber);
      return can;
    }
    if (pageNumber == 1) {                                // plot resolution
      can = plotResolution();
      return can;
    }
    cout << "-1: \t All" << endl;                         // print menu
    for (int i = 0; i < nNames; i++) {
      cout << i+1 << ":\t " << shortName[i] << endl;
    }
    cout << "     page number? ";
    fgets(tmp, sizeof(tmp), stdin);
    pageNumber = atoi(tmp);
  }
  if (pageNumber > 0 && pageNumber <= nSingles) {         // plot singles
    singleGraph = kTRUE;
    multiGraph  = kFALSE;
  }
  pageNumber--;
  cout << pageNumber+1 << ":  graph name= " << shortName[pageNumber] << endl;

  // set constants
  float twopi   = 2. * TMath::Pi();
  float phiMax  = twopi; 
  float Ycm     =   0.0;
  TString* histProjName = NULL;

  // set row and column numbers
  char* cp = strstr(shortName[pageNumber],"Subs");
  int columns = (cp) ? nSubs * nSels : nSels;
  int rows;
  rows = (strstr(shortName[pageNumber],"Phi_") ||
	  strstr(shortName[pageNumber],"Psi_Diff") ? 2 : nHars);
  //rows = (strstr(shortName[pageNumber],"Phi_") ? 2 : nHars);
  if (strcmp(shortName[pageNumber],"Flow_Phi_Corr")==0) rows = nHars;
  if (strcmp(shortName[pageNumber],"Flow_Psi_Sub_Corr_Diff")==0) rows = 3;
  int pads = rows*columns;
  cout << "pads = " << pads << endl;

  // make the graph page
  if (multiGraph) {
    int canvasWidth = 600, canvasHeight = 780;             // portrait
  } else {
    int canvasWidth = 780, canvasHeight = 600;             // landscape
  }
  TString* canName = new TString(shortName[pageNumber]);
  if (runNumber) {
    *canName += runNumber;
  } 
  cout << canName->Data() << endl;
  can = new TCanvas(canName->Data(), canName->Data(), canvasWidth, canvasHeight);
  can->ToggleEventStatus();
  if (multiGraph) {
    TPaveLabel* title = new TPaveLabel(0.1,0.96,0.9,0.99,shortName[pageNumber]);
    title->Draw();
  }
  TPaveLabel* run = new TPaveLabel(0.1,0.01,0.2,0.03,runName);  
  run->Draw();
  TDatime now;
  TPaveLabel* date = new TPaveLabel(0.7,0.01,0.9,0.03,now.AsString());
  date->Draw();
  TPad* graphPad = new TPad("Graphs","Graphs",0.01,0.05,0.97,0.95);
  graphPad->Draw();
  graphPad->cd();

  if (multiGraph) {               // many graphs on one page
    graphPad->Divide(columns,rows);
    int firstK = 0, firstJ = 0, lastK = columns, lastJ = rows;
  } else if (singleGraph) {       // single graph on a page
    int firstK = 0, firstJ = 0, lastK = 1, lastJ = 1;
  } else {                        // one graph from a multi graph page
    int firstK = selN -1, firstJ = harN -1, lastK = selN, lastJ = harN;
  }
  TLine* lineZeroY  = new TLine(-etaMax, 0., etaMax, 0.);
  TLine* lineYcm    = new TLine(Ycm, -10., Ycm, 10.);
  TLine* lineOnePhi = new TLine(0., 1., phiMax, 1.);
  for (int j = firstJ; j < lastJ; j++) {
    float order = (float)(j+1);
    for (int k = firstK ; k < lastK; k++) {
      int padN = j*columns + k + 1;                    // pad number

      // construct histName and histProjName
      char* temp = new char[35];
      strcpy(temp,shortName[pageNumber]);
      char* cproj = strstr(temp,".");
      if (cproj) {                                     // a projection
	*cproj = '\0';                                 // remove from "." on
	if (singleGraph) {
	  cproj = strstr(temp,"2");
	  if (cproj) {                                 // a 2D projection 
	    *cproj = '\0';                             // remove from "2D" on
	    strcat(temp,"3D");
	  } else {
	    strcat(temp,"2D");
	  }
	} else {
	  strcat(temp,"2D_Sel");
	}
	TString* histName = new TString(temp);
	histProjName = new TString(baseName[pageNumber]);
	if (!singleGraph) {
	  *histProjName += k+1;
	  histProjName->Append("_Har");
	  *histProjName += j+1;
	}
      } else {                                         // not projection
	TString* histName = new TString(baseName[pageNumber]);
      }
      if (!singleGraph) {
	*histName += k+1;
	histName->Append("_Har");
	*histName += j+1;
      }
      cout << " col= " << k+1 << " row= " << order << " pad= " << padN << "\t" 
	   << histName->Data() << endl;

      // get the histogram
      bool twoD = kFALSE;
      bool threeD = kFALSE;
      if (histProjName) {
	if (strstr(temp,"3D")) {                      // 2D projection
	  twoD = kTRUE;
	  TH3* hist3D = dynamic_cast<TH3*>(histFile->Get(histName->Data()));
	  if (!hist3D) {
	    cout << "### Can't find histogram " << histName->Data() << endl;
	    return can;
	  }
	} else {                                      // 1D projection
	  TH2* hist2D = dynamic_cast<TH2*>(histFile->Get(histName->Data()));
	  if (!hist2D) {
	    cout << "### Can't find histogram " << histName->Data() << endl;
	    return can;
	  }
	}
      } else {
	if (strstr(shortName[pageNumber],"3D")!=0) {  // 3D
	  threeD = kTRUE;
	  TH3* hist3D = dynamic_cast<TH3*>(histFile->Get(histName->Data()));
	  if (!hist3D) {
	    cout << "### Can't find histogram " << histName->Data() << endl;
	    return can;
	  }
	} else if (strstr(shortName[pageNumber],"2D")!=0) { // 2D
	  twoD = kTRUE;
	  TH2* hist2D = dynamic_cast<TH2*>(histFile->Get(histName->Data()));
	  if (!hist2D) {
	    cout << "### Can't find histogram " << histName->Data() << endl;
	    return can;
	  }
	} else {                                            // 1D
	  TH1* hist = dynamic_cast<TH1*>(histFile->Get(histName->Data()));
	  if (!hist) {
	    cout << "### Can't find histogram " << histName->Data() << endl;
	    return can;
	  }
	  float max = hist->GetXaxis()->GetXmax();
	}
      }
      
      // make the plots
      if (multiGraph) graphPad->cd(padN);
      if (threeD) {                                            // 3D
	gStyle->SetOptStat(10);
	hist3D->Draw("BOX");
      } else if (twoD) {                                       // 2D
	if (strstr(shortName[pageNumber],".PhiEta")!=0) {      // 3D Phi Eta proj.
	  TH2D* projZX = hist3D->Project3D("zx");
	  projZX->SetName(histProjName->Data());
	  projZX->SetYTitle("azimuthal angle (rad)");
	  projZX->SetXTitle("rapidity");
	  gStyle->SetOptStat(0);
	  if (projZX) projZX->Draw("COLZ");
	} else if (strstr(shortName[pageNumber],".PhiPt")!=0) { // 3D Phi Pt proj.
	  TH2D* projZY = hist3D->Project3D("zy");
	  projZY->SetName(histProjName->Data());
	  projZY->SetYTitle("azimuthal angle (rad");
	  projZY->SetXTitle("Pt (GeV/c)");
	  gStyle->SetOptStat(0);
	  if (projZY) projZY->Draw("COLZ");
	} else	if ((strstr(shortName[pageNumber],"QXY")!=0)) { // Q XY
	  TLine* lineZeroX = new TLine(-0.5, 0., 0.5, 0.);
	  TLine* lineZeroY = new TLine(0., -0.5, 0., 0.5);
	  gStyle->SetOptStat(100110);
	  hist2D->Draw("COLZ");
	  //hist2D->Draw("CONT");
	  lineZeroX->Draw();
	  lineZeroY->Draw();
	} else	if ((strstr(shortName[pageNumber],"TPCSubXY")!=0)) { // QSub XY
	  TLine* lineZeroX = new TLine(-1., 0., 1., 0.);
	  TLine* lineZeroY = new TLine(0., -1., 0., 1.);
	  gStyle->SetOptStat(100110);
	  hist2D->Draw("COLZ");
	  //hist2D->Draw("CONT");
	  lineZeroX->Draw();
	  lineZeroY->Draw();
	} else	if (strstr(shortName[pageNumber],"XY")!=0) {    // Vertex XY
	  TLine* lineZeroX = new TLine(-1., 0., 1., 0.);
	  TLine* lineZeroY = new TLine(0., -1., 0., 1.);
	  gStyle->SetOptStat(10);
	  hist2D->Draw("COLZ");
	  lineZeroX->Draw();
	  lineZeroY->Draw();
	} else if (strstr(shortName[pageNumber],"Dedx")!=0) {   // dE/dx
	  TVirtualPad::Pad()->SetLogz();
	  gStyle->SetOptStat(10);
	  hist2D->Draw("COLZ");
	} else if (strstr(shortName[pageNumber],"_v")!=0) {    // v
	  hist2D->SetMaximum(20.);
	  hist2D->SetMinimum(-20.);
	  gStyle->SetOptStat(0);
	  hist2D->Draw("COLZ");
	} else {                                               // other 2D
	  gStyle->SetOptStat(10);
	  hist2D->Draw("COLZ");
	}
      } else if (strstr(shortName[pageNumber],".Eta")!=0) { // 2D Eta projection
	if (singleGraph) {
	  TH1D* projX = hist2D->ProjectionX(histName->Data());
	} else {
	  //TH1D* projX = hist2D->ProjectionX(histName->Data(), -1, 9999, "E");
	  TH1D* projX = hist2D->ProjectionX(histName->Data());
	}
	projX->SetName(histProjName->Data());
	char* xTitle = hist2D->GetXaxis()->GetTitle();
	projX->SetXTitle(xTitle);
	projX->SetYTitle("Counts");
	gStyle->SetOptStat(0);
	if (projX) projX->Draw("H");
	if (!singleGraph) lineZeroY->Draw();
      } else if (strstr(shortName[pageNumber],".Pt")!=0) { // 2D Pt projection
	if (singleGraph) {
	  TH1D* projY = hist2D->ProjectionY(histName->Data()); 
	} else {
	  TH1D* projY = hist2D->ProjectionY(histName->Data(), -1, 9999, "E");
	}
	projY->SetName(histProjName->Data());
	projY->SetXTitle("Pt (GeV/c)");
	projY->SetYTitle("Counts");
	TVirtualPad::Pad()->SetLogy();
	gStyle->SetOptStat(0);
	if (projY) projY->Draw("H");
      } else if (strstr(shortName[pageNumber],".Phi")!=0) { // 2D Phi projection
	if (singleGraph) {
	  TH1D* projY = hist2D->ProjectionY(histName->Data()); 
	} else {
	  TH1D* projY = hist2D->ProjectionY(histName->Data(), -1, 9999, "E");
	}
	projY->SetName(histProjName->Data());
	projY->SetXTitle("azimuthal angle (rad");
	projY->SetYTitle("Counts");
	gStyle->SetOptStat(0);
	if (projY) projY->Draw("H");
      } else if (strstr(shortName[pageNumber],"Corr")!=0) { // azimuthal corr.
	float norm = (hist->Integral()) ? ((float)(hist->GetNbinsX()) / hist->Integral()) : 1.; 
	cout << "  Normalized by: " << norm << endl;
	hist->Scale(norm);                           // normalize height to one
	if (strstr(shortName[pageNumber],"Diff")!=0) { 
	  TF1* funcCos1 = new TF1("funcCos1",
				  "1+[0]*2/100*cos([1]*x)", 0., twopi);
	  funcCos1->SetParNames("k=1", "har");
	  funcCos1->SetParameters(0, order+1);              // initial values
	  funcCos1->SetParLimits(1, 1, 1);                  // har is fixed
	  //hist->Fit("funcCos1");
	  delete funcCos1;
	} else if (strstr(shortName[pageNumber],"Sub")!=0) { 
	  TF1* funcSubCorr = new TF1("SubCorr", SubCorr, 0., twopi/order, 2);
	  funcSubCorr->SetParNames("chi", "har");
	  funcSubCorr->SetParameters(1., order);             // initial value
	  funcSubCorr->SetParLimits(1, 1, 1);                // har is fixed
	  hist->Fit("SubCorr");
	  delete funcSubCorr;
	} else {
	  TF1* funcCos3 = new TF1("funcCos3",
	 "1+[0]*2/100*cos([3]*x)+[1]*2/100*cos(([3]+1)*x)+[2]*2/100*cos(([3]+2)*x)",
				  0., twopi/order);
	  funcCos3->SetParNames("n=har", "n=har+1", "n=har+2", "har");
	  funcCos3->SetParameters(0, 0, 0, order);          // initial values
	  funcCos3->SetParLimits(3, 1, 1);                  // har is fixed
	  //hist->Fit("funcCos3");
	  delete funcCos3;
	}
	if (strstr(shortName[pageNumber],"Phi")!=0)
	  hist->SetMinimum(0.9*(hist->GetMinimum()));
	gStyle->SetOptStat(10);
	gStyle->SetOptFit(111);
	hist->Draw("E1");
      } else if (strstr(shortName[pageNumber],"_q")!=0) {   // q distibution
	gStyle->SetOptStat(110);
	gStyle->SetOptFit(111);
	hist->Draw("E1");
	float n_qBins = (float)hist->GetNbinsX();
	double area = hist->Integral() * max / n_qBins; 
	TString* histName = new TString("Flow_Mul_Sel");
	*histName += k+1;
	histName->Append("_Har");
	*histName += j+1;
	TH1* histMult = dynamic_cast<TH1*>(histFile->Get(histName->Data()));
	if (!histMult) {
	  cout << "### Can't find histogram " << histName->Data() << endl;
	  return can;
	}
	delete histName;
	float mult = histMult->GetMean();
	TF1* fit_q = new TF1("qDist", qDist, 0., max, 4);
	fit_q->SetParNames("v", "mult", "area", "g");
	float qMean = hist->GetMean();
	float vGuess = (qMean > 1.) ? sqrt(2.*(qMean - 1.) / mult) : 0.03;
	// the 0.03 is a wild guess
	vGuess *= 100.;
	cout << "vGuess = " << vGuess << endl;
	fit_q->SetParameters(vGuess, mult, area, 0.3); // initial values
	fit_q->SetParLimits(1, 1, 1);             // mult is fixed
	fit_q->SetParLimits(2, 1, 1);             // area is fixed
	hist->Fit("qDist", "Q");
	//hist->Fit("qDist");
	fit_q->Draw();
 	fit_q->FixParameter(3, 0.);               // g is fixed
	fit_q->SetLineStyle(kDotted);
	hist->Fit("qDist", "Q+");
	fit_q->Draw("same");
 	fit_q->ReleaseParameter(3);               // g is unfixed
 	fit_q->SetParameter(3, 0.5);              // initial value
 	fit_q->FixParameter(0, 0.);               // v is fixed
	fit_q->SetLineStyle(kDashed);
	hist->Fit("qDist", "Q+");
	fit_q->Draw("same");
 	fit_q->ReleaseParameter(0);               // v is unfixed
      } else if (strstr(shortName[pageNumber],"PhiLab")!=0) {  // Phi lab distibution
	hist->Draw("E1");
	float norm = (hist->Integral()) ? ((float)(hist->GetNbinsX()) / hist->Integral()) : 1.; 
	cout << "  Normalized by: " << norm << endl;
	hist->Scale(norm);                           // normalize height to one
        hist->SetMinimum(0.);
        hist->SetMaximum(1.3);
	TF1* funcCosSin = new TF1("funcCosSin",
		   "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, order); // initial values
	funcCosSin->SetParLimits(2, 1, 1); // har is fixed
	hist->Fit("funcCosSin");
	delete funcCosSin;
	gStyle->SetOptFit(111);
      } else if (strstr(shortName[pageNumber],"Phi")!=0) {  // other Phi distibutions
       	//hist->SetMinimum(0.9*(hist->GetMinimum()));
       	hist->SetMinimum(0.);
	gStyle->SetOptStat(10);
	hist->Draw(); 
      } else if (strstr(shortName[pageNumber],"Psi_Diff")!=0) { // Psi_Diff distibutions
	hist->Draw("E1");
      } else if (strstr(shortName[pageNumber],"Psi")!=0) {    // Psi distibutions
 	float norm = (hist->Integral()) ? ((float)(hist->GetNbinsX()) / hist->Integral()) : 1.; 
	cout << "  Normalized by: " << norm << endl;
	hist->Scale(norm);                           // normalize height to one
      	hist->SetMinimum(0.);
	TF1* funcCosSin = new TF1("funcCosSin",
				"1.+[0]*2./100.*cos([2]*x)+[1]*2./100.*sin([2]*x)", 0., twopi/order);
	funcCosSin->SetParNames("100*cos", "100*sin", "har");
	funcCosSin->SetParameters(0, 0, order); // initial values
	funcCosSin->SetParLimits(2, 1, 1); // har is fixed
	hist->Fit("funcCosSin");
	delete funcCosSin;
	gStyle->SetOptFit(111);
	hist->Draw("E1");
      } else if (strstr(shortName[pageNumber],"Eta")!=0) {    // Eta distibutions
	gStyle->SetOptStat(100110);
 	if (strstr(shortName[pageNumber],"_v")!=0 ) {
	  hist->SetMaximum(10.);
	  hist->SetMinimum(-10.);
	  hist->Draw();
	  lineYcm->Draw();
	} else {
	  hist->Draw();
	}
	lineZeroY->Draw();
      } else if (strstr(shortName[pageNumber],"Pt")!=0) {     // Pt distibutions
 	if (strstr(shortName[pageNumber],"_v")!=0 ) {
	  hist->SetMaximum(25.);
	  hist->SetMinimum(-5.);
	}
	gStyle->SetOptStat(100110);
	hist->Draw();
	if (strstr(shortName[pageNumber],"v")!=0) {
	  TLine* lineZeroPt  = new TLine(0., 0., max, 0.);
	  lineZeroPt->Draw();
	}
      } else if (strstr(shortName[pageNumber],"Bin")!=0) {    // Bin hists
	if (strstr(shortName[pageNumber],"Pt")!=0) {
	  TLine* lineDiagonal = new TLine(0., 0., max, max);
	} else {
	  TLine* lineDiagonal = new TLine(-max, -max, max, max);
	}
	gStyle->SetOptStat(0);
	hist->SetMarkerStyle(21);
	hist->SetMarkerColor(2);
	hist->Draw();
	lineDiagonal->Draw();
      } else if (strstr(shortName[pageNumber],"PidMult")!=0) {  // PID Mult
	TVirtualPad::Pad()->SetLogy();
	gStyle->SetOptStat(0);
	hist->Draw();
      } else {                                              // all other 1D
	gStyle->SetOptStat(100110);
	hist->Draw(); 
      }
      delete [] temp;
      delete histName;
      if (histProjName) delete histProjName;
    }
  }
  delete [] shortName;

  return can;
}

// macro for the resolution plot
TCanvas* plotResolution() {
  char* resName[] = {
    "Flow_Cos_Sel",
    "Flow_Res_Sel",
    "Flow_v_Sel"
//     "Flow_v_ScalarProd_Sel",
//     "Flow_Cumul_v_Order2_Sel",
//     "Flow_Cumul_v_Order4_Sel"
  };
  int columns = nSels;
  int rows = sizeof(resName) / sizeof(char*);
  int pads = rows*columns;

  // make the graph page
  can = new TCanvas(resName[1], resName[1], 600, 780);
  can->ToggleEventStatus();
  TPaveLabel* run = new TPaveLabel(0.1,0.01,0.2,0.03,runName);  
  run->Draw();
  TDatime now;
  TPaveLabel* date = new TPaveLabel(0.7,0.01,0.9,0.03,now.AsString());
  date->Draw();
  TLine* lineZeroHar = new TLine(0.5, 0., nHars+0.5, 0.);
  TPad* graphPad = new TPad("Graphs","Graphs",0.01,0.05,0.97,0.99);
  graphPad->Draw();
  graphPad->cd();
  graphPad->Divide(columns,rows);

  // make the plots
  float v;
  float err;
  for (int j = 0; j < rows; j++) {
    int resNumber = j;
    cout << "resolution name= " << resName[resNumber] << endl;
    for (int k = 0; k < columns; k++) {
      int padN = j*columns + k +1;
      TString* histName = new TString(resName[resNumber]);
      *histName += k+1;
      cout << "row= " << j << " col= " << k << " pad= " << padN << "\t" 
	   << histName->Data() << endl;
      TH1* hist = dynamic_cast<TH1*>(histFile->Get(histName->Data()));
      if (!hist) {
	cout << "### Can't find histogram " << histName->Data() << endl;
	return can;
      }
      graphPad->cd(padN);
      gStyle->SetOptStat(0);
      if (strstr(resName[resNumber],"_v")!=0) {
	hist->SetMaximum(10.);
	hist->SetMinimum(-5.);
      } else {
	hist->SetMaximum(1.1);
	hist->SetMinimum(0.);
      }
      for (int n=1; n < nHars+1; n++) {
	v   = hist->GetBinContent(n);                       // output v values
	err = hist->GetBinError(n);
	cout << " " << n << ": " << setprecision(3) << v << " +/- " << err << endl;
	if (TMath::IsNaN(v)) {
	  hist->SetBinContent(n, 0.);
	  hist->SetBinError(n, 0.);
	}
      }
      hist->Draw();
      lineZeroHar->Draw();
      delete histName;
    }
  }

  return can;
} 

void plotAll(Int_t nNames, Int_t selN, Int_t harN, Int_t first = 1) {
  for (int i =  first; i < nNames + 1; i++) {
    can = plot(i, selN, harN);
    can->Update();
    cout << "save? y/[n], quit? q" << endl;
    fgets(tmp, sizeof(tmp), stdin);
    if (strstr(tmp,"y")!=0) can->Print(".pdf");
    else if (strstr(tmp,"q")!=0) return;
    else if (strstr(tmp," ")!=0) continue;
  }
  cout << "  plotAll Done" << endl;
}

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

static Double_t qDist(double* q, double* par) {
  // Calculates the q distribution given the parameters v, mult, area, g

  double sig2 = 0.5 * (1. + par[3]);
  double expo = (par[1]*par[0]*par[0]/10000. + q[0]*q[0]) / (2*sig2);
  Double_t dNdq = par[2] * (q[0]*exp(-expo)/sig2) * 
    TMath::BesselI0(q[0]*par[0]/100.*sqrt(par[1])/sig2);

  return dNdq;
}

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

static Double_t SubCorr(double* x, double* par) {
  // Calculates the n(Psi_a - Psi_b) distribution by fitting chi
  // From J.-Y. Ollitrault, Nucl. Phys. A590, 561c (1995), Eq. 6. with correc.

  double chi2 = par[0] * par[0] / 2;     // divide by two for SV chi
  double z = chi2 * cos(par[1]*x[0]);
  double TwoOverPi = 2./TMath::Pi();

  Double_t dNdPsi = exp(-chi2)/TwoOverPi * (TwoOverPi*(1.+chi2) 
                    + z*(TMath::BesselI0(z) + TMath::StruveL0(z))
                    + chi2*(TMath::BesselI1(z) + TMath::StruveL1(z)));

  return dNdPsi;
}

///////////////////////////////////////////////////////////////////////////////
//
// $Log: plot.C,v $
// Revision 1.70  2011/07/25 15:54:53  posk
// Added correction for non-flatness of event plane.
//
// Revision 1.69  2011/03/10 18:56:39  posk
// Added histogram for laboratory azimuthal distribution of particles.
//
// Revision 1.68  2010/09/30 19:28:23  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.67  2010/03/05 17:04:38  posk
// ROOT 5.22 compatable.
// Moved doFlowEvents.C here from StRoot/macros/analysis/
//
// Revision 1.66  2009/11/24 19:29:19  posk
// Added reCenter to remove acceptance correlations as an option instead of phiWgt.
//
// Revision 1.65  2006/02/22 19:35:18  posk
// Added graphs for the StFlowLeeYangZerosMaker
//
// Revision 1.64  2005/08/26 19:00:25  posk
// plot style back to bold
//
// Revision 1.63  2005/08/05 20:13:43  posk
// Improved first guess for qDist fit.
//
// Revision 1.62  2004/12/09 23:47:11  posk
// Minor changes in code formatting.
// Added hist for TPC primary dca to AnalysisMaker.
//
// Revision 1.61  2004/12/07 23:10:23  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.60  2004/12/02 16:10:55  posk
// Added  gInterpreter->ProcessLine(".O0");
//
// Revision 1.59  2004/11/19 16:54:40  posk
// Replaced gPad with (TVirtualPad::Pad()). Reverted to TMath::Struve functions.
//
// Revision 1.58  2004/11/11 18:25:43  posk
// Minor updates.
//
// Revision 1.57  2004/03/11 18:00:05  posk
// Added Random Subs analysis method.
//
// Revision 1.56  2004/03/01 22:43:43  posk
// Changed some "->" to ".".
//
// Revision 1.55  2003/08/26 21:10:13  posk
// Calculates v8 if nHars=8.
//
// Revision 1.54  2003/07/30 22:03:26  oldi
// Fit("pol4") taken out (which was in accidentally).
//
// Revision 1.53  2003/07/07 21:58:21  posk
// Made units of momentum GeV/c instead of GeV.
//
// Revision 1.52  2003/06/27 21:25:44  posk
// v4 and v6 are with repect to the 2nd harmonic event plane.
//
// Revision 1.51  2003/05/06 21:33:07  posk
// Removed some histograms.
//
// Revision 1.50  2003/05/02 21:11:13  posk
// Reduced the number of harmonics from 3 to 2.
//
// Revision 1.49  2003/03/18 17:58:36  posk
// Kirill Fillimonov's improved fit to the angle between subevent planes.
//
// Revision 1.48  2003/03/17 20:46:55  posk
// Improved fit to q dist.
//
// Revision 1.47  2003/03/11 23:04:30  posk
// Includes scalar product hists.
//
// Revision 1.46  2003/02/25 19:25:32  posk
// Improved plotting.
//
// Revision 1.45  2003/02/05 18:52:50  posk
// Added Bool_t includeFtpc
//
// Revision 1.44  2003/01/16 16:02:30  posk
// Some plotting changes.
//
// Revision 1.43  2003/01/10 16:40:53  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.42  2002/11/26 22:11:54  posk
// First use of doxygen.
//
// Revision 1.41  2002/06/11 21:54:16  posk
// Kirill's further correction to minBias.C for bins with one count.
//
// Revision 1.40  2002/05/21 18:42:18  posk
// Kirill's correction to minBias.C for bins with one count.
//
// Revision 1.39  2002/02/13 22:31:50  posk
// Pt Weight now also weights Phi Weight. Added Eta Weught, default=FALSE.
//
// Revision 1.38  2002/01/14 23:42:57  posk
// Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
//
// Revision 1.37  2001/12/18 19:27:37  posk
// "proton" and "antiproton" replaced by "pr+" and "pr-".
//
// Revision 1.36  2001/12/11 22:04:13  posk
// Four sets of phiWgt histograms.
// StFlowMaker StFlowEvent::PhiWeight() changes.
// Cumulant histogram names changed.
//
// Revision 1.35  2001/11/13 22:47:35  posk
// Documentation updated. Fit to q function moved to macro.
//
// Revision 1.34  2001/11/09 21:15:04  posk
// Switched from CERNLIB to TMath. Using global dca instead of dca.
//
// Revision 1.33  2001/05/22 20:11:20  posk
// Changed dEdx graphs.
//
// Revision 1.32  2000/12/12 15:01:11  posk
// Put log comments at end of file.
//
// Revision 1.31  2000/12/08 17:04:09  oldi
// Phi weights for both FTPCs included.
//
// Revision 1.27  2000/09/26 20:54:12  posk
// Updated documentation.
//
// Revision 1.26  2000/09/15 22:52:56  posk
// Added Pt weighting for event plane calculation.
//
// Revision 1.25  2000/08/31 18:50:32  posk
// Added plotCen.C to plot from a series of files with different centralities.
//
// Revision 1.24  2000/08/12 20:20:15  posk
// More centrality bins.
//
// Revision 1.23  2000/08/01 21:51:20  posk
// Added doubly integrated v.
//
// Revision 1.22  2000/07/12 17:49:39  posk
// Changed EtaSym plots.
//
// Revision 1.21  2000/06/30 14:51:20  posk
// Using MessageMgr. Added graph for Eta Symmetry vs. Vertex Z.
//
// Revision 1.20  2000/05/26 21:25:23  posk
// Use TProfile2D class and profile projection methods.
// Correction needed for >2 subevents.
//
// Revision 1.18  2000/04/13 22:34:16  posk
// Resolution correction is now made.
//
// Revision 1.17  2000/04/10 18:49:10  posk
// Asks for the histogram file number.
//
// Revision 1.16  2000/03/28 23:25:37  posk
// Allow multiple instances.
//
// Revision 1.15  2000/03/21 00:24:45  posk
// Added GetCVS and changed some plot names.
//
// Revision 1.12  2000/02/18 23:44:54  posk
// Added PID and centrality.
//
// Revision 1.11  2000/02/04 16:26:43  posk
// Added correct calculation of event plane resolution for large flow.
//
// Revision 1.2  1999/10/05 16:54:14  posk
// Added getPhiWeight method for making the event plane isotropic.
//
///////////////////////////////////////////////////////////////////////////////