StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StGlauberPlotMaker.cxx
1 /******************************************************************************
2  * $Id: StGlauberPlotMaker.cxx,v 1.2 2012/04/25 05:02:16 hmasui Exp $
3  * $Log: StGlauberPlotMaker.cxx,v $
4  * Revision 1.2 2012/04/25 05:02:16 hmasui
5  * 5% increment for centrality bins. Added 3rd harmonic eccentricity
6  *
7  ******************************************************************************/
8 
9 #include <assert.h>
10 #include <fstream>
11 
12 #include "TColor.h"
13 #include "TCanvas.h"
14 #include "TError.h"
15 #include "TGraph.h"
16 #include "TGraphErrors.h"
17 #include "TH1.h"
18 #include "TLegend.h"
19 #include "TLine.h"
20 #include "TMath.h"
21 #include "TStyle.h"
22 
23 #include "StMessMgr.h"
24 
25 #include "StGlauberConstUtilities.h"
26 #include "StGlauberPlotMaker.h"
27 
28 ClassImp(StGlauberPlotMaker)
29 
30  using std::ifstream ;
31  using std::ofstream ;
32  using std::endl ;
33  using std::vector ;
34 
35  UInt_t StGlauberPlotMaker::mCanvasId = 0 ;
36  UInt_t StGlauberPlotMaker::mGraphId = 0 ;
37 
38  //____________________________________________________________________________________________________
39  // Default constructor
40  StGlauberPlotMaker::StGlauberPlotMaker(const TString name)
41 : mName(name)
42 {
43  mGraph.clear() ;
44  mGraphDraw.clear() ;
45  mSystematicError = 0 ;
46 }
47 
48 //____________________________________________________________________________________________________
49 // Default destructor
50 StGlauberPlotMaker::~StGlauberPlotMaker()
51 {
52  mGraph.clear() ;
53  mGraphDraw.clear() ;
54  if(mSystematicError) delete mSystematicError ;
55 }
56 
57 //____________________________________________________________________________________________________
58 Int_t StGlauberPlotMaker::Read(const TString filename, const TString type)
59 {
60  ifstream fin(filename.Data());
61  if(!fin){
62  Error("StGlauberPlotMaker::Read", "can't open %s", filename.Data());
63  return 0 ;
64  }
65 
66  LOG_INFO << Form("StGlauberPlotMaker::Read OPEN %s (%s)", filename.Data(), type.Data()) << endm;
67 
68  TGraphErrors* gall = new TGraphErrors() ;
69  TGraphErrors* gdraw = new TGraphErrors() ;
70  const TString name("g" + mName) ;
71  gall ->SetName(Form("%s_all_%d", name.Data(), mGraphId));
72  gdraw->SetName(Form("%s_draw_%d", name.Data(), mGraphId));
73  gall ->SetTitle(type);
74  gdraw->SetTitle(type);
75  mGraphId++; // increment graph id
76 
77  LOG_INFO << "StGlauberPlotMaker::Read Init graphs: "
78  << gall->GetName() << " and " << gdraw->GetName()
79  << endm;
80 
81  const UInt_t ncent = mNCentrality ; // 16 centrality bins for 0-80% (5% increment in 0-80%)
82  UInt_t centralityId ;
83  Double_t centMin, centMax ;
84  Double_t val, error ;
85 
86  for(UInt_t ic=0; ic<StGlauberConstUtilities::GetCentralityBin(); ic++)
87  {
88  fin >> centralityId >> centMin >> centMax >> val >> error ;
89 
90  const Double_t cent = (centMin + centMax)/2.0 ;
91  if( centralityId < ncent ){
92  gdraw->SetPoint(centralityId, cent, val);
93  gdraw->SetPointError(centralityId, 0.0, error);
94  }
95 
96  LOG_INFO << "StGlauberPlotMaker::Read " << mName << ": centrality (min,max) = ("
97  << centMin << "," << centMax << "), val = "
98  << val << " +/- " << error
99  << endm;
100 
101  gall->SetPoint(centralityId, cent, val);
102  gall->SetPointError(centralityId, 0.0, error);
103  }
104 
106  mGraph.push_back( gall );
107  mGraphDraw.push_back( gdraw );
108 
109  return 1 ;
110 }
111 
112 //____________________________________________________________________________________________________
113 Double_t StGlauberPlotMaker::GetYMinimum() const
114 {
115  Double_t ymin = 0.0 ;
116  if( mName.CompareTo("impactparameter", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
117  else if ( mName.CompareTo("npart", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
118  else if ( mName.CompareTo("ncoll", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
119  else if ( mName.CompareTo("multiplicity", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
120  else if ( mName.CompareTo("arearp", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
121  else if ( mName.CompareTo("areapp", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
122  else if ( mName.CompareTo("eccrp", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
123  else if ( mName.CompareTo("eccrpm", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
124  else if ( mName.CompareTo("eccpp_0", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
125  else if ( mName.CompareTo("eccpp_0_2", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
126  else if ( mName.CompareTo("eccpp_1", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
127  else if ( mName.CompareTo("eccpp_1_2", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
128  else if ( mName.CompareTo("eccppm_0", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
129  else if ( mName.CompareTo("eccppm_0_2", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
130  else if ( mName.CompareTo("eccppm_1", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
131  else if ( mName.CompareTo("eccppm_1_2", TString::kIgnoreCase) == 0 ) ymin = 0.0 ;
132  else{
133  Error("StGlauberPlotMaker::GetYMinimum", Form("Cannot find %s in the list. return 0", mName.Data()));
134  return 0.0 ;
135  }
136 
137  return ymin ;
138 }
139 
140 //____________________________________________________________________________________________________
141 Double_t StGlauberPlotMaker::GetYMaximum() const
142 {
143  Double_t ymax = 0.0 ;
144  if( mName.CompareTo("impactparameter", TString::kIgnoreCase) == 0 ) ymax = 18.0 ;
145  else if ( mName.CompareTo("npart", TString::kIgnoreCase) == 0 ) ymax = 420.0 ;
146  else if ( mName.CompareTo("ncoll", TString::kIgnoreCase) == 0 ) ymax = 1200.0 ;
147  else if ( mName.CompareTo("multiplicity", TString::kIgnoreCase) == 0 ) ymax = 1200.0 ;
148  else if ( mName.CompareTo("arearp", TString::kIgnoreCase) == 0 ) ymax = 42.0 ;
149  else if ( mName.CompareTo("areapp", TString::kIgnoreCase) == 0 ) ymax = 42.0 ;
150  else if ( mName.CompareTo("eccrp", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
151  else if ( mName.CompareTo("eccrpm", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
152  else if ( mName.CompareTo("eccpp_0", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
153  else if ( mName.CompareTo("eccpp_0_2", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
154  else if ( mName.CompareTo("eccpp_1", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
155  else if ( mName.CompareTo("eccpp_1_2", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
156  else if ( mName.CompareTo("eccppm_0", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
157  else if ( mName.CompareTo("eccppm_0_2", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
158  else if ( mName.CompareTo("eccppm_1", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
159  else if ( mName.CompareTo("eccppm_1_2", TString::kIgnoreCase) == 0 ) ymax = 0.98 ;
160  else{
161  Error("StGlauberPlotMaker::GetYMaximum", Form("Cannot find %s in the list. return 0", mName.Data()));
162  return 0.0 ;
163  }
164 
165  return ymax ;
166 }
167 
168 //____________________________________________________________________________________________________
169 TString StGlauberPlotMaker::GetYTitle() const
170 {
171  TString title("");
172 
173  if( mName.CompareTo("impactparameter", TString::kIgnoreCase) == 0 ) title = "#LTb#GT (fm)" ;
174  else if ( mName.CompareTo("npart", TString::kIgnoreCase) == 0 ) title = "N_{part}" ;
175  else if ( mName.CompareTo("ncoll", TString::kIgnoreCase) == 0 ) title = "N_{coll}" ;
176  else if ( mName.CompareTo("multiplicity", TString::kIgnoreCase) == 0 ) title = "Multiplicity" ;
177  else if ( mName.CompareTo("arearp", TString::kIgnoreCase) == 0 ) title = "#LTS_{RP}#GT (fm^{2})" ;
178  else if ( mName.CompareTo("areapp", TString::kIgnoreCase) == 0 ) title = "#LTS_{PP}#GT (fm^{2})" ;
179  else if ( mName.CompareTo("eccrp", TString::kIgnoreCase) == 0 ) title = "#LT#varepsilon_{RP}#GT" ;
180  else if ( mName.CompareTo("eccrpm", TString::kIgnoreCase) == 0 ) title = "#LT#varepsilon_{RP}#GT" ;
181  else if ( mName.CompareTo("eccpp_0", TString::kIgnoreCase) == 0 ) title = "#LT#varepsilon_{part}#GT" ;
182  else if ( mName.CompareTo("eccpp_0_2", TString::kIgnoreCase) == 0 ) title = "#varepsilon_{part}{2}" ;
183  else if ( mName.CompareTo("eccpp_1", TString::kIgnoreCase) == 0 ) title = "#LT#varepsilon_{3,part}#GT" ;
184  else if ( mName.CompareTo("eccpp_1_2", TString::kIgnoreCase) == 0 ) title = "#varepsilon_{3,part}{2}" ;
185  else if ( mName.CompareTo("eccppm_0", TString::kIgnoreCase) == 0 ) title = "#LT#varepsilon_{part}#GT" ;
186  else if ( mName.CompareTo("eccppm_0_2", TString::kIgnoreCase) == 0 ) title = "#varepsilon_{part}{2}" ;
187  else if ( mName.CompareTo("eccppm_1", TString::kIgnoreCase) == 0 ) title = "#LT#varepsilon_{3,part}#GT" ;
188  else if ( mName.CompareTo("eccppm_1_2", TString::kIgnoreCase) == 0 ) title = "#varepsilon_{3,part}{2}" ;
189  else{
190  Error("StGlauberPlotMaker::GetYTitle", Form("Cannot find %s in the list. return 0", mName.Data()));
191  return "";
192  }
193 
194  return title ;
195 }
196 
197 //____________________________________________________________________________________________________
198 TGraphErrors* StGlauberPlotMaker::Divide(const TGraphErrors& g0, const TGraphErrors& g1) const
199 {
200  TGraphErrors* g = new TGraphErrors();
201  g->SetMarkerSize( g0.GetMarkerSize() );
202  g->SetMarkerStyle( g0.GetMarkerStyle() );
203  g->SetMarkerColor( g0.GetMarkerColor() );
204  g->SetLineColor( g0.GetLineColor() );
205 
206  for(Int_t i=0;i<g0.GetN();i++){
207  const Double_t y0 = g0.GetY()[i] ;
208  const Double_t y1 = g1.GetY()[i] ;
209  const Double_t y0e = g0.GetEY()[i] ;
210  const Double_t y1e = g1.GetEY()[i] ;
211 
212  if( y1 != 0.0 ){
213  const Double_t r = y0/y1 ;
214  const Double_t re = TMath::Abs(r) * TMath::Sqrt(TMath::Power(y0e/y0,2.0)+TMath::Power(y1e/y1,2.0)) ;
215  const Int_t n = g->GetN() ;
216 
217  g->SetPoint(n, g0.GetX()[i], r) ;
218  g->SetPointError(n, g0.GetEX()[i], re) ;
219  }
220  }
221 
222  return g ;
223 }
224 
225 //____________________________________________________________________________________________________
226 TGraphErrors* StGlauberPlotMaker::SystematicErrors(const UInt_t mode)
227 {
229  if( mGraph.size() < 2 ){
230  Error("StGlauberPlotMaker::SystematicErrors", "Number of graphs < 2, at least 2 graphs are needed to evaluate sys. errors. abort");
231  assert(0);
232  }
233 
235  LOG_INFO << "StGlauberPlotMaker::SystematicErrors Start evaluating systematic errors" << endm;
236 
239  const UInt_t nCent = mGraph[0]->GetN() ;
240  Double_t qSum[nCent] ;
241  for(UInt_t ic=0; ic<nCent; ic++){
242  qSum[ic] = 0.0 ;
243  }
244 
245  const UInt_t nSize = mGraph.size() ;
246  for(UInt_t id=1; id<nSize; id++){
247  LOG_INFO << "StGlauberPlotMaker::SystematicErrors Use " << mGraph[id]->GetTitle() << endm;
248 
249  for(UInt_t ic=0; ic<nCent; ic++){
250  const Double_t val = mGraph[id]->GetY()[ic];
251  const Double_t ref = mGraph[0]->GetY()[ic];
252  const Double_t diff = val - ref;
253  qSum[ic] += diff * diff ;
254  }
255  }
256 
257  // qSum = sqrt(sum{diff})
258  mSystematicError = new TGraph();
259 
260  Double_t factor = 1.0/TMath::Sqrt(12.0) ;
261  if( mode == 1 ) factor = 1.0 ;
262 
263  for(UInt_t ic=0; ic<nCent; ic++){
264  qSum[ic] = TMath::Sqrt(qSum[ic]) * factor ;
265  mSystematicError->SetPoint(ic, ic+0.5, qSum[ic]);
266  }
267 
268  // Fill graph
269  TGraphErrors* g = new TGraphErrors() ;
270  g->SetFillColor(kYellow);
271  g->SetLineColor(kYellow);
272 
273  for(UInt_t ic=0; ic<mNCentrality; ic++){
274  g->SetPoint(ic, mGraph[0]->GetX()[ic], mGraph[0]->GetY()[ic]);
275  g->SetPointError(ic, 0.0, qSum[ic]);
276  }
277 
278  return g ;
279 }
280 
281 //____________________________________________________________________________________________________
282 void StGlauberPlotMaker::Draw(const UInt_t mode)
283 {
285  TColor::CreateColorWheel();
286  gStyle->SetPadRightMargin(0.05);
287 
289  TGraphErrors* gSysError = SystematicErrors(mode) ;
290 
292  TCanvas* c1 = new TCanvas(Form("c%d", mCanvasId), Form("c%d", mCanvasId++), 600, 800);
293  c1->Divide(1, 2);
294 
295  // Styles, and colors
296  const UInt_t style[] = {20, 21, 25, 22, 26, 29, 21, 25, 22, 26, 27, 28};
297  const UInt_t color[] = {kBlack, kRed, kRed, kBlue, kBlue, kMagenta+1,
298  kGreen+2, kGreen+2,
299  kOrange+1, kOrange+1,
300  kCyan-3, kCyan-3
301  };
302 
303  //----------------------------------------------------------------------------------------------------
304  // Draw
305  //----------------------------------------------------------------------------------------------------
306  c1->cd(1);
307 
309  const Double_t ymin = GetYMinimum() ;
310  const Double_t ymax = GetYMaximum() ;
311 
312  TH1* frame = c1->GetPad(1)->DrawFrame(0, ymin, 80, ymax);
313  frame->SetXTitle("% Most central");
314  frame->SetYTitle(GetYTitle());
315 
316  TLegend* leg = new TLegend(0.05, 0.1, 0.95, 0.9);
317  leg->SetTextSize(0.05);
318  leg->SetFillColor(10);
319  // leg->SetBorderSize(1);
320 
321  const UInt_t nSize = mGraphDraw.size() ;
322  TGraphErrors* gRatio[nSize] ;
323  for(UInt_t id=0; id<nSize; id++){
324  mGraphDraw[id]->SetMarkerSize(1.2);
325  mGraphDraw[id]->SetMarkerStyle(style[id]);
326  mGraphDraw[id]->SetMarkerColor(color[id]);
327  mGraphDraw[id]->SetLineColor(color[id]);
328  mGraphDraw[id]->Draw("P");
329 
330  // Calculate ratio
331  gRatio[id] = Divide(*mGraphDraw[id], *mGraphDraw[0]) ;
332 
333  // Add entry in legend
334  leg->AddEntry(mGraphDraw[id], gRatio[id]->GetTitle(), "P");
335  }
336  mGraphDraw[0]->Draw("P");
337 
338  //----------------------------------------------------------------------------------------------------
339  // Draw ratio
340  //----------------------------------------------------------------------------------------------------
341  TPad* pad = (TPad*) c1->GetPad(2);
342  pad->Divide(2, 1);
343  pad->cd(1);
344 
345  const Double_t rmin = 1.0 - 0.3 ;
346  const Double_t rmax = 1.0 + 0.3 ;
347  TH1* frameRatio = pad->GetPad(1)->DrawFrame(0, rmin, 80, rmax);
348  frameRatio->SetXTitle("% Most central");
349  frameRatio->SetYTitle(Form("Ratio of %s", GetYTitle().Data()));
350 
351  TLine* lzero = new TLine(0, 1.0, 80, 1.0);
352  lzero->SetLineStyle(3);
353  lzero->Draw() ;
354 
355  for(UInt_t id=1; id<nSize; id++){
356  gRatio[id]->Draw("P");
357  }
358 
359  //----------------------------------------------------------------------------------------------------
360  // Draw legend
361  //----------------------------------------------------------------------------------------------------
362  pad->cd(2);
363  leg->Draw() ;
364 
365  c1->cd();
366  c1->Update();
367 
368  c1->Print(Form("figure/systematicerror_%s.eps", mName.Data()));
369  // c1->Print(Form("%s_systematic_error.png", mName.Data()));
370 
371  //----------------------------------------------------------------------------------------------------
372  // Draw default vs centrality with systematic error (eps)
373  //----------------------------------------------------------------------------------------------------
374  TCanvas* c2 = new TCanvas(Form("c%d", mCanvasId), Form("c%d", mCanvasId++));
375  TH1* frame2 = c2->DrawFrame(0, ymin, 80, ymax);
376  frame2->SetXTitle("% Most central");
377  frame2->SetYTitle(GetYTitle());
378 
379  gSysError->Draw("E3");
380  mGraphDraw[0]->Draw("P");
381 
382  c2->cd();
383  c2->Update();
384 
385  c2->Print(Form("figure/%s_vs_centrality_with_systematicerror.eps", mName.Data()));
386 
387  //----------------------------------------------------------------------------------------------------
388  // Draw default vs centrality with systematic error (png for web, small TCanvas)
389  //----------------------------------------------------------------------------------------------------
390  TCanvas* c3 = new TCanvas(Form("c%d", mCanvasId), Form("c%d", mCanvasId++), 700*0.5, 500*0.5);
391  TH1* frame3 = c3->DrawFrame(0, ymin, 80, ymax);
392  frame3->SetXTitle("% Most central");
393  frame3->SetYTitle(GetYTitle());
394 
395  gSysError->Draw("E3");
396  mGraphDraw[0]->Draw("P");
397 
398  c3->cd();
399  c3->Update();
400 
401  c3->Print(Form("figure/%s_vs_centrality_with_systematicerror.png", mName.Data()));
402 
403  //----------------------------------------------------------------------------------------------------
404  // Write table
405  //----------------------------------------------------------------------------------------------------
406  const TString tableName(Form("table_%s_vs_centrality_systematicerror.txt", mName.Data()));
407  ofstream fout(tableName.Data()) ;
408  LOG_INFO << Form("StGlauberPlotMaker::Draw Write table %s in the current directory", tableName.Data()) << endm;
409 
410  for(Int_t ic=0; ic<mGraph[0]->GetN(); ic++){
411  const Double_t val = mGraph[0]->GetY()[ic] ;
412  const Double_t sys = mSystematicError->GetY()[ic] ;
413  const Double_t centMin = StGlauberConstUtilities::GetCentralityMin(ic) ;
414  const Double_t centMax = StGlauberConstUtilities::GetCentralityMax(ic) ;
415 
416  fout << Form("%10d %1.1f %1.1f %1.5f %1.5f",
417  ic, centMin, centMax, val, sys) << endl;
418  }
419  fout << endl << endl ;
420  fout << "# <centrality bin> <minimum centrality> <maximum centrality> <value> <sys. error>" << endl;
421 }
422 
423 
Error
Definition: FJcore.h:367
StGlauberPlotMaker::Draw
void Draw(const UInt_t mode=0)
Definition: StGlauberPlotMaker.cxx:282
StGlauberPlotMaker::~StGlauberPlotMaker
virtual ~StGlauberPlotMaker()
Default constructor.
Definition: StGlauberPlotMaker.cxx:50
StGlauberPlotMaker::Read
Int_t Read(const TString filename, const TString type)
Default destructor.
Definition: StGlauberPlotMaker.cxx:58
Generated by   doxygen 1.8.5