StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plotCen.C
1 //
3 // $Id: plotCen.C,v 1.29 2010/09/30 19:28:25 posk Exp $
4 //
5 // Author: Art Poskanzer, LBNL, July 2000
6 // FTPC added by Markus Oldenburg, MPI, Dec 2000
7 // Description: Macro to plot histograms made by StFlowAnalysisMaker.
8 // Plots a set of histograms with different centrality
9 // starting with anaXX.root given by first run number XX.
10 // Run Number appended to "ana" is entered in the bottom, left box.
11 // Default selN = 2 and harN = 2.
12 // First time type .x plotCen.C() to see the menu.
13 // After the first execution, just type plotCen(N) .
14 // A negative N plots all pages starting with page N.
15 // Place a symbolic link to this file in StRoot/macros/analysis .
16 //
17 //
19 #include <iomanip.h>
20 #include <math.h>
21 #include "TMath.h"
22 
23 const Int_t nCens = 10; // min bias + 9 centralities
24 //const Int_t nCens = 9; // 9 centralities
25 int runNumber = 0;
26 char runName[60];
27 char fileName[60];
28 char histTitle[30];
29 TFile* histFile[nCens];
30 char tmp[10];
31 TCanvas* can;
32 
33 TCanvas* plotCen(Int_t pageNumber=0, Int_t selN=2, Int_t harN=2){
34  gInterpreter->ProcessLine(".O0");
35 
36  TCanvas* cOld = (TCanvas*)gROOT->GetListOfCanvases(); // delete old canvas
37  if (cOld) cOld->Delete();
38 
39  //gROOT->SetStyle("Pub"); // set style
40  gROOT->SetStyle("Bold"); // set style
41  gStyle->SetPalette(1);
42  gROOT->ForceStyle();
43  gStyle->SetLabelSize(.1,"X");
44 
45  int canvasWidth = 600, canvasHeight = 780; // portrait
46  int columns = 2;
47  int rows;
48  bool oddPads = (nCens) % 2;
49  if (oddPads) {
50  rows = nCens/columns + 1;
51  } else {
52  rows = nCens/columns;
53  }
54  int pads = nCens;
55 
56  // names of histograms made by StFlowAnalysisMaker
57  // also projections of some of these histograms
58  const char* baseName[] = {
59  "Flow_Cos_Sel",
60  "Flow_Res_Sel",
61  "Flow_v_Sel",
62  "FlowLYZ_V_Sel",
63  "FlowLYZ_vr0_Sel",
64  "FlowLYZ_v_Sel",
65  "Flow_v_ScalarProd_Sel",
66  "Flow_Cumul_v_Order2_Sel",
67  "Flow_Cumul_v_Order4_Sel",
68  "Flow_VertexZ",
69  "Flow_VertexXY2D",
70  "Flow_EtaSymVerZ2D_Tpc",
71  "Flow_EtaSym_Tpc",
72  "Flow_EtaSymVerZ2D_Ftpc",
73  "Flow_EtaSym_Ftpc",
74  "Flow_CTBvsZDC2D",
75  "Flow_Cent",
76  "Flow_OrigMult",
77  "Flow_Mult",
78  "FlowLYZ_Mult",
79  "Flow_MultOverOrig",
80  "Flow_MultEta",
81  "Flow_MultPart",
82  "Flow_Charge_Ftpc",
83  "Flow_DcaGlobal_Tpc",
84  "Flow_Dca_Ftpc",
85  "Flow_DcaGlobal_Ftpc",
86  "Flow_Chi2_Tpc",
87  "Flow_Chi2_Ftpc",
88  "Flow_FitPts_Tpc",
89  "Flow_MaxPts_Tpc",
90  "Flow_FitOverMax_Tpc",
91  "Flow_FitPts_Ftpc",
92  "Flow_MaxPts_Ftpc",
93  "Flow_FitOverMax_Ftpc",
94  "Flow_YieldAll2D",
95  "Flow_YieldAll.Eta",
96  "Flow_YieldAll.Pt",
97  "Flow_YieldPart2D",
98  "Flow_YieldPart.Eta",
99  "Flow_YieldPart.Pt",
100  "Flow_MeanDedxPos2D",
101  "Flow_MeanDedxNeg2D",
102 // "Flow_PidPiPlusPart",
103 // "Flow_PidPiMinusPart",
104 // "Flow_PidProtonPart",
105 // "Flow_PidAntiProtonPart",
106 // "Flow_PidKplusPart",
107 // "Flow_PidKminusPart",
108 // "Flow_PidDeuteronPart",
109 // "Flow_PidAntiDeuteronPart",
110 // "Flow_PidElectronPart",
111 // "Flow_PidPositronPart",
112  "Flow_PidMult",
113 // "Flow_CosPhiLab",
114  "Flow_Yield2D_Sel",
115  "Flow_Yield.Eta_Sel",
116  "Flow_Yield.Pt_Sel",
117  "Flow_Mul_Sel",
118 // "Flow_Phi_East_Sel",
119 // "Flow_Phi_Flat_East_Sel",
120 // "Flow_Phi_West_Sel",
121 // "Flow_Phi_Flat_West_Sel",
122 // "Flow_Phi_FtpcEast_Sel",
123 // "Flow_Phi_Flat_FtpcEast_Sel",
124 // "Flow_Phi_FtpcWest_Sel",
125 // "Flow_Phi_Flat_FtpcWest_Sel",
126  "Flow_Psi_Subs",
127  "Flow_Psi_Sel",
128  "Flow_Psi_Sub_Corr_Sel",
129  "Flow_Psi_Sub_Corr_Diff_Sel",
130  "Flow_Phi_Corr_Sel",
131  "Flow_vObs2D_Sel",
132  "Flow_vObsEta_Sel",
133  "Flow_vObsPt_Sel",
134  "Flow_v2D_Sel",
135  "Flow_vEta_Sel",
136  "Flow_vPt_Sel",
137  "Flow_q_Sel",
138  "FlowLYZ_r0theta_Sel",
139  "FlowLYZ_Vtheta_Sel",
140  "FlowLYZ_vEta_Sel",
141  "FlowLYZ_vPt_Sel",
142  "FlowLYZ_Gtheta0_Sel",
143  "Flow_vObs2D_ScalarProd_Sel",
144  "Flow_v2D_ScalarProd_Sel",
145  "Flow_vEta_ScalarProd_Sel",
146  "Flow_vPt_ScalarProd_Sel",
147  "Flow_Cumul_vEta_Order2_Sel",
148  "Flow_Cumul_vPt_Order2_Sel",
149  "Flow_Cumul_vEta_Order4_Sel",
150  "Flow_Cumul_vPt_Order4_Sel"
151 };
152  int nName = sizeof(baseName) / sizeof(char*);
153  const Int_t nNames = nName;
154  const int nDoubles = 9;
155  const int nSingles = 46 + nDoubles;
156 
157  // construct array of short names
158  char* shortName[nNames];
159  for (int n = 0; n < nNames; n++) {
160  shortName[n] = new char[35];
161  strcpy(shortName[n], baseName[n]);
162  char* cp = strstr(shortName[n],"_Sel");
163  if (cp) *cp = '\0'; // truncate
164  }
165 
166  cout << "Harmonic = " << harN << endl << endl;
167  // input the first run number
168  if (runNumber == 0) {
169  cout << " first run number? ";
170  fgets(tmp, sizeof(tmp), stdin);
171  runNumber = atoi(tmp);
172  sprintf(runName, "ana%2d", runNumber); // add ana prefix
173  cout << " first run name = " << runName << endl;
174 
175  // open the files
176  for (int n = 0; n < nCens; n++) {
177  sprintf(fileName, "ana%2d.root", runNumber + n);
178  cout << " file name = " << fileName << endl;
179  histFile[n] = new TFile(fileName);
180  }
181  }
182 
183  // input the page number
184  while (pageNumber <= 0 || pageNumber > nNames) {
185  if (pageNumber < 0) { // plot all
186  plotCenAll(nNames, selN, harN, -pageNumber);
187  return can;
188  }
189  cout << "-1: \t All" << endl; // print menu
190  for (int i = 0; i < nNames; i++) {
191  cout << i+1 << ":\t " << shortName[i] << endl;
192  }
193  cout << " page number? ";
194  fgets(tmp, sizeof(tmp), stdin);
195  pageNumber = atoi(tmp);
196  }
197 
198  // set flags
199  bool multiGraph = kFALSE;
200  bool doubleGraph = kFALSE;
201  bool singleGraph = kFALSE;
202  if (pageNumber > 0 && pageNumber <= nDoubles) {
203  doubleGraph = kTRUE;
204  } else if (pageNumber > nDoubles && pageNumber <= nSingles) {
205  singleGraph = kTRUE;
206  } else {
207  multiGraph = kTRUE;
208  }
209  pageNumber--;
210 
211  // set constants
212  float twopi = 2. * TMath::Pi();
213  float etaMax = 1.5;
214  float yMin = -4.5;
215  float yMax = 4.5;
216  float qMax = 3.5;
217  float phiMax = twopi;
218  int n_qBins = 50;
219  float Ycm = 0.0;
220  TString* histProjName = NULL;
221 
222  // construct histName and histProjName
223  char sel[2];
224  sprintf(sel,"%d",selN);
225  char har[2];
226  sprintf(har,"%d",harN);
227  float order = (float)harN;
228  char* temp = new char[35]; // construct histName
229  strcpy(temp,shortName[pageNumber]);
230  char* cproj = strstr(temp,".");
231  if (cproj) { // a projection
232  *cproj = '\0'; // remove from "." on
233  if (singleGraph) {
234  cproj = strstr(temp,"2");
235  if (cproj) { // a 2D projection
236  *cproj = '\0'; // remove from "2D" on
237  strcat(temp,"3D");
238  } else {
239  strcat(temp,"2D");
240  }
241  } else {
242  strcat(temp,"2D_Sel");
243  }
244  TString* histName = new TString(temp);
245  histProjName = new TString(baseName[pageNumber]);
246  if (multiGraph) {
247  histProjName->Append(*sel);
248  histProjName->Append("_Har");
249  histProjName->Append(*har);
250  }
251  } else { // not projection
252  TString* histName = new TString(baseName[pageNumber]);
253  }
254  if (!singleGraph) histName->Append(*sel);
255  if (multiGraph) {
256  histName->Append("_Har");
257  histName->Append(*har);
258  }
259  cout << pageNumber+1 << ": graph name= " << histName->Data() << endl;
260 
261  // make the graph page
262  can = new TCanvas(shortName[pageNumber], shortName[pageNumber],
263  canvasWidth, canvasHeight);
264  can->ToggleEventStatus();
265  TPaveLabel* title = new TPaveLabel(0.1,0.96,0.9,0.99,histName->Data());
266  title->Draw();
267  TPaveLabel* run = new TPaveLabel(0.1,0.01,0.2,0.03,runName);
268  run->Draw();
269  TDatime now;
270  TPaveLabel* date = new TPaveLabel(0.7,0.01,0.9,0.03,now.AsString());
271  date->Draw();
272  TPad* graphPad = new TPad("Graphs","Graphs",0.01,0.05,0.97,0.95);
273  graphPad->Draw();
274  graphPad->cd();
275  graphPad->Divide(columns, rows);
276  TLine* lineZeroY = new TLine(yMin, 0., yMax, 0.);
277  TLine* lineYcm = new TLine(Ycm, -10., Ycm, 10.);
278  float v;
279  float err;
280  int centr;
281  for (int i = 0; i < pads; i++) {
282  int fileN = i; // file number
283  int padN = fileN + 1; // pad number
284  centr = oddPads ? padN : padN-1;
285  sprintf(histTitle,"Centrality %d",centr);
286  cout << "centrality= " << centr << endl;
287 
288  // get the histogram
289  bool twoD;
290  bool threeD;
291  if (histProjName) {
292  if (strstr(temp,"3D")) { // 2D projection
293  twoD = kTRUE;
294  TH3* hist3D = (TH3*)(histFile[fileN]->Get(histName->Data()));
295  if (!hist3D) {
296  cout << "### Can't find histogram " << histName->Data() << endl;
297  return can;
298  }
299  hist3D->SetTitle(histTitle);
300  } else { // 1D projection
301  TH2* hist2D = (TH2*)(histFile[fileN]->Get(histName->Data()));
302  if (!hist2D) {
303  cout << "### Can't find histogram " << histName->Data() << endl;
304  return can;
305  }
306  hist2D->SetTitle(histTitle);
307  }
308  } else {
309  if (strstr(shortName[pageNumber],"3D")!=0) { // 3D
310  threeD = kTRUE;
311  TH3* hist3D = (TH3*)(histFile[fileN]->Get(histName->Data()));
312  if (!hist3D) {
313  cout << "### Can't find histogram " << histName->Data() << endl;
314  return can;
315  }
316  hist3D->SetTitle(histTitle);
317  } else if (strstr(shortName[pageNumber],"2D")!=0) { // 2D
318  twoD = kTRUE;
319  TH2* hist2D = (TH2*)(histFile[fileN]->Get(histName->Data()));
320  if (!hist2D) {
321  cout << "### Can't find histogram " << histName->Data() << endl;
322  return can;
323  }
324  hist2D->SetTitle(histTitle);
325  } else { // 1D
326  TH1* hist = (TH1*)(histFile[fileN]->Get(histName->Data()));
327  if (!hist) {
328  cout << "### Can't find histogram " << histName->Data() << endl;
329  return can;
330  }
331  float ptMax = hist->GetXaxis()->GetXmax();
332  hist->SetTitle(histTitle);
333  }
334  }
335 
336  // make the plots
337  graphPad->cd(padN);
338 
339  if (threeD) { // 3D
340  gStyle->SetOptStat(10);
341  hist3D->Draw("BOX");
342  } else if (twoD) { // 2D
343  if (strstr(shortName[pageNumber],".PhiEta")!=0) { // 3D Phi Eta proj.
344  TH2D* projZX = (TH2*)hist3D->Project3D("zx");
345  projZX->SetName(histProjName->Data());
346  projZX->SetYTitle("azimuthal angle (rad)");
347  projZX->SetXTitle("rapidity");
348  gStyle->SetOptStat(0);
349  if (projZX) projZX->Draw("COLZ");
350  } else if (strstr(shortName[pageNumber],".PhiPt")!=0) { // 3D Phi Pt proj.
351  TH2D* projZY = (TH2*)hist3D->Project3D("zy");
352  projZY->SetName(histProjName->Data());
353  projZY->SetYTitle("azimuthal angle (rad");
354  projZY->SetXTitle("Pt (GeV)");
355  gStyle->SetOptStat(0);
356  if (projZY) projZY->Draw("COLZ");
357  } else if (strstr(shortName[pageNumber],"XY")!=0) { // Vertex XY
358  TLine* lineZeroX = new TLine(-1., 0., 1., 0.);
359  TLine* lineZero = new TLine(0., -1., 0., 1.);
360  gStyle->SetOptStat(10);
361  hist2D->Draw("COLZ");
362  lineZeroX->Draw();
363  lineZero->Draw();
364  } else if (strstr(shortName[pageNumber],"Dedx")!=0) { // dE/dx
365  gStyle->SetOptStat(10);
366  (TVirtualPad::Pad()) ->SetLogz();
367  hist2D->Draw("COLZ");
368  } else if (strstr(shortName[pageNumber],"_v")!=0) { // v
369  hist2D->SetMaximum(20.);
370  hist2D->SetMinimum(-20.);
371  gStyle->SetOptStat(0);
372  hist2D->Draw("COLZ");
373  } else { // other 2D
374  gStyle->SetOptStat(10);
375  hist2D->Draw("COLZ");
376  }
377  } else if (strstr(shortName[pageNumber],".Eta")!=0) { // 2D Eta projection
378  if (singleGraph) {
379  TH1D* projX = hist2D->ProjectionX(histName->Data(), 0, 9999);
380  } else {
381  TH1D* projX = hist2D->ProjectionX(histName->Data(), 0, 9999);
382  }
383  projX->SetName(histProjName->Data());
384  char* xTitle = hist2D->GetXaxis()->GetTitle();
385  projX->SetXTitle(xTitle);
386  projX->SetYTitle("Counts");
387  gStyle->SetOptStat(0);
388  if (projX) projX->Draw("H");
389  if (!singleGraph) lineZeroY->Draw();
390  } else if (strstr(shortName[pageNumber],".Pt")!=0) { // 2D Pt projection
391  if (singleGraph) {
392  TH1D* projY = hist2D->ProjectionY(histName->Data(), 0, 9999);
393  } else {
394  TH1D* projY = hist2D->ProjectionY(histName->Data(), 0, 9999, "E");
395  }
396  projY->SetName(histProjName->Data());
397  projY->SetXTitle("Pt (GeV)");
398  projY->SetYTitle("Counts");
399  (TVirtualPad::Pad())->SetLogy();
400  gStyle->SetOptStat(0);
401  if (projY) projY->Draw("H");
402  } else if (strstr(shortName[pageNumber],"Corr")!=0) { // azimuthal corr.
403  float norm = (float)(hist->GetNbinsX()) / hist->Integral();
404  cout << " Normalized by: " << norm << endl;
405  hist->Scale(norm); // normalize height to one
406  if (strstr(shortName[pageNumber],"Sub")!=0) {
407  TF1* funcSubCorr = new TF1("SubCorr", SubCorr, 0., twopi/order, 2);
408  funcSubCorr->SetParNames("chi", "har");
409  funcSubCorr->SetParameters(1., order); // initial value
410  funcSubCorr->SetParLimits(1, 1, 1); // har is fixed
411  hist->Fit("SubCorr");
412  delete funcSubCorr;
413  } else {
414  TF1* funcCos2 = new TF1("funcCos2",
415  "1+[0]*2/100*cos([2]*x)+[1]*2/100*cos(([2]+1)*x)", 0., twopi/order);
416  funcCos2->SetParNames("k=1", "k=2", "har");
417  funcCos2->SetParameters(0, 0, order); // initial values
418  funcCos2->SetParLimits(2, 1, 1); // har is fixed
419  hist->Fit("funcCos2");
420  delete funcCos2;
421  }
422  if (strstr(shortName[pageNumber],"Phi")!=0)
423  hist->SetMinimum(0.9*(hist->GetMinimum()));
424  gStyle->SetOptStat(10);
425  gStyle->SetOptFit(111);
426  hist->Draw("E1");
427  } else if (strstr(shortName[pageNumber],"_q")!=0) { // q distibution
428  gStyle->SetOptStat(110);
429  gStyle->SetOptFit(111);
430  double area = hist->Integral() * qMax / (float)n_qBins;
431  TString* histMulName = new TString("Flow_Mul_Sel");
432  histMulName->Append(*sel);
433  histMulName->Append("_Har");
434  histMulName->Append(*har);
435  TH1* histMult = (TH1*)(histFile[fileN]->Get(histMulName->Data()));
436  if (!histMult) {
437  cout << "### Can't find histogram " << histMulName->Data() << endl;
438  return can;
439  }
440  delete histMulName;
441  float mult = histMult->GetMean();
442  TF1* fit_q = new TF1("qDist", qDist, 0., qMax, 4);
443  fit_q->SetParNames("v", "mult", "area", "g");
444  float qMean = hist->GetMean();
445  float vGuess = (qMean > 1.) ? sqrt(2.*(qMean - 1.) / mult) : 0.03;
446  // the 0.03 is a wild guess
447  vGuess *= 100.;
448  cout << "vGuess = " << vGuess << endl;
449  fit_q->SetParameters(vGuess, mult, area, 0.3); // initial values
450  fit_q->SetParLimits(1, 1, 1); // mult is fixed
451  fit_q->SetParLimits(2, 1, 1); // area is fixed
452  //fit_q->FixParameter(3, 0.6); // g is fixed
453  hist->Fit("qDist");
454  fit_q->Draw("same");
455  } else if (strstr(shortName[pageNumber],"Phi")!=0) { // Phi distibutions
456  hist->SetMinimum(0.9*(hist->GetMinimum()));
457  if (strstr(shortName[pageNumber],"Weight")!=0) {
458  TLine* lineOnePhi = new TLine(0., 1., phiMax, 1.);
459  gStyle->SetOptStat(0);
460  hist->Draw();
461  lineOnePhi->Draw();
462  } else {
463  gStyle->SetOptStat(10);
464  hist->Draw();
465  }
466  } else if (strstr(shortName[pageNumber],"Psi")!=0) { // Psi distibutions
467  gStyle->SetOptStat(10);
468  hist->Draw("E1");
469  } else if (strstr(shortName[pageNumber],"Eta")!=0) { // Eta distibutions
470  if (strstr(shortName[pageNumber],"_v")!=0 ) {
471  hist->SetMaximum(10.);
472  hist->SetMinimum(-10.);
473  }
474  gStyle->SetOptStat(100110);
475  hist->Draw();
476  lineZeroY->Draw();
477  lineYcm->Draw();
478  } else if (strstr(shortName[pageNumber],"Pt")!=0) { // Pt distibutions
479  if (strstr(shortName[pageNumber],"_v")!=0 ) {
480  hist->SetMaximum(30.);
481  hist->SetMinimum(-5.);
482  }
483  gStyle->SetOptStat(100110);
484  hist->Draw();
485  if (strstr(shortName[pageNumber],"v")!=0) {
486  TLine* lineZeroPt = new TLine(0., 0., ptMax, 0.);
487  lineZeroPt->Draw();
488  }
489  } else if (strstr(shortName[pageNumber],"Bin")!=0) { // Bin hists
490  if (strstr(shortName[pageNumber],"Pt")!=0) {
491  TLine* lineDiagonal = new TLine(0., 0., ptMax, ptMax);
492  } else {
493  TLine* lineDiagonal = new TLine(-etaMax, -etaMax, etaMax, etaMax);
494  }
495  gStyle->SetOptStat(0);
496  hist->SetMarkerStyle(21);
497  hist->SetMarkerColor(2);
498  hist->Draw();
499  lineDiagonal->Draw();
500 // } else if (strstr(shortName[pageNumber],"CosPhi")!=0) { // CosPhiLab
501 // TLine* lineZeroHar = new TLine(0.5, 0., 3.5, 0.);
502 // gStyle->SetOptStat(0);
503 // hist->Draw();
504 // lineZeroHar->Draw();
505  } else if (strstr(shortName[pageNumber],"Mult")!=0) { // Mult
506  float mult = hist->GetMean();
507  cout << centr << ": " << mult << endl;
508  (TVirtualPad::Pad())->SetLogy();
509  gStyle->SetOptStat(0);
510  hist->Draw();
511  } else if (strstr(shortName[pageNumber],"MultPart")!=0) { // Part Mult
512  float multPart = hist->GetMean();
513  cout << centr << ": " << multPart << endl;
514  (TVirtualPad::Pad())->SetLogy();
515  gStyle->SetOptStat(0);
516  hist->Draw();
517  } else if (strstr(shortName[pageNumber],"PidMult")!=0) { // PID Mult
518  (TVirtualPad::Pad())->SetLogy();
519  gStyle->SetOptStat(0);
520  hist->Draw();
521  } else if (strstr(shortName[pageNumber],"LYZ_G")!=0) { // LYZ G
522  (TVirtualPad::Pad())->SetLogy();
523  gStyle->SetOptStat(0);
524  TString hist_r0Name("FlowLYZ_r0theta_Sel"); // get r0
525  hist_r0Name += selN;
526  hist_r0Name += "_Har";
527  hist_r0Name += harN;
528  cout << hist_r0Name << endl;
529  TH1D* hist_r0 = (TH1D*)histFile[fileN]->Get(hist_r0Name);
530  float r0 = hist_r0->GetBinContent(1);
531  hist->SetAxisRange(0., 2*r0, "X");
532  hist->Draw();
533  TLine* r0Line = new TLine(r0, 0., r0, 1.);
534  r0Line->SetLineColor(kBlue);
535  r0Line->Draw("same");
536  } else if (strstr(shortName[pageNumber],"LYZ_M")!=0) { // LYZ mult
537  hist->SetMinimum(0.);
538  gStyle->SetOptStat(0);
539  hist->Draw();
540  Float_t mult = hist->GetMean();
541  TString* multChar = new TString("mult= ");
542  *multChar += (int)mult;
543  TLatex l;
544  l.SetNDC();
545  l.SetTextSize(0.1);
546  l.DrawLatex(0.65,0.8,multChar->Data());
547  } else if (strstr(shortName[pageNumber],"LYZ")!=0) { // LYZ
548  hist->SetMinimum(0.);
549  hist->Draw();
550  } else if (strstr(shortName[pageNumber],"_v")!=0 ) { // v 1D
551  TLine* lineZeroHar = new TLine(0.5, 0., 4.5, 0.);
552  hist->SetMaximum(10.);
553  gStyle->SetOptStat(0);
554  hist->Draw();
555  lineZeroHar->Draw();
556  for (int n=1; n <= 4; n++) {
557  v = hist->GetBinContent(n); // output v values
558  err = hist->GetBinError(n);
559  if (n==2) cout << " v2 = " << setprecision(3) << v << " +/- " <<
560  setprecision(2) << err << endl;
561  if (n==4) cout << " v4 = " << setprecision(3) << v << " +/- " <<
562  setprecision(2) << err << endl;
563  if (TMath::IsNaN(v)) {
564  hist->SetBinContent(n, 0.);
565  hist->SetBinError(n, 0.);
566  }
567  }
568  } else if (strstr(shortName[pageNumber],"_Res")!=0 ) { // res
569  for (int n=1; n < 4; n++) {
570  double res = hist->GetBinContent(n); // output res values
571  err = hist->GetBinError(n);
572  if (n==2) cout << " res = " << setprecision(3) << res << " +/- " <<
573  setprecision(2) << err << endl;
574  if (TMath::IsNaN(v)) {
575  hist->SetBinContent(n, 0.);
576  hist->SetBinError(n, 0.);
577  }
578  }
579  hist->Draw();
580  } else { // all other 1D
581  gStyle->SetOptStat(100110);
582  hist->Draw();
583  }
584  }
585 
586  delete [] temp;
587  delete histName;
588  if (histProjName) delete histProjName;
589  for (int m = 0; m < nNames; m++) {
590  delete [] shortName[m];
591  }
592  delete shortName[];
593 
594  return can;
595 }
596 
597 void plotCenAll(Int_t nNames, Int_t selN, Int_t harN, Int_t first = 1) {
598  for (int i = first; i < nNames + 1; i++) {
599  can = plotCen(i, selN, harN);
600  can->Update();
601  cout << "save? y/[n], quit? q" << endl;
602  fgets(tmp, sizeof(tmp), stdin);
603  if (strstr(tmp,"y")!=0) can->Print(".pdf");
604  else if (strstr(tmp,"q")!=0) return;
605  }
606  cout << " Done" << endl;
607 }
608 
609 //-----------------------------------------------------------------------
610 
611 static Double_t qDist(double* q, double* par) {
612  // Calculates the q distribution given the parameters v, mult, area, g
613 
614  double sig2 = 0.5 * (1. + par[3]);
615  double expo = (par[1]*par[0]*par[0]/10000. + q[0]*q[0]) / (2*sig2);
616  Double_t dNdq = par[2] * (q[0]*exp(-expo)/sig2) *
617  TMath::BesselI0(q[0]*par[0]/100.*sqrt(par[1])/sig2);
618 
619  return dNdq;
620 }
621 
622 //-----------------------------------------------------------------------
623 
624 static Double_t SubCorr(double* x, double* par) {
625  // Calculates the n(Psi_a - Psi_b) distribution by fitting chi
626  // From J.-Y. Ollitrault, Nucl. Phys. A590, 561c (1995), Eq. 6. with correc.
627  // The Struve functions are included.
628 
629  double chi2 = par[0] * par[0] / 2; // divide by two for SV chi
630  double z = chi2 * cos(par[1]*x[0]);
631  double TwoOverPi = 2./TMath::Pi();
632 
633  Double_t dNdPsi = exp(-chi2)/TwoOverPi * (TwoOverPi*(1.+chi2)
634  + z*(TMath::BesselI0(z) + TMath::StruveL0(z))
635  + chi2*(TMath::BesselI1(z) + TMath::StruveL1(z)));
636 
637  return dNdPsi;
638 }
639 
641 //
642 // $Log: plotCen.C,v $
643 // Revision 1.29 2010/09/30 19:28:25 posk
644 // Instead of reversing the weight for negative pseudrapidity for odd harmonics,
645 // it is now done only for the first harmonic.
646 // Recentering is now done for all harmonics.
647 //
648 // Revision 1.28 2010/03/05 17:04:40 posk
649 // ROOT 5.22 compatable.
650 // Moved doFlowEvents.C here from StRoot/macros/analysis/
651 //
652 // Revision 1.26 2006/03/22 22:02:12 posk
653 // Updates to macros.
654 //
655 // Revision 1.25 2006/02/22 19:35:20 posk
656 // Added graphs for the StFlowLeeYangZerosMaker
657 //
658 // Revision 1.24 2005/08/26 19:00:26 posk
659 // plot style back to bold
660 //
661 // Revision 1.23 2005/08/05 20:13:44 posk
662 // Improved first guess for qDist fit.
663 //
664 // Revision 1.22 2005/02/08 22:37:55 posk
665 // Fixed trigger histogram for year=4.
666 //
667 // Revision 1.21 2004/11/19 16:54:41 posk
668 // Replaced gPad with (TVirtualPad::Pad()). Reverted to TMath::Struve functions.
669 //
670 // Revision 1.20 2004/11/11 18:25:55 posk
671 // Minor updates.
672 //
673 // Revision 1.19 2004/03/11 18:00:06 posk
674 // Added Random Subs analysis method.
675 //
676 // Revision 1.18 2004/03/01 22:43:44 posk
677 // Changed some "->" to ".".
678 //
679 // Revision 1.17 2003/06/27 21:25:45 posk
680 // v4 and v6 are with repect to the 2nd harmonic event plane.
681 //
682 // Revision 1.16 2003/05/06 21:33:08 posk
683 // Removed some histograms.
684 //
685 // Revision 1.15 2003/03/18 17:58:38 posk
686 // Kirill Fillimonov's improved fit to the angle between subevent planes.
687 //
688 // Revision 1.14 2003/03/17 20:46:57 posk
689 // Improved fit to q dist.
690 //
691 // Revision 1.13 2003/03/11 23:04:31 posk
692 // Includes scalar product hists.
693 //
694 // Revision 1.12 2003/02/25 19:25:33 posk
695 // Improved plotting.
696 //
697 // Revision 1.11 2003/01/10 16:40:55 oldi
698 // Several changes to comply with FTPC tracks:
699 // - Switch to include/exclude FTPC tracks introduced.
700 // The same switch changes the range of the eta histograms.
701 // - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
702 // - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
703 // West, FarWest (depending on vertex.z()).
704 // - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
705 // - Cut to exclude mu-events with no primary vertex introduced.
706 // (This is possible for UPC events and FTPC tracks.)
707 // - Global DCA cut for FTPC tracks added.
708 // - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
709 // - Charge cut for FTPC tracks added.
710 //
711 // Revision 1.10 2001/12/11 22:04:17 posk
712 // Four sets of phiWgt histograms.
713 // StFlowMaker StFlowEvent::PhiWeight() changes.
714 // Cumulant histogram names changed.
715 //
716 // Revision 1.9 2001/11/09 21:15:11 posk
717 // Switched from CERNLIB to TMath. Using global dca instead of dca.
718 //
719 // Revision 1.8 2001/05/22 20:11:24 posk
720 // Changed dEdx graphs.
721 //
722 // Revision 1.7 2000/12/12 15:01:11 posk
723 // Put log comments at end of file.
724 //
725 // Revision 1.6 2000/12/08 17:04:09 oldi
726 // Phi weights for both FTPCs included.
727 //
728 // Revision 1.3 2000/09/15 22:52:56 posk
729 // Added Pt weighting for event plane calculation.
730 //
731 // Revision 1.1 2000/08/31 18:50:33 posk
732 // Added plotCen.C to plot from a series of files with different centralities.
733 //