StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
plot.C
1 //
3 // $Id: plot.C,v 1.70 2011/07/25 15:54:53 posk Exp $
4 //
5 // Author: Art Poskanzer, LBNL, Aug 1999
6 // FTPC added by Markus Oldenburg, MPI, Dec 2000
7 // Description: Macro to plot histograms made by StFlowAnalysisMaker,
8 // StFlowLeeYangZerosMaker, and StFlowScalarProdMaker.
9 // If selN = 0 plot all selections and harmonics.
10 // First time type .x plot.C() to see the menu.
11 // Run Number appended to "ana" is entered in the bottom, left box.
12 // Hist file is anaXX.root where XX is the number.
13 // Default hist file is flow.hist.root .
14 // After the first execution, just type plot(N) .
15 // A negative N plots all pages starting with page N.
16 // Place a symbolic link to this file in StRoot/macros/analysis .
17 //
19 
20 #include <math.h>
21 #include "TMath.h"
22 #include <iostream.h>
23 #include <iomanip.h>
24 
25 
26 const Int_t nHars = 4; // 4
27 const Int_t nSels = 2;
28 const Int_t nSubs = 2;
29 Int_t runNumber = 0;
30 char runName[6];
31 char fileNumber[4] = "x";
32 char fileName[30];
33 TFile* histFile;
34 char tmp[10];
35 char runNo[10];
36 TCanvas* can;
37 
38 TCanvas* plot(Int_t pageNumber=0, Int_t selN=0, Int_t harN=0){
39  gInterpreter->ProcessLine(".O0");
40 
41  Bool_t includeFtpc = kTRUE;
42  Bool_t reCent = kFALSE;
43 
44  bool multiGraph = kFALSE; // set flags
45  bool singleGraph = kFALSE;
46  if (selN == 0) multiGraph = kTRUE;
47 
48  TCanvas* cOld = (TCanvas*)gROOT->GetListOfCanvases(); // delete old canvas
49  if (cOld) cOld->Delete();
50 
51  //gROOT->SetStyle("Pub"); // set style
52  gROOT->SetStyle("Bold"); // set style
53  gROOT->ForceStyle();
54  gStyle->SetOptStat(10);
55  gStyle->SetPalette(1);
56 
57  // names of histograms made by StFlowAnalysisMaker
58  // also projections of some of these histograms
59  const char* baseName1[] = { // with FTPCs
60  "Flow_Res_Sel",
61  "Flow_Trigger",
62  "Flow_VertexZ",
63  "Flow_VertexXY2D",
64  "Flow_EtaSymVerZ2D_Tpc",
65  "Flow_EtaSym_Tpc",
66  "Flow_EtaSymVerZ2D_Ftpc",
67  "Flow_EtaSym_Ftpc",
68  "Flow_CTBvsZDC2D",
69  "Flow_Cent",
70  "Flow_Mult",
71  "Flow_OrigMult",
72  "Flow_MultOverOrig",
73  "Flow_MultEta",
74  "Flow_MultPart",
75  "Flow_Charge_Ftpc",
76  "Flow_DcaGlobal_Tpc",
77  "Flow_DcaGlobal_Ftpc",
78  "Flow_Dca_Tpc",
79  "Flow_Dca_Ftpc",
80  "Flow_Chi2_Tpc",
81  "Flow_Chi2_Ftpc",
82  "Flow_FitPts_Tpc",
83  "Flow_MaxPts_Tpc",
84  "Flow_FitOverMax_Tpc",
85  "Flow_FitPts_Ftpc",
86  "Flow_MaxPts_Ftpc",
87  "Flow_FitOverMax_Ftpc",
88  //"Flow_EtaPtPhi3D",
89  "Flow_EtaPtPhi2D.PhiEta",
90  "Flow_EtaPtPhi2D.PhiPt",
91  "Flow_YieldAll2D",
92  "Flow_YieldAll.Eta",
93  "Flow_YieldAll.Pt",
94  "Flow_YieldPart2D",
95  "Flow_YieldPart.Eta",
96  "Flow_YieldPart.Pt",
97  //"Flow_Bin_Eta",
98  //"Flow_Bin_Pt",
99  "Flow_MeanDedxPos2D",
100  "Flow_MeanDedxNeg2D",
101 // "Flow_PidPiPlusPart",
102 // "Flow_PidPiMinusPart",
103 // "Flow_PidProtonPart",
104 // "Flow_PidAntiProtonPart",
105 // "Flow_PidKplusPart",
106 // "Flow_PidKminusPart",
107 // "Flow_PidDeuteronPart",
108 // "Flow_PidAntiDeuteronPart",
109 // "Flow_PidElectronPart",
110 // "Flow_PidPositronPart",
111  "Flow_PidMult",
112  "Flow_Phi_FarEast_Sel", // first multi graph hist
113  "Flow_Phi_Flat_FarEast_Sel",
114  "Flow_Phi_East_Sel",
115  "Flow_Phi_Flat_East_Sel",
116  "Flow_Phi_West_Sel",
117  "Flow_Phi_Flat_West_Sel",
118  "Flow_Phi_FarWest_Sel",
119  "Flow_Phi_Flat_FarWest_Sel",
120  "Flow_Phi_FtpcFarEast_Sel",
121  "Flow_Phi_Flat_FtpcFarEast_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_Phi_FtpcFarWest_Sel",
127  "Flow_Phi_Flat_FtpcFarWest_Sel",
128  "Flow_Mul_Sel",
129  "Flow_Yield2D_Sel",
130  "Flow_Yield.Eta_Sel",
131  "Flow_Yield.Pt_Sel",
132  "Flow_Psi_Subs",
133  "Flow_Psi_Sel",
134  "Flow_PhiLab_Sel",
135  "Flow_Psi_Diff_Sel",
136  "Flow_Psi_Sub_Corr_Sel",
137  "Flow_Psi_Sub_Corr_Diff_Sel",
138  "Flow_Phi_Corr_Sel",
139  "Flow_vObs2D_Sel",
140  "Flow_vObsEta_Sel",
141  "Flow_vObsPt_Sel",
142  "Flow_v2D_Sel",
143  "Flow_vEta_Sel",
144  "Flow_vPt_Sel",
145  "Flow_q_Sel",
146  "FlowLYZ_vEta_Sel",
147  "FlowLYZ_vPt_Sel",
148  "Flow_vObs2D_ScalarProd_Sel",
149  "Flow_v2D_ScalarProd_Sel",
150  "Flow_vEta_ScalarProd_Sel",
151  "Flow_vPt_ScalarProd_Sel"
152  };
153 
154  const char* baseName2[] = { // without FTPCs
155  "Flow_Res_Sel",
156  "Flow_Trigger",
157  "Flow_VertexZ",
158  "Flow_VertexXY2D",
159  "Flow_EtaSymVerZ2D_Tpc",
160  "Flow_EtaSym_Tpc",
161  "Flow_CTBvsZDC2D",
162  "Flow_Cent",
163  "Flow_Mult",
164  "Flow_OrigMult",
165  "Flow_MultOverOrig",
166  "Flow_MultEta",
167  "Flow_MultPart",
168  "Flow_Dca_Tpc",
169  "Flow_DcaGlobal_Tpc",
170  "Flow_Chi2_Tpc",
171  "Flow_FitPts_Tpc",
172  "Flow_MaxPts_Tpc",
173  "Flow_FitOverMax_Tpc",
174  //"Flow_EtaPtPhi3D",
175  "Flow_EtaPtPhi2D.PhiEta",
176  "Flow_EtaPtPhi2D.PhiPt",
177  "Flow_YieldAll2D",
178  "Flow_YieldAll.Eta",
179  "Flow_YieldAll.Pt",
180  "Flow_YieldPart2D",
181  "Flow_YieldPart.Eta",
182  "Flow_YieldPart.Pt",
183  //"Flow_Bin_Eta",
184  //"Flow_Bin_Pt",
185 // "Flow_MeanDedxPos2D",
186 // "Flow_MeanDedxNeg2D",
187 // "Flow_PidPiPlusPart",
188 // "Flow_PidPiMinusPart",
189 // "Flow_PidProtonPart",
190 // "Flow_PidAntiProtonPart",
191 // "Flow_PidKplusPart",
192 // "Flow_PidKminusPart",
193 // "Flow_PidDeuteronPart",
194 // "Flow_PidAntiDeuteronPart",
195 // "Flow_PidElectronPart",
196 // "Flow_PidPositronPart",
197  "Flow_PidMult",
198  "Flow_Phi_FarEast_Sel", // first multi graph hist
199  "Flow_Phi_Flat_FarEast_Sel",
200  "Flow_Phi_East_Sel",
201  "Flow_Phi_Flat_East_Sel",
202  "Flow_Phi_West_Sel",
203  "Flow_Phi_Flat_West_Sel",
204  "Flow_Phi_FarWest_Sel",
205  "Flow_Phi_Flat_FarWest_Sel",
206  "Flow_Mul_Sel",
207  "Flow_Yield2D_Sel",
208  "Flow_Yield.Eta_Sel",
209  "Flow_Yield.Pt_Sel",
210  "Flow_Psi_Subs",
211  "Flow_Psi_Sel",
212  "Flow_PhiLab_Sel",
213  "Flow_Psi_Diff_Sel",
214  "Flow_Psi_Sub_Corr_Sel",
215  "Flow_Psi_Sub_Corr_Diff_Sel",
216  "Flow_Phi_Corr_Sel",
217  "Flow_vObs2D_Sel",
218  "Flow_vObsEta_Sel",
219  "Flow_vObsPt_Sel",
220  "Flow_v2D_Sel",
221  "Flow_vEta_Sel",
222  "Flow_vPt_Sel",
223  "Flow_q_Sel",
224  "FlowLYZ_vEta_Sel",
225  "FlowLYZ_vPt_Sel",
226  "Flow_vObs2D_ScalarProd_Sel",
227  "Flow_v2D_ScalarProd_Sel",
228  "Flow_vEta_ScalarProd_Sel",
229  "Flow_vPt_ScalarProd_Sel"
230  };
231 
232  const char* baseName3[] = { // for reCent
233  "Flow_Res_Sel",
234  "Flow_Trigger",
235  "Flow_VertexZ",
236  "Flow_VertexXY2D",
237  "Flow_EtaSymVerZ2D_Tpc",
238  "Flow_EtaSym_Tpc",
239  "Flow_EtaSymVerZ2D_Ftpc",
240  "Flow_EtaSym_Ftpc",
241  "Flow_CTBvsZDC2D",
242  "Flow_Cent",
243  "Flow_Mult",
244  "Flow_OrigMult",
245  "Flow_MultOverOrig",
246  "Flow_MultEta",
247  "Flow_MultPart",
248  "Flow_Charge_Ftpc",
249  "Flow_DcaGlobal_Tpc",
250  "Flow_DcaGlobal_Ftpc",
251  "Flow_Dca_Tpc",
252  "Flow_Dca_Ftpc",
253  "Flow_Chi2_Tpc",
254  "Flow_Chi2_Ftpc",
255  "Flow_FitPts_Tpc",
256  "Flow_MaxPts_Tpc",
257  "Flow_FitOverMax_Tpc",
258  "Flow_FitPts_Ftpc",
259  "Flow_MaxPts_Ftpc",
260  "Flow_FitOverMax_Ftpc",
261  //"Flow_EtaPtPhi3D",
262  "Flow_EtaPtPhi2D.PhiEta",
263  "Flow_EtaPtPhi2D.PhiPt",
264  "Flow_YieldAll2D",
265  "Flow_YieldAll.Eta",
266  "Flow_YieldAll.Pt",
267  "Flow_YieldPart2D",
268  "Flow_YieldPart.Eta",
269  "Flow_YieldPart.Pt",
270  //"Flow_Bin_Eta",
271  //"Flow_Bin_Pt",
272  "Flow_MeanDedxPos2D",
273  "Flow_MeanDedxNeg2D",
274 // "Flow_PidPiPlusPart",
275 // "Flow_PidPiMinusPart",
276 // "Flow_PidProtonPart",
277 // "Flow_PidAntiProtonPart",
278 // "Flow_PidKplusPart",
279 // "Flow_PidKminusPart",
280 // "Flow_PidDeuteronPart",
281 // "Flow_PidAntiDeuteronPart",
282 // "Flow_PidElectronPart",
283 // "Flow_PidPositronPart",
284  "Flow_PidMult",
285  "Flow_Mul_Sel", // first multi graph hist
286  "Flow_Yield2D_Sel",
287  "Flow_Yield.Eta_Sel",
288  "Flow_Yield.Pt_Sel",
289  "Flow_Psi_Subs",
290  "Flow_Psi_Sel",
291  "Flow_PhiLab_Sel",
292  "Flow_Psi_Diff_Sel",
293  "Flow_Psi_Sub_Corr_Sel",
294  "Flow_Psi_Sub_Corr_Diff_Sel",
295  "Flow_Phi_Corr_Sel",
296  "Flow_QXY2D_Sel",
297  "Flow_QFTPCSubXY2D_Sel",
298  "Flow_QTPCSubXY2D_Sel",
299  "Flow_vObs2D_Sel",
300  "Flow_vObsEta_Sel",
301  "Flow_vObsPt_Sel",
302  "Flow_v2D_Sel",
303  "Flow_vEta_Sel",
304  "Flow_vPt_Sel",
305  "Flow_q_Sel",
306  "FlowLYZ_vEta_Sel",
307  "FlowLYZ_vPt_Sel",
308  "Flow_vObs2D_ScalarProd_Sel",
309  "Flow_v2D_ScalarProd_Sel",
310  "Flow_vEta_ScalarProd_Sel",
311  "Flow_vPt_ScalarProd_Sel"
312  };
313 
314  const Int_t firstMultFtpc = 38; // these must be set by hand
315  const Int_t firstMultTpc = 27;
316  Int_t nSingles;
317  Int_t nName;
318  if (reCent) {
319  nName = sizeof(baseName3) / sizeof(char*);
320  nSingles = firstMultFtpc + 1;
321  } else if (includeFtpc) {
322  nName = sizeof(baseName1) / sizeof(char*);
323  nSingles = firstMultFtpc + 1;
324  } else {
325  nName = sizeof(baseName2) / sizeof(char*);
326  nSingles = firstMultTpc + 1;
327  }
328  const Int_t nNames = nName;
329 
330  // construct arrays of base and short names
331  char* baseName[nNames];
332  char* shortName[nNames];
333  for (int n = 0; n < nNames; n++) {
334  baseName[n] = new char[35];
335  shortName[n] = new char[35];
336  float etaMax;
337  if (reCent) {
338  strcpy(baseName[n], baseName3[n]);
339  etaMax = 4.5;
340  } else if (includeFtpc) {
341  strcpy(baseName[n], baseName1[n]);
342  etaMax = 4.5;
343  } else {
344  strcpy(baseName[n], baseName2[n]);
345  etaMax = 1.5;
346  }
347  strcpy(shortName[n], baseName[n]);
348  char* cp = strstr(shortName[n],"_Sel");
349  if (cp) *cp = '\0'; // truncate
350  }
351 
352  // input the run number
353  if (runNumber == 0) {
354  cout << " run number? ";
355  fgets(runNo, sizeof(runNo), stdin);
356  runNumber = atoi(runNo);
357  sprintf(runName, "ana%2d", runNumber); // add ana prefix
358  cout << " run name = " << runName << endl;
359  }
360 
361  // input the file number (0 opens flow.hist.root)
362  if (strstr(fileNumber, "x")!=0) {
363  cout << " anaXX.root file number? [0= flow.hist.root] ";
364  fgets(fileNumber, sizeof(fileNumber), stdin);
365  fileNumber[strlen(fileNumber)-1] = '\0';
366  if (strlen(fileNumber) == 1 && strstr(fileNumber,"0")) {
367  sprintf(fileName, "flow.hist.root");
368  } else {
369  sprintf(fileName, "ana%s.root", fileNumber); // insert
370  }
371  cout << " file name = " << fileName << endl;
372  histFile = new TFile(fileName);
373  }
374 
375  // input the page number
376  while (pageNumber <= 1 || pageNumber > nNames) {
377  if (pageNumber < 0) { // plot all
378  plotAll(nNames, selN, harN, -pageNumber);
379  return can;
380  }
381  if (pageNumber == 1) { // plot resolution
382  can = plotResolution();
383  return can;
384  }
385  cout << "-1: \t All" << endl; // print menu
386  for (int i = 0; i < nNames; i++) {
387  cout << i+1 << ":\t " << shortName[i] << endl;
388  }
389  cout << " page number? ";
390  fgets(tmp, sizeof(tmp), stdin);
391  pageNumber = atoi(tmp);
392  }
393  if (pageNumber > 0 && pageNumber <= nSingles) { // plot singles
394  singleGraph = kTRUE;
395  multiGraph = kFALSE;
396  }
397  pageNumber--;
398  cout << pageNumber+1 << ": graph name= " << shortName[pageNumber] << endl;
399 
400  // set constants
401  float twopi = 2. * TMath::Pi();
402  float phiMax = twopi;
403  float Ycm = 0.0;
404  TString* histProjName = NULL;
405 
406  // set row and column numbers
407  char* cp = strstr(shortName[pageNumber],"Subs");
408  int columns = (cp) ? nSubs * nSels : nSels;
409  int rows;
410  rows = (strstr(shortName[pageNumber],"Phi_") ||
411  strstr(shortName[pageNumber],"Psi_Diff") ? 2 : nHars);
412  //rows = (strstr(shortName[pageNumber],"Phi_") ? 2 : nHars);
413  if (strcmp(shortName[pageNumber],"Flow_Phi_Corr")==0) rows = nHars;
414  if (strcmp(shortName[pageNumber],"Flow_Psi_Sub_Corr_Diff")==0) rows = 3;
415  int pads = rows*columns;
416  cout << "pads = " << pads << endl;
417 
418  // make the graph page
419  if (multiGraph) {
420  int canvasWidth = 600, canvasHeight = 780; // portrait
421  } else {
422  int canvasWidth = 780, canvasHeight = 600; // landscape
423  }
424  TString* canName = new TString(shortName[pageNumber]);
425  if (runNumber) {
426  *canName += runNumber;
427  }
428  cout << canName->Data() << endl;
429  can = new TCanvas(canName->Data(), canName->Data(), canvasWidth, canvasHeight);
430  can->ToggleEventStatus();
431  if (multiGraph) {
432  TPaveLabel* title = new TPaveLabel(0.1,0.96,0.9,0.99,shortName[pageNumber]);
433  title->Draw();
434  }
435  TPaveLabel* run = new TPaveLabel(0.1,0.01,0.2,0.03,runName);
436  run->Draw();
437  TDatime now;
438  TPaveLabel* date = new TPaveLabel(0.7,0.01,0.9,0.03,now.AsString());
439  date->Draw();
440  TPad* graphPad = new TPad("Graphs","Graphs",0.01,0.05,0.97,0.95);
441  graphPad->Draw();
442  graphPad->cd();
443 
444  if (multiGraph) { // many graphs on one page
445  graphPad->Divide(columns,rows);
446  int firstK = 0, firstJ = 0, lastK = columns, lastJ = rows;
447  } else if (singleGraph) { // single graph on a page
448  int firstK = 0, firstJ = 0, lastK = 1, lastJ = 1;
449  } else { // one graph from a multi graph page
450  int firstK = selN -1, firstJ = harN -1, lastK = selN, lastJ = harN;
451  }
452  TLine* lineZeroY = new TLine(-etaMax, 0., etaMax, 0.);
453  TLine* lineYcm = new TLine(Ycm, -10., Ycm, 10.);
454  TLine* lineOnePhi = new TLine(0., 1., phiMax, 1.);
455  for (int j = firstJ; j < lastJ; j++) {
456  float order = (float)(j+1);
457  for (int k = firstK ; k < lastK; k++) {
458  int padN = j*columns + k + 1; // pad number
459 
460  // construct histName and histProjName
461  char* temp = new char[35];
462  strcpy(temp,shortName[pageNumber]);
463  char* cproj = strstr(temp,".");
464  if (cproj) { // a projection
465  *cproj = '\0'; // remove from "." on
466  if (singleGraph) {
467  cproj = strstr(temp,"2");
468  if (cproj) { // a 2D projection
469  *cproj = '\0'; // remove from "2D" on
470  strcat(temp,"3D");
471  } else {
472  strcat(temp,"2D");
473  }
474  } else {
475  strcat(temp,"2D_Sel");
476  }
477  TString* histName = new TString(temp);
478  histProjName = new TString(baseName[pageNumber]);
479  if (!singleGraph) {
480  *histProjName += k+1;
481  histProjName->Append("_Har");
482  *histProjName += j+1;
483  }
484  } else { // not projection
485  TString* histName = new TString(baseName[pageNumber]);
486  }
487  if (!singleGraph) {
488  *histName += k+1;
489  histName->Append("_Har");
490  *histName += j+1;
491  }
492  cout << " col= " << k+1 << " row= " << order << " pad= " << padN << "\t"
493  << histName->Data() << endl;
494 
495  // get the histogram
496  bool twoD = kFALSE;
497  bool threeD = kFALSE;
498  if (histProjName) {
499  if (strstr(temp,"3D")) { // 2D projection
500  twoD = kTRUE;
501  TH3* hist3D = dynamic_cast<TH3*>(histFile->Get(histName->Data()));
502  if (!hist3D) {
503  cout << "### Can't find histogram " << histName->Data() << endl;
504  return can;
505  }
506  } else { // 1D projection
507  TH2* hist2D = dynamic_cast<TH2*>(histFile->Get(histName->Data()));
508  if (!hist2D) {
509  cout << "### Can't find histogram " << histName->Data() << endl;
510  return can;
511  }
512  }
513  } else {
514  if (strstr(shortName[pageNumber],"3D")!=0) { // 3D
515  threeD = kTRUE;
516  TH3* hist3D = dynamic_cast<TH3*>(histFile->Get(histName->Data()));
517  if (!hist3D) {
518  cout << "### Can't find histogram " << histName->Data() << endl;
519  return can;
520  }
521  } else if (strstr(shortName[pageNumber],"2D")!=0) { // 2D
522  twoD = kTRUE;
523  TH2* hist2D = dynamic_cast<TH2*>(histFile->Get(histName->Data()));
524  if (!hist2D) {
525  cout << "### Can't find histogram " << histName->Data() << endl;
526  return can;
527  }
528  } else { // 1D
529  TH1* hist = dynamic_cast<TH1*>(histFile->Get(histName->Data()));
530  if (!hist) {
531  cout << "### Can't find histogram " << histName->Data() << endl;
532  return can;
533  }
534  float max = hist->GetXaxis()->GetXmax();
535  }
536  }
537 
538  // make the plots
539  if (multiGraph) graphPad->cd(padN);
540  if (threeD) { // 3D
541  gStyle->SetOptStat(10);
542  hist3D->Draw("BOX");
543  } else if (twoD) { // 2D
544  if (strstr(shortName[pageNumber],".PhiEta")!=0) { // 3D Phi Eta proj.
545  TH2D* projZX = hist3D->Project3D("zx");
546  projZX->SetName(histProjName->Data());
547  projZX->SetYTitle("azimuthal angle (rad)");
548  projZX->SetXTitle("rapidity");
549  gStyle->SetOptStat(0);
550  if (projZX) projZX->Draw("COLZ");
551  } else if (strstr(shortName[pageNumber],".PhiPt")!=0) { // 3D Phi Pt proj.
552  TH2D* projZY = hist3D->Project3D("zy");
553  projZY->SetName(histProjName->Data());
554  projZY->SetYTitle("azimuthal angle (rad");
555  projZY->SetXTitle("Pt (GeV/c)");
556  gStyle->SetOptStat(0);
557  if (projZY) projZY->Draw("COLZ");
558  } else if ((strstr(shortName[pageNumber],"QXY")!=0)) { // Q XY
559  TLine* lineZeroX = new TLine(-0.5, 0., 0.5, 0.);
560  TLine* lineZeroY = new TLine(0., -0.5, 0., 0.5);
561  gStyle->SetOptStat(100110);
562  hist2D->Draw("COLZ");
563  //hist2D->Draw("CONT");
564  lineZeroX->Draw();
565  lineZeroY->Draw();
566  } else if ((strstr(shortName[pageNumber],"TPCSubXY")!=0)) { // QSub XY
567  TLine* lineZeroX = new TLine(-1., 0., 1., 0.);
568  TLine* lineZeroY = new TLine(0., -1., 0., 1.);
569  gStyle->SetOptStat(100110);
570  hist2D->Draw("COLZ");
571  //hist2D->Draw("CONT");
572  lineZeroX->Draw();
573  lineZeroY->Draw();
574  } else if (strstr(shortName[pageNumber],"XY")!=0) { // Vertex XY
575  TLine* lineZeroX = new TLine(-1., 0., 1., 0.);
576  TLine* lineZeroY = new TLine(0., -1., 0., 1.);
577  gStyle->SetOptStat(10);
578  hist2D->Draw("COLZ");
579  lineZeroX->Draw();
580  lineZeroY->Draw();
581  } else if (strstr(shortName[pageNumber],"Dedx")!=0) { // dE/dx
582  TVirtualPad::Pad()->SetLogz();
583  gStyle->SetOptStat(10);
584  hist2D->Draw("COLZ");
585  } else if (strstr(shortName[pageNumber],"_v")!=0) { // v
586  hist2D->SetMaximum(20.);
587  hist2D->SetMinimum(-20.);
588  gStyle->SetOptStat(0);
589  hist2D->Draw("COLZ");
590  } else { // other 2D
591  gStyle->SetOptStat(10);
592  hist2D->Draw("COLZ");
593  }
594  } else if (strstr(shortName[pageNumber],".Eta")!=0) { // 2D Eta projection
595  if (singleGraph) {
596  TH1D* projX = hist2D->ProjectionX(histName->Data());
597  } else {
598  //TH1D* projX = hist2D->ProjectionX(histName->Data(), -1, 9999, "E");
599  TH1D* projX = hist2D->ProjectionX(histName->Data());
600  }
601  projX->SetName(histProjName->Data());
602  char* xTitle = hist2D->GetXaxis()->GetTitle();
603  projX->SetXTitle(xTitle);
604  projX->SetYTitle("Counts");
605  gStyle->SetOptStat(0);
606  if (projX) projX->Draw("H");
607  if (!singleGraph) lineZeroY->Draw();
608  } else if (strstr(shortName[pageNumber],".Pt")!=0) { // 2D Pt projection
609  if (singleGraph) {
610  TH1D* projY = hist2D->ProjectionY(histName->Data());
611  } else {
612  TH1D* projY = hist2D->ProjectionY(histName->Data(), -1, 9999, "E");
613  }
614  projY->SetName(histProjName->Data());
615  projY->SetXTitle("Pt (GeV/c)");
616  projY->SetYTitle("Counts");
617  TVirtualPad::Pad()->SetLogy();
618  gStyle->SetOptStat(0);
619  if (projY) projY->Draw("H");
620  } else if (strstr(shortName[pageNumber],".Phi")!=0) { // 2D Phi projection
621  if (singleGraph) {
622  TH1D* projY = hist2D->ProjectionY(histName->Data());
623  } else {
624  TH1D* projY = hist2D->ProjectionY(histName->Data(), -1, 9999, "E");
625  }
626  projY->SetName(histProjName->Data());
627  projY->SetXTitle("azimuthal angle (rad");
628  projY->SetYTitle("Counts");
629  gStyle->SetOptStat(0);
630  if (projY) projY->Draw("H");
631  } else if (strstr(shortName[pageNumber],"Corr")!=0) { // azimuthal corr.
632  float norm = (hist->Integral()) ? ((float)(hist->GetNbinsX()) / hist->Integral()) : 1.;
633  cout << " Normalized by: " << norm << endl;
634  hist->Scale(norm); // normalize height to one
635  if (strstr(shortName[pageNumber],"Diff")!=0) {
636  TF1* funcCos1 = new TF1("funcCos1",
637  "1+[0]*2/100*cos([1]*x)", 0., twopi);
638  funcCos1->SetParNames("k=1", "har");
639  funcCos1->SetParameters(0, order+1); // initial values
640  funcCos1->SetParLimits(1, 1, 1); // har is fixed
641  //hist->Fit("funcCos1");
642  delete funcCos1;
643  } else if (strstr(shortName[pageNumber],"Sub")!=0) {
644  TF1* funcSubCorr = new TF1("SubCorr", SubCorr, 0., twopi/order, 2);
645  funcSubCorr->SetParNames("chi", "har");
646  funcSubCorr->SetParameters(1., order); // initial value
647  funcSubCorr->SetParLimits(1, 1, 1); // har is fixed
648  hist->Fit("SubCorr");
649  delete funcSubCorr;
650  } else {
651  TF1* funcCos3 = new TF1("funcCos3",
652  "1+[0]*2/100*cos([3]*x)+[1]*2/100*cos(([3]+1)*x)+[2]*2/100*cos(([3]+2)*x)",
653  0., twopi/order);
654  funcCos3->SetParNames("n=har", "n=har+1", "n=har+2", "har");
655  funcCos3->SetParameters(0, 0, 0, order); // initial values
656  funcCos3->SetParLimits(3, 1, 1); // har is fixed
657  //hist->Fit("funcCos3");
658  delete funcCos3;
659  }
660  if (strstr(shortName[pageNumber],"Phi")!=0)
661  hist->SetMinimum(0.9*(hist->GetMinimum()));
662  gStyle->SetOptStat(10);
663  gStyle->SetOptFit(111);
664  hist->Draw("E1");
665  } else if (strstr(shortName[pageNumber],"_q")!=0) { // q distibution
666  gStyle->SetOptStat(110);
667  gStyle->SetOptFit(111);
668  hist->Draw("E1");
669  float n_qBins = (float)hist->GetNbinsX();
670  double area = hist->Integral() * max / n_qBins;
671  TString* histName = new TString("Flow_Mul_Sel");
672  *histName += k+1;
673  histName->Append("_Har");
674  *histName += j+1;
675  TH1* histMult = dynamic_cast<TH1*>(histFile->Get(histName->Data()));
676  if (!histMult) {
677  cout << "### Can't find histogram " << histName->Data() << endl;
678  return can;
679  }
680  delete histName;
681  float mult = histMult->GetMean();
682  TF1* fit_q = new TF1("qDist", qDist, 0., max, 4);
683  fit_q->SetParNames("v", "mult", "area", "g");
684  float qMean = hist->GetMean();
685  float vGuess = (qMean > 1.) ? sqrt(2.*(qMean - 1.) / mult) : 0.03;
686  // the 0.03 is a wild guess
687  vGuess *= 100.;
688  cout << "vGuess = " << vGuess << endl;
689  fit_q->SetParameters(vGuess, mult, area, 0.3); // initial values
690  fit_q->SetParLimits(1, 1, 1); // mult is fixed
691  fit_q->SetParLimits(2, 1, 1); // area is fixed
692  hist->Fit("qDist", "Q");
693  //hist->Fit("qDist");
694  fit_q->Draw();
695  fit_q->FixParameter(3, 0.); // g is fixed
696  fit_q->SetLineStyle(kDotted);
697  hist->Fit("qDist", "Q+");
698  fit_q->Draw("same");
699  fit_q->ReleaseParameter(3); // g is unfixed
700  fit_q->SetParameter(3, 0.5); // initial value
701  fit_q->FixParameter(0, 0.); // v is fixed
702  fit_q->SetLineStyle(kDashed);
703  hist->Fit("qDist", "Q+");
704  fit_q->Draw("same");
705  fit_q->ReleaseParameter(0); // v is unfixed
706  } else if (strstr(shortName[pageNumber],"PhiLab")!=0) { // Phi lab distibution
707  hist->Draw("E1");
708  float norm = (hist->Integral()) ? ((float)(hist->GetNbinsX()) / hist->Integral()) : 1.;
709  cout << " Normalized by: " << norm << endl;
710  hist->Scale(norm); // normalize height to one
711  hist->SetMinimum(0.);
712  hist->SetMaximum(1.3);
713  TF1* funcCosSin = new TF1("funcCosSin",
714  "1.+[0]*2./100.*cos([2]*x)+[1]*2./100.*sin([2]*x)", 0., twopi);
715  funcCosSin->SetParNames("100*cos", "100*sin", "har");
716  funcCosSin->SetParameters(0, 0, order); // initial values
717  funcCosSin->SetParLimits(2, 1, 1); // har is fixed
718  hist->Fit("funcCosSin");
719  delete funcCosSin;
720  gStyle->SetOptFit(111);
721  } else if (strstr(shortName[pageNumber],"Phi")!=0) { // other Phi distibutions
722  //hist->SetMinimum(0.9*(hist->GetMinimum()));
723  hist->SetMinimum(0.);
724  gStyle->SetOptStat(10);
725  hist->Draw();
726  } else if (strstr(shortName[pageNumber],"Psi_Diff")!=0) { // Psi_Diff distibutions
727  hist->Draw("E1");
728  } else if (strstr(shortName[pageNumber],"Psi")!=0) { // Psi distibutions
729  float norm = (hist->Integral()) ? ((float)(hist->GetNbinsX()) / hist->Integral()) : 1.;
730  cout << " Normalized by: " << norm << endl;
731  hist->Scale(norm); // normalize height to one
732  hist->SetMinimum(0.);
733  TF1* funcCosSin = new TF1("funcCosSin",
734  "1.+[0]*2./100.*cos([2]*x)+[1]*2./100.*sin([2]*x)", 0., twopi/order);
735  funcCosSin->SetParNames("100*cos", "100*sin", "har");
736  funcCosSin->SetParameters(0, 0, order); // initial values
737  funcCosSin->SetParLimits(2, 1, 1); // har is fixed
738  hist->Fit("funcCosSin");
739  delete funcCosSin;
740  gStyle->SetOptFit(111);
741  hist->Draw("E1");
742  } else if (strstr(shortName[pageNumber],"Eta")!=0) { // Eta distibutions
743  gStyle->SetOptStat(100110);
744  if (strstr(shortName[pageNumber],"_v")!=0 ) {
745  hist->SetMaximum(10.);
746  hist->SetMinimum(-10.);
747  hist->Draw();
748  lineYcm->Draw();
749  } else {
750  hist->Draw();
751  }
752  lineZeroY->Draw();
753  } else if (strstr(shortName[pageNumber],"Pt")!=0) { // Pt distibutions
754  if (strstr(shortName[pageNumber],"_v")!=0 ) {
755  hist->SetMaximum(25.);
756  hist->SetMinimum(-5.);
757  }
758  gStyle->SetOptStat(100110);
759  hist->Draw();
760  if (strstr(shortName[pageNumber],"v")!=0) {
761  TLine* lineZeroPt = new TLine(0., 0., max, 0.);
762  lineZeroPt->Draw();
763  }
764  } else if (strstr(shortName[pageNumber],"Bin")!=0) { // Bin hists
765  if (strstr(shortName[pageNumber],"Pt")!=0) {
766  TLine* lineDiagonal = new TLine(0., 0., max, max);
767  } else {
768  TLine* lineDiagonal = new TLine(-max, -max, max, max);
769  }
770  gStyle->SetOptStat(0);
771  hist->SetMarkerStyle(21);
772  hist->SetMarkerColor(2);
773  hist->Draw();
774  lineDiagonal->Draw();
775  } else if (strstr(shortName[pageNumber],"PidMult")!=0) { // PID Mult
776  TVirtualPad::Pad()->SetLogy();
777  gStyle->SetOptStat(0);
778  hist->Draw();
779  } else { // all other 1D
780  gStyle->SetOptStat(100110);
781  hist->Draw();
782  }
783  delete [] temp;
784  delete histName;
785  if (histProjName) delete histProjName;
786  }
787  }
788  delete [] shortName;
789 
790  return can;
791 }
792 
793 // macro for the resolution plot
794 TCanvas* plotResolution() {
795  char* resName[] = {
796  "Flow_Cos_Sel",
797  "Flow_Res_Sel",
798  "Flow_v_Sel"
799 // "Flow_v_ScalarProd_Sel",
800 // "Flow_Cumul_v_Order2_Sel",
801 // "Flow_Cumul_v_Order4_Sel"
802  };
803  int columns = nSels;
804  int rows = sizeof(resName) / sizeof(char*);
805  int pads = rows*columns;
806 
807  // make the graph page
808  can = new TCanvas(resName[1], resName[1], 600, 780);
809  can->ToggleEventStatus();
810  TPaveLabel* run = new TPaveLabel(0.1,0.01,0.2,0.03,runName);
811  run->Draw();
812  TDatime now;
813  TPaveLabel* date = new TPaveLabel(0.7,0.01,0.9,0.03,now.AsString());
814  date->Draw();
815  TLine* lineZeroHar = new TLine(0.5, 0., nHars+0.5, 0.);
816  TPad* graphPad = new TPad("Graphs","Graphs",0.01,0.05,0.97,0.99);
817  graphPad->Draw();
818  graphPad->cd();
819  graphPad->Divide(columns,rows);
820 
821  // make the plots
822  float v;
823  float err;
824  for (int j = 0; j < rows; j++) {
825  int resNumber = j;
826  cout << "resolution name= " << resName[resNumber] << endl;
827  for (int k = 0; k < columns; k++) {
828  int padN = j*columns + k +1;
829  TString* histName = new TString(resName[resNumber]);
830  *histName += k+1;
831  cout << "row= " << j << " col= " << k << " pad= " << padN << "\t"
832  << histName->Data() << endl;
833  TH1* hist = dynamic_cast<TH1*>(histFile->Get(histName->Data()));
834  if (!hist) {
835  cout << "### Can't find histogram " << histName->Data() << endl;
836  return can;
837  }
838  graphPad->cd(padN);
839  gStyle->SetOptStat(0);
840  if (strstr(resName[resNumber],"_v")!=0) {
841  hist->SetMaximum(10.);
842  hist->SetMinimum(-5.);
843  } else {
844  hist->SetMaximum(1.1);
845  hist->SetMinimum(0.);
846  }
847  for (int n=1; n < nHars+1; n++) {
848  v = hist->GetBinContent(n); // output v values
849  err = hist->GetBinError(n);
850  cout << " " << n << ": " << setprecision(3) << v << " +/- " << err << endl;
851  if (TMath::IsNaN(v)) {
852  hist->SetBinContent(n, 0.);
853  hist->SetBinError(n, 0.);
854  }
855  }
856  hist->Draw();
857  lineZeroHar->Draw();
858  delete histName;
859  }
860  }
861 
862  return can;
863 }
864 
865 void plotAll(Int_t nNames, Int_t selN, Int_t harN, Int_t first = 1) {
866  for (int i = first; i < nNames + 1; i++) {
867  can = plot(i, selN, harN);
868  can->Update();
869  cout << "save? y/[n], quit? q" << endl;
870  fgets(tmp, sizeof(tmp), stdin);
871  if (strstr(tmp,"y")!=0) can->Print(".pdf");
872  else if (strstr(tmp,"q")!=0) return;
873  else if (strstr(tmp," ")!=0) continue;
874  }
875  cout << " plotAll Done" << endl;
876 }
877 
878 //-----------------------------------------------------------------------
879 
880 static Double_t qDist(double* q, double* par) {
881  // Calculates the q distribution given the parameters v, mult, area, g
882 
883  double sig2 = 0.5 * (1. + par[3]);
884  double expo = (par[1]*par[0]*par[0]/10000. + q[0]*q[0]) / (2*sig2);
885  Double_t dNdq = par[2] * (q[0]*exp(-expo)/sig2) *
886  TMath::BesselI0(q[0]*par[0]/100.*sqrt(par[1])/sig2);
887 
888  return dNdq;
889 }
890 
891 //-----------------------------------------------------------------------
892 
893 static Double_t SubCorr(double* x, double* par) {
894  // Calculates the n(Psi_a - Psi_b) distribution by fitting chi
895  // From J.-Y. Ollitrault, Nucl. Phys. A590, 561c (1995), Eq. 6. with correc.
896 
897  double chi2 = par[0] * par[0] / 2; // divide by two for SV chi
898  double z = chi2 * cos(par[1]*x[0]);
899  double TwoOverPi = 2./TMath::Pi();
900 
901  Double_t dNdPsi = exp(-chi2)/TwoOverPi * (TwoOverPi*(1.+chi2)
902  + z*(TMath::BesselI0(z) + TMath::StruveL0(z))
903  + chi2*(TMath::BesselI1(z) + TMath::StruveL1(z)));
904 
905  return dNdPsi;
906 }
907 
909 //
910 // $Log: plot.C,v $
911 // Revision 1.70 2011/07/25 15:54:53 posk
912 // Added correction for non-flatness of event plane.
913 //
914 // Revision 1.69 2011/03/10 18:56:39 posk
915 // Added histogram for laboratory azimuthal distribution of particles.
916 //
917 // Revision 1.68 2010/09/30 19:28:23 posk
918 // Instead of reversing the weight for negative pseudrapidity for odd harmonics,
919 // it is now done only for the first harmonic.
920 // Recentering is now done for all harmonics.
921 //
922 // Revision 1.67 2010/03/05 17:04:38 posk
923 // ROOT 5.22 compatable.
924 // Moved doFlowEvents.C here from StRoot/macros/analysis/
925 //
926 // Revision 1.66 2009/11/24 19:29:19 posk
927 // Added reCenter to remove acceptance correlations as an option instead of phiWgt.
928 //
929 // Revision 1.65 2006/02/22 19:35:18 posk
930 // Added graphs for the StFlowLeeYangZerosMaker
931 //
932 // Revision 1.64 2005/08/26 19:00:25 posk
933 // plot style back to bold
934 //
935 // Revision 1.63 2005/08/05 20:13:43 posk
936 // Improved first guess for qDist fit.
937 //
938 // Revision 1.62 2004/12/09 23:47:11 posk
939 // Minor changes in code formatting.
940 // Added hist for TPC primary dca to AnalysisMaker.
941 //
942 // Revision 1.61 2004/12/07 23:10:23 posk
943 // Only odd and even phiWgt hists. If the old phiWgt file contains more than
944 // two harmonics, only the first two are read. Now writes only the first two.
945 //
946 // Revision 1.60 2004/12/02 16:10:55 posk
947 // Added gInterpreter->ProcessLine(".O0");
948 //
949 // Revision 1.59 2004/11/19 16:54:40 posk
950 // Replaced gPad with (TVirtualPad::Pad()). Reverted to TMath::Struve functions.
951 //
952 // Revision 1.58 2004/11/11 18:25:43 posk
953 // Minor updates.
954 //
955 // Revision 1.57 2004/03/11 18:00:05 posk
956 // Added Random Subs analysis method.
957 //
958 // Revision 1.56 2004/03/01 22:43:43 posk
959 // Changed some "->" to ".".
960 //
961 // Revision 1.55 2003/08/26 21:10:13 posk
962 // Calculates v8 if nHars=8.
963 //
964 // Revision 1.54 2003/07/30 22:03:26 oldi
965 // Fit("pol4") taken out (which was in accidentally).
966 //
967 // Revision 1.53 2003/07/07 21:58:21 posk
968 // Made units of momentum GeV/c instead of GeV.
969 //
970 // Revision 1.52 2003/06/27 21:25:44 posk
971 // v4 and v6 are with repect to the 2nd harmonic event plane.
972 //
973 // Revision 1.51 2003/05/06 21:33:07 posk
974 // Removed some histograms.
975 //
976 // Revision 1.50 2003/05/02 21:11:13 posk
977 // Reduced the number of harmonics from 3 to 2.
978 //
979 // Revision 1.49 2003/03/18 17:58:36 posk
980 // Kirill Fillimonov's improved fit to the angle between subevent planes.
981 //
982 // Revision 1.48 2003/03/17 20:46:55 posk
983 // Improved fit to q dist.
984 //
985 // Revision 1.47 2003/03/11 23:04:30 posk
986 // Includes scalar product hists.
987 //
988 // Revision 1.46 2003/02/25 19:25:32 posk
989 // Improved plotting.
990 //
991 // Revision 1.45 2003/02/05 18:52:50 posk
992 // Added Bool_t includeFtpc
993 //
994 // Revision 1.44 2003/01/16 16:02:30 posk
995 // Some plotting changes.
996 //
997 // Revision 1.43 2003/01/10 16:40:53 oldi
998 // Several changes to comply with FTPC tracks:
999 // - Switch to include/exclude FTPC tracks introduced.
1000 // The same switch changes the range of the eta histograms.
1001 // - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
1002 // - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
1003 // West, FarWest (depending on vertex.z()).
1004 // - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
1005 // - Cut to exclude mu-events with no primary vertex introduced.
1006 // (This is possible for UPC events and FTPC tracks.)
1007 // - Global DCA cut for FTPC tracks added.
1008 // - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
1009 // - Charge cut for FTPC tracks added.
1010 //
1011 // Revision 1.42 2002/11/26 22:11:54 posk
1012 // First use of doxygen.
1013 //
1014 // Revision 1.41 2002/06/11 21:54:16 posk
1015 // Kirill's further correction to minBias.C for bins with one count.
1016 //
1017 // Revision 1.40 2002/05/21 18:42:18 posk
1018 // Kirill's correction to minBias.C for bins with one count.
1019 //
1020 // Revision 1.39 2002/02/13 22:31:50 posk
1021 // Pt Weight now also weights Phi Weight. Added Eta Weught, default=FALSE.
1022 //
1023 // Revision 1.38 2002/01/14 23:42:57 posk
1024 // Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
1025 //
1026 // Revision 1.37 2001/12/18 19:27:37 posk
1027 // "proton" and "antiproton" replaced by "pr+" and "pr-".
1028 //
1029 // Revision 1.36 2001/12/11 22:04:13 posk
1030 // Four sets of phiWgt histograms.
1031 // StFlowMaker StFlowEvent::PhiWeight() changes.
1032 // Cumulant histogram names changed.
1033 //
1034 // Revision 1.35 2001/11/13 22:47:35 posk
1035 // Documentation updated. Fit to q function moved to macro.
1036 //
1037 // Revision 1.34 2001/11/09 21:15:04 posk
1038 // Switched from CERNLIB to TMath. Using global dca instead of dca.
1039 //
1040 // Revision 1.33 2001/05/22 20:11:20 posk
1041 // Changed dEdx graphs.
1042 //
1043 // Revision 1.32 2000/12/12 15:01:11 posk
1044 // Put log comments at end of file.
1045 //
1046 // Revision 1.31 2000/12/08 17:04:09 oldi
1047 // Phi weights for both FTPCs included.
1048 //
1049 // Revision 1.27 2000/09/26 20:54:12 posk
1050 // Updated documentation.
1051 //
1052 // Revision 1.26 2000/09/15 22:52:56 posk
1053 // Added Pt weighting for event plane calculation.
1054 //
1055 // Revision 1.25 2000/08/31 18:50:32 posk
1056 // Added plotCen.C to plot from a series of files with different centralities.
1057 //
1058 // Revision 1.24 2000/08/12 20:20:15 posk
1059 // More centrality bins.
1060 //
1061 // Revision 1.23 2000/08/01 21:51:20 posk
1062 // Added doubly integrated v.
1063 //
1064 // Revision 1.22 2000/07/12 17:49:39 posk
1065 // Changed EtaSym plots.
1066 //
1067 // Revision 1.21 2000/06/30 14:51:20 posk
1068 // Using MessageMgr. Added graph for Eta Symmetry vs. Vertex Z.
1069 //
1070 // Revision 1.20 2000/05/26 21:25:23 posk
1071 // Use TProfile2D class and profile projection methods.
1072 // Correction needed for >2 subevents.
1073 //
1074 // Revision 1.18 2000/04/13 22:34:16 posk
1075 // Resolution correction is now made.
1076 //
1077 // Revision 1.17 2000/04/10 18:49:10 posk
1078 // Asks for the histogram file number.
1079 //
1080 // Revision 1.16 2000/03/28 23:25:37 posk
1081 // Allow multiple instances.
1082 //
1083 // Revision 1.15 2000/03/21 00:24:45 posk
1084 // Added GetCVS and changed some plot names.
1085 //
1086 // Revision 1.12 2000/02/18 23:44:54 posk
1087 // Added PID and centrality.
1088 //
1089 // Revision 1.11 2000/02/04 16:26:43 posk
1090 // Added correct calculation of event plane resolution for large flow.
1091 //
1092 // Revision 1.2 1999/10/05 16:54:14 posk
1093 // Added getPhiWeight method for making the event plane isotropic.
1094 //
Generated by   doxygen 1.8.5