23 const Int_t nCens = 10;
29 TFile* histFile[nCens];
33 TCanvas* plotCen(Int_t pageNumber=0, Int_t selN=2, Int_t harN=2){
34 gInterpreter->ProcessLine(
".O0");
36 TCanvas* cOld = (TCanvas*)gROOT->GetListOfCanvases();
37 if (cOld) cOld->Delete();
40 gROOT->SetStyle(
"Bold");
41 gStyle->SetPalette(1);
43 gStyle->SetLabelSize(.1,
"X");
45 int canvasWidth = 600, canvasHeight = 780;
48 bool oddPads = (nCens) % 2;
50 rows = nCens/columns + 1;
58 const char* baseName[] = {
65 "Flow_v_ScalarProd_Sel",
66 "Flow_Cumul_v_Order2_Sel",
67 "Flow_Cumul_v_Order4_Sel",
70 "Flow_EtaSymVerZ2D_Tpc",
72 "Flow_EtaSymVerZ2D_Ftpc",
85 "Flow_DcaGlobal_Ftpc",
90 "Flow_FitOverMax_Tpc",
93 "Flow_FitOverMax_Ftpc",
100 "Flow_MeanDedxPos2D",
101 "Flow_MeanDedxNeg2D",
115 "Flow_Yield.Eta_Sel",
128 "Flow_Psi_Sub_Corr_Sel",
129 "Flow_Psi_Sub_Corr_Diff_Sel",
138 "FlowLYZ_r0theta_Sel",
139 "FlowLYZ_Vtheta_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"
152 int nName =
sizeof(baseName) /
sizeof(
char*);
153 const Int_t nNames = nName;
154 const int nDoubles = 9;
155 const int nSingles = 46 + nDoubles;
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");
166 cout <<
"Harmonic = " << harN << endl << endl;
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);
173 cout <<
" first run name = " << runName << endl;
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);
184 while (pageNumber <= 0 || pageNumber > nNames) {
185 if (pageNumber < 0) {
186 plotCenAll(nNames, selN, harN, -pageNumber);
189 cout <<
"-1: \t All" << endl;
190 for (
int i = 0; i < nNames; i++) {
191 cout << i+1 <<
":\t " << shortName[i] << endl;
193 cout <<
" page number? ";
194 fgets(tmp,
sizeof(tmp), stdin);
195 pageNumber = atoi(tmp);
199 bool multiGraph = kFALSE;
200 bool doubleGraph = kFALSE;
201 bool singleGraph = kFALSE;
202 if (pageNumber > 0 && pageNumber <= nDoubles) {
204 }
else if (pageNumber > nDoubles && pageNumber <= nSingles) {
212 float twopi = 2. * TMath::Pi();
217 float phiMax = twopi;
220 TString* histProjName = NULL;
224 sprintf(sel,
"%d",selN);
226 sprintf(har,
"%d",harN);
227 float order = (float)harN;
228 char* temp =
new char[35];
229 strcpy(temp,shortName[pageNumber]);
230 char* cproj = strstr(temp,
".");
234 cproj = strstr(temp,
"2");
242 strcat(temp,
"2D_Sel");
244 TString* histName =
new TString(temp);
245 histProjName =
new TString(baseName[pageNumber]);
247 histProjName->Append(*sel);
248 histProjName->Append(
"_Har");
249 histProjName->Append(*har);
252 TString* histName =
new TString(baseName[pageNumber]);
254 if (!singleGraph) histName->Append(*sel);
256 histName->Append(
"_Har");
257 histName->Append(*har);
259 cout << pageNumber+1 <<
": graph name= " << histName->Data() << endl;
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());
267 TPaveLabel* run =
new TPaveLabel(0.1,0.01,0.2,0.03,runName);
270 TPaveLabel* date =
new TPaveLabel(0.7,0.01,0.9,0.03,now.AsString());
272 TPad* graphPad =
new TPad(
"Graphs",
"Graphs",0.01,0.05,0.97,0.95);
275 graphPad->Divide(columns, rows);
276 TLine* lineZeroY =
new TLine(yMin, 0., yMax, 0.);
277 TLine* lineYcm =
new TLine(Ycm, -10., Ycm, 10.);
281 for (
int i = 0; i < pads; i++) {
283 int padN = fileN + 1;
284 centr = oddPads ? padN : padN-1;
285 sprintf(histTitle,
"Centrality %d",centr);
286 cout <<
"centrality= " << centr << endl;
292 if (strstr(temp,
"3D")) {
294 TH3* hist3D = (TH3*)(histFile[fileN]->Get(histName->Data()));
296 cout <<
"### Can't find histogram " << histName->Data() << endl;
299 hist3D->SetTitle(histTitle);
301 TH2* hist2D = (TH2*)(histFile[fileN]->Get(histName->Data()));
303 cout <<
"### Can't find histogram " << histName->Data() << endl;
306 hist2D->SetTitle(histTitle);
309 if (strstr(shortName[pageNumber],
"3D")!=0) {
311 TH3* hist3D = (TH3*)(histFile[fileN]->Get(histName->Data()));
313 cout <<
"### Can't find histogram " << histName->Data() << endl;
316 hist3D->SetTitle(histTitle);
317 }
else if (strstr(shortName[pageNumber],
"2D")!=0) {
319 TH2* hist2D = (TH2*)(histFile[fileN]->Get(histName->Data()));
321 cout <<
"### Can't find histogram " << histName->Data() << endl;
324 hist2D->SetTitle(histTitle);
326 TH1* hist = (TH1*)(histFile[fileN]->Get(histName->Data()));
328 cout <<
"### Can't find histogram " << histName->Data() << endl;
331 float ptMax = hist->GetXaxis()->GetXmax();
332 hist->SetTitle(histTitle);
340 gStyle->SetOptStat(10);
343 if (strstr(shortName[pageNumber],
".PhiEta")!=0) {
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) {
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) {
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");
364 }
else if (strstr(shortName[pageNumber],
"Dedx")!=0) {
365 gStyle->SetOptStat(10);
366 (TVirtualPad::Pad()) ->SetLogz();
367 hist2D->Draw(
"COLZ");
368 }
else if (strstr(shortName[pageNumber],
"_v")!=0) {
369 hist2D->SetMaximum(20.);
370 hist2D->SetMinimum(-20.);
371 gStyle->SetOptStat(0);
372 hist2D->Draw(
"COLZ");
374 gStyle->SetOptStat(10);
375 hist2D->Draw(
"COLZ");
377 }
else if (strstr(shortName[pageNumber],
".Eta")!=0) {
379 TH1D* projX = hist2D->ProjectionX(histName->Data(), 0, 9999);
381 TH1D* projX = hist2D->ProjectionX(histName->Data(), 0, 9999);
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) {
392 TH1D* projY = hist2D->ProjectionY(histName->Data(), 0, 9999);
394 TH1D* projY = hist2D->ProjectionY(histName->Data(), 0, 9999,
"E");
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) {
403 float norm = (float)(hist->GetNbinsX()) / hist->Integral();
404 cout <<
" Normalized by: " << norm << endl;
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);
410 funcSubCorr->SetParLimits(1, 1, 1);
411 hist->Fit(
"SubCorr");
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);
418 funcCos2->SetParLimits(2, 1, 1);
419 hist->Fit(
"funcCos2");
422 if (strstr(shortName[pageNumber],
"Phi")!=0)
423 hist->SetMinimum(0.9*(hist->GetMinimum()));
424 gStyle->SetOptStat(10);
425 gStyle->SetOptFit(111);
427 }
else if (strstr(shortName[pageNumber],
"_q")!=0) {
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()));
437 cout <<
"### Can't find histogram " << histMulName->Data() << endl;
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;
448 cout <<
"vGuess = " << vGuess << endl;
449 fit_q->SetParameters(vGuess, mult, area, 0.3);
450 fit_q->SetParLimits(1, 1, 1);
451 fit_q->SetParLimits(2, 1, 1);
455 }
else if (strstr(shortName[pageNumber],
"Phi")!=0) {
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);
463 gStyle->SetOptStat(10);
466 }
else if (strstr(shortName[pageNumber],
"Psi")!=0) {
467 gStyle->SetOptStat(10);
469 }
else if (strstr(shortName[pageNumber],
"Eta")!=0) {
470 if (strstr(shortName[pageNumber],
"_v")!=0 ) {
471 hist->SetMaximum(10.);
472 hist->SetMinimum(-10.);
474 gStyle->SetOptStat(100110);
478 }
else if (strstr(shortName[pageNumber],
"Pt")!=0) {
479 if (strstr(shortName[pageNumber],
"_v")!=0 ) {
480 hist->SetMaximum(30.);
481 hist->SetMinimum(-5.);
483 gStyle->SetOptStat(100110);
485 if (strstr(shortName[pageNumber],
"v")!=0) {
486 TLine* lineZeroPt =
new TLine(0., 0., ptMax, 0.);
489 }
else if (strstr(shortName[pageNumber],
"Bin")!=0) {
490 if (strstr(shortName[pageNumber],
"Pt")!=0) {
491 TLine* lineDiagonal =
new TLine(0., 0., ptMax, ptMax);
493 TLine* lineDiagonal =
new TLine(-etaMax, -etaMax, etaMax, etaMax);
495 gStyle->SetOptStat(0);
496 hist->SetMarkerStyle(21);
497 hist->SetMarkerColor(2);
499 lineDiagonal->Draw();
505 }
else if (strstr(shortName[pageNumber],
"Mult")!=0) {
506 float mult = hist->GetMean();
507 cout << centr <<
": " << mult << endl;
508 (TVirtualPad::Pad())->SetLogy();
509 gStyle->SetOptStat(0);
511 }
else if (strstr(shortName[pageNumber],
"MultPart")!=0) {
512 float multPart = hist->GetMean();
513 cout << centr <<
": " << multPart << endl;
514 (TVirtualPad::Pad())->SetLogy();
515 gStyle->SetOptStat(0);
517 }
else if (strstr(shortName[pageNumber],
"PidMult")!=0) {
518 (TVirtualPad::Pad())->SetLogy();
519 gStyle->SetOptStat(0);
521 }
else if (strstr(shortName[pageNumber],
"LYZ_G")!=0) {
522 (TVirtualPad::Pad())->SetLogy();
523 gStyle->SetOptStat(0);
524 TString hist_r0Name(
"FlowLYZ_r0theta_Sel");
526 hist_r0Name +=
"_Har";
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");
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) {
537 hist->SetMinimum(0.);
538 gStyle->SetOptStat(0);
540 Float_t mult = hist->GetMean();
541 TString* multChar =
new TString(
"mult= ");
542 *multChar += (int)mult;
546 l.DrawLatex(0.65,0.8,multChar->Data());
547 }
else if (strstr(shortName[pageNumber],
"LYZ")!=0) {
548 hist->SetMinimum(0.);
550 }
else if (strstr(shortName[pageNumber],
"_v")!=0 ) {
551 TLine* lineZeroHar =
new TLine(0.5, 0., 4.5, 0.);
552 hist->SetMaximum(10.);
553 gStyle->SetOptStat(0);
556 for (
int n=1; n <= 4; n++) {
557 v = hist->GetBinContent(n);
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.);
568 }
else if (strstr(shortName[pageNumber],
"_Res")!=0 ) {
569 for (
int n=1; n < 4; n++) {
570 double res = hist->GetBinContent(n);
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.);
581 gStyle->SetOptStat(100110);
588 if (histProjName)
delete histProjName;
589 for (
int m = 0; m < nNames; m++) {
590 delete [] shortName[m];
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);
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;
606 cout <<
" Done" << endl;
611 static Double_t qDist(
double* q,
double* par) {
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);
624 static Double_t SubCorr(
double* x,
double* par) {
629 double chi2 = par[0] * par[0] / 2;
630 double z = chi2 * cos(par[1]*x[0]);
631 double TwoOverPi = 2./TMath::Pi();
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)));