StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
postProcess.C
1 #include <stdio.h>
2 #include <string.h>
3 #include "TCanvas.h"
4 #include "StRoot/StEStructPool/Fluctuations/StEStructSigAnal.h"
5 #include "TStyle.h"
6 #include "TSystem.h"
7 #include "TFile.h"
8 #include "TH2D.h"
9 
10 /***********************************************************
11  * NOTE: In exportForMatlab I have typed in centrality and pt range names.
12  * This will likely need to be modified sometime somewhere.
13  ***********************************************************/
14 
15 /*
16 // gROOT->LoadMacro("load2ptLibs.C");
17 // load2ptLibs();
18 // .L postProcess.C+
19  */
20 
21  /*
22 // postProcess *pp = new postProcess("/common/star/stardata/estruct/prindle/Hijing/auau200/QuenchOff/testForJamie");
23 // pp->createSigmas();
24 // pp->centralityMults(1200,1);
25 // pp->multiplicity(1200,1);
26 // pp->deltaSigma2(1);
27  */
28 class postProcess {
29  public:
30  char mDir[2048];
31  StEStructSigAnal *msigA;
32  TFile *mFile;
33  TH2D *mSig[16*11];
34  TCanvas *mCan;
35  char *mCanType;
36  int mnEta, mnPhi;
37  int mnCents, mnPts, mnPtCents, mnx, mny;
38  double iEta, iPhi, sigN, sigPt, error;
39 
40  postProcess( char *dir );
41  virtual ~postProcess();
42  void createSigmas();
43  void plotCentrality( int iCent, int ptBin=0, int print=0 );
44  void plotType( char *pType, int ptBin=0, int print=0 );
45  void centralityMults( int N0, int print=0 );
46  void multiplicity( int N0, int print=0 );
47  void pt( int P0, int print=0 );
48  void deltaSigma2( int print=0 );
49  void exportForMatlab();
50  void importFromMatlab();
51  int getNCents();
52  int getNPts();
53  int getNPtCents();
54  int getNX( int nCent );
55  int getNY( int nCent );
56  // A few methods only useful for Hijing data.
57  void plotMCInfo( int print=0 );
58  void plotBinaryCollisions( int print=0 );
59  void plotParticipants( int print=0 );
60  void exportMCInfo();
61 };
62 
63 postProcess::postProcess( char *dir ) {
64  gStyle->SetPadTopMargin(0.2);
65  gStyle->SetPadLeftMargin(0.15);
66  gStyle->SetCanvasColor(0);
67  gStyle->SetTitleBorderSize(0);
68  gStyle->SetTitleBorderSize(0);
69  gStyle->SetPadBorderSize(0);
70  gStyle->SetTitleColor(0);
71  gStyle->SetTitleFillColor(0);
72  gStyle->SetOptStat(0);
73  gStyle->SetPalette(1);
74  gStyle->SetTitleX(0.2);
75  gStyle->SetTitleH(0.08);
76  gStyle->SetPadBorderMode(0);
77 
78  strcpy( mDir, dir );
79 
80  mnEta = 16;
81  mnPhi = 24;
82  mnCents = 0;
83  mnPts = 0;
84  mnPtCents = 0;
85  mnx = 0;
86  mny = 0;
87  mFile = 0;
88  msigA = 0;
89  mCan = 0;
90  mCanType = 0;
91 }
92 postProcess::~postProcess() {
93  free(mDir);
94  if (mFile) {
95  delete mFile;
96  }
97  if (msigA) {
98  delete msigA;
99  }
100  if (mCan) {
101  delete mCan;
102  }
103 }
104 void postProcess::createSigmas() {
105 
106  char buffer[1024];
107  sprintf(buffer,"%s/data/Data.root",mDir);
108  if (msigA) {
109  delete msigA;
110  }
111  msigA = new StEStructSigAnal(buffer);
112  msigA->fillHistograms();
113  sprintf(buffer,"%s/data/Sigmas.root",mDir);
114  TFile *out = new TFile(buffer,"RECREATE");
115  msigA->writeHistograms(out);
116  delete out;
117 }
118 void postProcess::plotCentrality( int iCent, int ptBin, int print ) {
119  if (iCent < 0) {
120  printf("I don't understand the meaning of centrality = %i\n",iCent);
121  return;
122  }
123  char buffer[1024];
124  sprintf(buffer,"%s/data/Sigmas.root",mDir);
125  if (mFile) {
126  delete mFile;
127  }
128  mFile = new TFile(buffer);
129 
130  char hist[1024];
131 
132  mnPts = getNPts();
133  if (ptBin > mnPts) {
134  printf("You asked for pt bin %i. There are only %i available.\n",ptBin, mnPts);
135  return;
136  }
137  int nCents, nx, ny;
138  if (0 == ptBin) {
139  nCents = getNCents();
140  nx = getNX(nCents);
141  ny = getNY(nCents);
142  } else {
143  nCents = getNPtCents();
144  nx = getNX(nCents);
145  ny = getNY(nCents);
146  }
147  if (iCent >= nCents) {
148  printf("You asked for centrality %i. There are only %i available (and we start counting at 0).\n",iCent, nCents);
149  return;
150  }
151 
152  char *type[] = {"NSig", "NDel", "NPlus", "NMinus", "NPlusMinus",
153  "PSig", "PDel", "PPlus", "PMinus", "PPlusMinus",
154  "PNSig", "PNDel", "PNPlus", "PNMinus", "PNPlusMinus",
155  "PNMinusPlus"};
156  char *index[] = {"", "", "", "", "",
157  "1", "1", "1", "1", "1",
158  "1", "1", "1", "1", "1",
159  "1"};
160  int order[] = {1, 4, 7, 10, 13, 2, 5, 8, 11, 14, 3, 6, 9, 12, 15};
161  TH2D *tmp;
162 
163  if (!mCanType) {
164  mCan = new TCanvas("c","c",3*250,5*250);
165  mCan->Divide(3,5);
166  } else if (strcmp(mCanType,"plotCentrality")) {
167  mCan->SetWindowSize(3*250,5*250);
168  mCan->Clear();
169  mCan->Divide(3,5);
170  }
171  mCanType = "plotCentrality";
172  // Omit plotting PNMinusPlus. Should be same as PNPlusMinus.
173  // (Hope I don't miss a big surprise.)
174  for (int iT=0;iT<15;iT++) {
175  mCan->cd(order[iT]);
176  if (0 == ptBin) {
177  sprintf( hist, "%s%s_%i", type[iT], index[iT], iCent );
178  } else {
179  sprintf( hist, "%s%s_%i_%i", type[iT], index[iT], iCent, ptBin-1 );
180  }
181  tmp = (TH2D *) gDirectory->Get(hist);
182  double val0 = tmp->GetBinContent(1,1);
183  double val1 = tmp->GetBinContent(24,16);
184  if (val1 < val0) {
185  gPad->SetPhi(-155);
186  } else {
187  gPad->SetPhi();
188  }
189  if (0 == ptBin) {
190  sprintf(buffer,"%s, centrality %i",type[iT],iCent);
191  } else {
192  sprintf(buffer,"%s, centrality %i, pt Bin %i",type[iT],iCent,ptBin-1);
193  }
194  tmp->SetTitle(buffer);
195  TAxis *x = tmp->GetXaxis();
196  TAxis *y = tmp->GetYaxis();
197  TAxis *z = tmp->GetZaxis();
198  x->SetTitleSize(0.07);
199  x->SetTitleColor(1);
200  x->SetTitleOffset(1);
201  x->SetNdivisions(505);
202  x->SetLabelSize(0.05);
203  x->SetTitle("#delta_{#phi}");
204  y->SetTitleSize(0.07);
205  y->SetTitleOffset(1);
206  y->SetTitleColor(1);
207  y->SetNdivisions(505);
208  y->SetLabelSize(0.05);
209  y->SetTitle("#delta_{#eta}");
210  z->SetTitleSize(0.07);
211  z->SetTitleOffset(0.6);
212  z->SetTitleColor(13);
213  z->SetNdivisions(505);
214  z->SetLabelSize(0.05);
215  tmp->Draw("surf1");
216  }
217  if (print) {
218  if (0 == ptBin) {
219  sprintf(buffer,"%s/txt/centrality%i.eps",mDir,iCent);
220  } else {
221  sprintf(buffer,"%s/txt/centrality%i_ptBin%i.eps",mDir,iCent,ptBin-1);
222  }
223  mCan->Print(buffer);
224  }
225 }
226 void postProcess::plotType( char *pType, int ptBin, int print ) {
227  char *type[] = {"NSig", "NDel", "NPlus", "NMinus", "NPlusMinus",
228  "PSig", "PDel", "PPlus", "PMinus", "PPlusMinus",
229  "PNSig", "PNDel", "PNPlus", "PNMinus", "PNPlusMinus",
230  "PNMinusPlus"};
231  char *index[] = {"", "", "", "", "",
232  "1", "1", "1", "1", "1",
233  "1", "1", "1", "1", "1",
234  "1"};
235  int iT;
236  int found = 0;
237  for (iT=0;iT<16;iT++) {
238  if (!strcmp(pType,type[iT])) {
239  found = 1;
240  break;
241  }
242  }
243  if (!found) {
244  printf("I didn't understand type of %s.\n",pType);
245  return;
246  }
247  char buffer[1024];
248  sprintf(buffer,"%s/data/Sigmas.root",mDir);
249  if (mFile) {
250  delete mFile;
251  }
252  mFile = new TFile(buffer);
253 
254  mnPts = getNPts();
255  if (ptBin > mnPts) {
256  printf("I only found %i ptBins. You asked for number %i\n",mnPts,ptBin);
257  return;
258  }
259  // Check for number of centrality bins by trying to read histograms until failure.
260  int nCents, nx, ny;
261  if (0 == ptBin) {
262  nCents = getNCents();
263  nx = getNX(nCents);
264  ny = getNY(nCents);
265  } else {
266  nCents = getNPtCents();
267  nx = getNX(nCents);
268  ny = getNY(nCents);
269  }
270 
271  TH2D *tmp;
272  if (!mCanType) {
273  mCan = new TCanvas("c","c",nx*250,ny*250);
274  mCan->Divide(mnx,mny);
275  } else if (strcmp(mCanType,"plotType") || nx != mnx || ny != mny) {
276  mCan->SetWindowSize(nx*250,ny*250);
277  mCan->Clear();
278  mCan->Divide(nx,ny);
279  }
280  mnx = nx;
281  mny = ny;
282  mCanType = "plotType";
283 
284  char hist[1024];
285  for (int iB=0;iB<nCents;iB++) {
286  mCan->cd(iB+1);
287  if (0 == ptBin) {
288  sprintf( hist, "%s%s_%i", pType, index[iT], iB );
289  } else {
290  sprintf( hist, "%s%s_%i_%i", pType, index[iT], iB, ptBin-1 );
291  }
292  tmp = (TH2D *) gDirectory->Get(hist);
293  double val0 = tmp->GetBinContent(1,1);
294  double val1 = tmp->GetBinContent(24,16);
295  if (val1 < val0) {
296  gPad->SetPhi(-155);
297  } else {
298  gPad->SetPhi();
299  }
300  if (0 == ptBin) {
301  sprintf(buffer,"%s, centrality %i",pType,iB);
302  } else {
303  sprintf(buffer,"%s, centrality %i, pt Bin %i",pType,iB,ptBin-1);
304  }
305  tmp->SetTitle(buffer);
306  TAxis *x = tmp->GetXaxis();
307  TAxis *y = tmp->GetYaxis();
308  TAxis *z = tmp->GetZaxis();
309  x->SetTitleSize(0.07);
310  x->SetTitleColor(1);
311  x->SetTitleOffset(1);
312  x->SetNdivisions(505);
313  x->SetLabelSize(0.05);
314  x->SetTitle("#delta_{#phi}");
315  y->SetTitleSize(0.07);
316  y->SetTitleOffset(1);
317  y->SetTitleColor(1);
318  y->SetNdivisions(505);
319  y->SetLabelSize(0.05);
320  y->SetTitle("#delta_{#eta}");
321  z->SetTitleSize(0.07);
322  z->SetTitleOffset(0.6);
323  z->SetTitleColor(13);
324  z->SetNdivisions(505);
325  z->SetLabelSize(0.05);
326  tmp->Draw("surf1");
327  }
328  if (print) {
329  if (0 == ptBin) {
330  sprintf(buffer,"%s/txt/%s.eps",mDir,pType);
331  } else {
332  sprintf(buffer,"%s/txt/%s_%i.eps",mDir,pType,ptBin-1);
333  }
334  mCan->Print(buffer);
335  }
336 }
337 void postProcess::centralityMults( int N0, int print ) {
338  char buffer[1024];
339  sprintf(buffer,"%s/QA/QA.root",mDir);
340  if (mFile) {
341  delete mFile;
342  }
343  mFile = new TFile(buffer);
344 
345  // Check for number of centrality bins by trying to read histograms until failure.
346  int mnCents = 0;
347  for (int iB=0;iB<100;iB++) {
348  sprintf(buffer,"multNSum_%i", iB );
349  if (gDirectory->Get(buffer)) {
350  mnCents++;
351  } else{
352  break;
353  }
354  }
355  int mnPtCents = 0;
356  for (int iB=0;iB<100;iB++) {
357  sprintf(buffer,"multNSum_%i_0", iB );
358  if (gDirectory->Get(buffer)) {
359  mnPtCents++;
360  } else{
361  break;
362  }
363  }
364  mnx = getNX(mnCents+1);
365  mny = getNY(mnCents+1);
366 
367  // Assume all histograms also have low and high pt version.
368  // All other routines in this file check not only for number of pt
369  // bins but how many centrality bins there are fro pt selected histograms.
370  // In this routine I am lazy (partly because I want to super-pose
371  // low and high pt on all-pt histogram versions.
372  float nBar[mnCents];
373  int num[mnCents];
374  TH1F *tmp = (TH1F *) gDirectory->Get("multNSum_0");
375  nBar[0] = tmp->GetMean();
376  num[0] = (int) tmp->GetEntries();
377  TH1F *mult = (TH1F *) tmp->Clone();
378  TH1F *mult34 = (TH1F *) tmp->Clone();
379  TH1F *multLow = (TH1F *) tmp->Clone();
380  TH1F *mult34Low = (TH1F *) tmp->Clone();
381  TH1F *multHigh = (TH1F *) tmp->Clone();
382  TH1F *mult34High = (TH1F *) tmp->Clone();
383 
384  for (int i=1;i<mnCents;i++) {
385  sprintf(buffer,"multNSum_%i",i);
386  TH1F *tmp = (TH1F *) gDirectory->Get(buffer);
387  nBar[i] = tmp->GetMean();
388  num[i] = (int) tmp->GetEntries();
389  mult->Add(tmp);
390  }
391  for (int i=1;i<mnPtCents;i++) {
392  sprintf(buffer,"multNSum_%i_0",i);
393  tmp = (TH1F *) gDirectory->Get(buffer);
394  multLow->Add(tmp);
395  sprintf(buffer,"multNSum_%i_1",i);
396  tmp = (TH1F *) gDirectory->Get(buffer);
397  multHigh->Add(tmp);
398  }
399 
400  for (int i=1;i<=mult->GetNbinsX();i++) {
401  float val = pow(i+0.0,0.75) * mult->GetBinContent(i);
402  mult34->SetBinContent(i,val);
403  val = pow(i+0.0,0.75) * multLow->GetBinContent(i);
404  mult34Low->SetBinContent(i,val);
405  val = pow(i+0.0,0.75) * multHigh->GetBinContent(i);
406  mult34High->SetBinContent(i,val);
407  }
408 
409  FILE *out;
410  char outFileName[1024];
411  sprintf( outFileName, "%s/txt/NBar.txt", mDir );
412  out = fopen(outFileName,"w");
413 
414  for (int i=0;i<mnCents;i++) {
415  fprintf(out,"Centrality %i has mean multiplicity of %f with %i events\n",i,nBar[i],num[i]);
416  }
417  fclose(out);
418 
419  if (!mCanType) {
420  mCan = new TCanvas("c","c",mnx*250,mny*250);
421  mCan->Divide(mnx,mny);
422  } else if (strcmp(mCanType,"centralityMults")) {
423  mCan->SetWindowSize(mnx*250,mny*250);
424  mCan->Clear();
425  mCan->Divide(mnx,mny);
426  }
427  mCanType = "centralityMults";
428 
429  for (int i=0;i<mnCents;i++) {
430  mCan->cd(i+1);
431  sprintf(buffer,"multNSum_%i",i);
432  tmp = (TH1F *) gDirectory->Get(buffer);
433  tmp->SetAxisRange(0,N0);
434  tmp->SetLineColor(1);
435  sprintf(buffer,"centrality bin %i",i);
436  tmp->SetTitle(buffer);
437  tmp->Draw();
438  if (mnPtCents == mnCents) {
439  sprintf(buffer,"multNSum_%i_0",i);
440  tmp = (TH1F *) gDirectory->Get(buffer);
441  tmp->SetAxisRange(0,N0);
442  tmp->SetLineColor(2);
443  tmp->Draw("same");
444 
445  sprintf(buffer,"multNSum_%i_1",i);
446  TH1F *tmp = (TH1F *) gDirectory->Get(buffer);
447  sprintf(buffer,"centrality bin %i",i);
448  tmp->SetTitle(buffer);
449  tmp->SetAxisRange(0,N0);
450  tmp->SetLineColor(3);
451  tmp->Draw("same");
452  }
453  }
454  mCan->cd(mnCents+1);
455  mult34Low->SetAxisRange(0,N0);
456  mult34Low->SetLineColor(3);
457  mult34Low->SetTitle("m^3/4*dN/dm");
458  mult34Low->Draw();
459  mult34High->SetAxisRange(0,N0);
460  mult34High->SetLineColor(2);
461  mult34High->Draw("same");
462 
463  mult34->SetTitle("m^3/4*dN/dm");
464  mult34->SetAxisRange(0,N0);
465  mult34->SetLineColor(1);
466  mult34->Draw("same");
467 
468  if (print) {
469  sprintf(buffer,"%s/txt/MultBins.eps",mDir);
470  mCan->Print(buffer);
471  }
472 }
473 void postProcess::multiplicity( int N0, int print ) {
474  char buffer[1024];
475  sprintf(buffer,"%s/data/Data.root",mDir);
476  if (mFile) {
477  delete mFile;
478  }
479  mFile = new TFile(buffer);
480 
481  TH1F *mult, *mult34;
482  mult = (TH1F *) gDirectory->Get("Multiplicity");
483  mult34 = (TH1F *) mult->Clone();
484  TAxis *x = mult->GetXaxis();
485  TAxis *y = mult->GetYaxis();
486  x->SetTitle("Mult");
487  y->SetTitle("#frac{dN}{dM}");
488  x = mult34->GetXaxis();
489  y = mult34->GetYaxis();
490  x->SetTitle("Mult");
491  y->SetTitle("#frac{dN}{dM} * M^{3/4}");
492  y->SetTitleOffset(1.75);
493 
494  FILE *out;
495  char outFileName[1024];
496  sprintf( outFileName, "%s/txt/Multiplicities.txt", mDir );
497  out = fopen(outFileName,"w");
498  for (int i=1;i<=mult->GetNbinsX();i++) {
499  float val = pow(i+0.0,0.75) * mult->GetBinContent(i);
500  mult34->SetBinContent(i,val);
501  fprintf(out,"%i %f\n",i,mult->GetBinContent(i));
502  }
503  fclose(out);
504 
505  if (!mCanType) {
506  mCan = new TCanvas("c","c",800,400);
507  mCan->Divide(2,1);
508  } else if (strcmp(mCanType,"multiplicity")) {
509  mCan->SetWindowSize(800,400);
510  mCan->Clear();
511  mCan->Divide(2,1);
512  }
513  mCanType = "multiplicity";
514 
515  mCan->cd(1);
516  gPad->SetLogy(1);
517  mult->SetAxisRange(0,N0);
518  mult->SetTitle("Multiplicity");
519  mult->Draw();
520  mCan->cd(2);
521  mult34->SetAxisRange(0,N0);
522  mult34->SetTitle("Multiplicity");
523  mult34->Draw();
524  if (print) {
525  sprintf(buffer,"%s/txt/Multiplicity.eps",mDir);
526  mCan->Print(buffer);
527  }
528 }
529 void postProcess::pt( int P0, int print ) {
530  char buffer[1024];
531  sprintf(buffer,"%s/data/Data.root",mDir);
532  if (mFile) {
533  delete mFile;
534  }
535  mFile = new TFile(buffer);
536 
537  TH1F *pt, *pt34;
538  pt = (TH1F *) gDirectory->Get("Pt");
539  pt34 = (TH1F *) pt->Clone();
540  TAxis *x = pt->GetXaxis();
541  TAxis *y = pt->GetYaxis();
542  x->SetTitle("P_{#perp}");
543  y->SetTitle("#frac{dN}{dP_{#perp}}");
544  x = pt34->GetXaxis();
545  y = pt34->GetYaxis();
546  x->SetTitle("P_{#perp}");
547  y->SetTitle("#frac{dN}{dP_{#perp}} * P_{#perp}^{3/4}");
548  y->SetTitleOffset(1.75);
549 
550  FILE *out;
551  char outFileName[1024];
552  sprintf( outFileName, "%s/txt/Pt.txt", mDir );
553  out = fopen(outFileName,"w");
554  for (int i=1;i<=pt->GetNbinsX();i++) {
555  float val = pow(i+0.0,0.75) * pt->GetBinContent(i);
556  pt34->SetBinContent(i,val);
557  fprintf(out,"%i %f\n",i,pt->GetBinContent(i));
558  }
559  fclose(out);
560 
561  if (!mCanType) {
562  mCan = new TCanvas("c","c",800,400);
563  mCan->Divide(2,1);
564  }
565  if (!mCanType || strcmp(mCanType,"pt")) {
566  mCan->SetWindowSize(800,400);
567  mCan->Clear();
568  mCan->Divide(2,1);
569  }
570  mCanType = "pt";
571 
572  mCan->cd(1);
573  gPad->SetLogy(1);
574  pt->SetAxisRange(0,P0);
575  pt->SetTitle("Pt");
576  pt->Draw();
577  mCan->cd(2);
578  pt34->SetAxisRange(0,P0);
579  pt34->SetTitle("Pt");
580  pt34->Draw();
581  if (print) {
582  sprintf(buffer,"%s/txt/Pt.eps",mDir);
583  mCan->Print(buffer);
584  }
585 }
586 void postProcess::deltaSigma2( int print ) {
587  char buffer[1024];
588  sprintf(buffer,"%s/data/Sigmas.root",mDir);
589  if (mFile) {
590  delete mFile;
591  }
592  mFile = new TFile(buffer);
593 
594  // Check for number of centrality bins by trying to read histograms until failure.
595  char hist[1024];
596  mnCents = getNCents();
597  mnx = getNX(mnCents);
598  mny = getNY(mnCents);
599 
600  FILE *out;
601  char outFileName[1024];
602  sprintf( outFileName, "%s/txt/NBar.txt", mDir );
603  out = fopen(outFileName,"a+");
604 
605  // Print out ptHat values for full scale.
606  TH2D *pthat;
607  fprintf(out,"\n");
608  fprintf(out,"\n");
609  fprintf(out,"\n");
610  for (int iB=0;iB<mnCents;iB++) {
611  sprintf(hist,"SptHat_%i", iB );
612  pthat = (TH2D *) gDirectory->Get(hist);
613  fprintf(out," Centrality %i has ptHat of %f\n",iB,pthat->GetBinContent(24,16));
614  }
615  fprintf(out,"\n");
616  fprintf(out,"\n");
617  fprintf(out,"\n");
618  for (int iB=0;iB<mnCents;iB++) {
619  sprintf(hist,"PptHat_%i", iB );
620  pthat = (TH2D *) gDirectory->Get(hist);
621  fprintf(out," Centrality %i has ptHat+ of %f\n",iB,pthat->GetBinContent(24,16));
622  }
623  fprintf(out,"\n");
624  fprintf(out,"\n");
625  fprintf(out,"\n");
626  for (int iB=0;iB<mnCents;iB++) {
627  sprintf(hist,"MptHat_%i", iB );
628  pthat = (TH2D *) gDirectory->Get(hist);
629  fprintf(out," Centrality %i has ptHat- of %f\n",iB,pthat->GetBinContent(24,16));
630  }
631  // Print out \sigma_{hat_p_perp} values for full scale.
632  fprintf(out,"\n");
633  fprintf(out,"\n");
634  fprintf(out,"\n");
635  for (int iB=0;iB<mnCents;iB++) {
636  sprintf(hist,"SsigPtHat_%i", iB );
637  pthat = (TH2D *) gDirectory->Get(hist);
638  fprintf(out," Centrality %i has sigma_{hat_p_perp} of %f\n",iB,sqrt(pthat->GetBinContent(24,16)));
639  }
640  fclose(out);
641 
642  // Plot \Delta\sigma^2 histograms.
643  TH2D *psig;
644  if (!mCanType || strcmp(mCanType,"deltaSigma2")) {
645  mCan = new TCanvas("c","c",mnx*250,mny*250);
646  mCan->SetWindowPosition(1,1);
647  mCan->Divide(mnx,mny);
648  mCanType = "deltaSigma2";
649  }
650  char title[1024];
651  for (int iB=0;iB<mnCents;iB++) {
652  sprintf(hist,"PSig1_%i", iB );
653  psig = (TH2D *) gDirectory->Get(hist);
654 
655  TAxis *x = psig->GetXaxis();
656  TAxis *y = psig->GetYaxis();
657  TAxis *z = psig->GetZaxis();
658  x->SetTitleSize(0.075);
659  x->SetTitleColor(1);
660  x->SetTitleOffset(1);
661  x->SetNdivisions(505);
662  x->SetLabelSize(0.05);
663  x->SetTitle("#delta_{#phi}");
664  y->SetTitleSize(0.075);
665  y->SetTitleOffset(1);
666  y->SetTitleColor(1);
667  y->SetNdivisions(505);
668  y->SetLabelSize(0.05);
669  y->SetTitle("#delta_{#eta}");
670  z->SetTitleSize(0.08);
671  z->SetTitleOffset(0.6);
672  z->SetTitleColor(13);
673  z->SetNdivisions(505);
674  z->SetLabelSize(0.05);
675  sprintf(title,"#Delta#sigma^{2}, bin %i",iB);
676  psig->SetTitle(title);
677 
678  mCan->cd(iB+1);
679  psig->Draw("surf1");
680  }
681  if (print) {
682  sprintf(buffer,"%s/txt/Deltasigma2.eps",mDir);
683  mCan->Print(buffer);
684  }
685 }
686 void postProcess::exportForMatlab() {
687  char buffer[1024];
688  sprintf( buffer, "%s/data/Sigmas.root", mDir );
689  if (mFile) {
690  delete mFile;
691  }
692  mFile = new TFile(buffer);
693 
694  char *measName[] = {"NSig", "PSig", "PNSig",
695  "NDel", "PDel", "PNDel",
696  "NPlus", "PPlus", "PNPlus",
697  "NMinus", "PMinus", "PNMinus",
698  "NPlusMinus", "PPlusMinus", "PNPlusMinus", "PNMinusPlus"};
699  char *outName[] = {"N_CI", "P_CI", "PN_CI",
700  "N_CD", "P_CD", "PN_CD",
701  "N_Plus", "P_Plus", "PN_Plus",
702  "N_Minus", "P_Minus", "PN_Minus",
703  "N_PlusMinus", "P_PlusMinus", "PN_PlusMinus", "PN_MinusPlus"};
704  char *index[] = {"", "1", "1",
705  "", "1", "1",
706  "", "1", "1",
707  "", "1", "1",
708  "", "1", "1", "1"};
709  char *centName[11];
710  char *ptRange[] = {"0.15-0.5", "0.5-2.0"};
711 
712  FILE *out;
713  char outFileName[1024];
714  TH2D *hist;
715 
716  int nCents, nPtCents, nPts;
717  nCents = getNCents();
718  nPtCents = getNPtCents();
719  nPts = getNPts();
720 
721  if (1 == nCents) {
722  centName[0] = "b0-3";
723  } else if (2 == nCents) {
724  centName[0] = "50-100";
725  centName[1] = "0-50";
726  } else if (4 == nCents) {
727  centName[0] = "0-25";
728  centName[1] = "25-50";
729  centName[2] = "50-75";
730  centName[3] = "75-100";
731  } else if (5 == nCents) {
732  centName[0] = "1-3";
733  centName[1] = "4-6";
734  centName[2] = "7-9";
735  centName[3] = "10-12";
736  centName[4] = "13-24";
737  } else if (6 == nCents) {
738  centName[0] = "1-17";
739  centName[1] = "17-33";
740  centName[2] = "33-50";
741  centName[3] = "50-67";
742  centName[4] = "67-83";
743  centName[5] = "83-100";
744  } else if (7 == nCents) {
745  centName[0] = "80-100";
746  centName[1] = "64-80";
747  centName[2] = "48-64";
748  centName[3] = "32-48";
749  centName[4] = "16-32";
750  centName[5] = "8-16";
751  centName[6] = "0-8";
752  } else if (11 == nCents) {
753  centName[0] = "90-100";
754  centName[1] = "80-90";
755  centName[2] = "70-80";
756  centName[3] = "60-70";
757  centName[4] = "50-60";
758  centName[5] = "40-50";
759  centName[6] = "30-40";
760  centName[7] = "20-30";
761  centName[8] = "10-20";
762  centName[9] = "5-10";
763  centName[10] = "0-5";
764  }
765  for (int iD=0;iD<16;iD++) {
766  for (int iC=0;iC<nCents;iC++) {
767  sprintf( outFileName, "%s/txt/fluct_%s_%s.txt", mDir, outName[iD], centName[iC] );
768  sprintf(buffer,"%s%s_%i", measName[iD], index[iD], iC );
769  out = fopen(outFileName,"w");
770  hist = (TH2D *) gDirectory->Get(buffer);
771  for (int n=1;n<=hist->GetNbinsX();n++) {
772  for (int m=1;m<=hist->GetNbinsY();m++) {
773  fprintf(out,"%f\n",hist->GetBinContent(n,m));
774  }
775  }
776  fclose(out);
777 
778  sprintf( outFileName, "%s/txt/fluct_%s_%s_err.txt", mDir, outName[iD], centName[iC] );
779  sprintf(buffer,"%sErrors_%i", measName[iD], iC );
780  out = fopen(outFileName,"w");
781  hist = (TH2D *) gDirectory->Get(buffer);
782  for (int n=1;n<=hist->GetNbinsX();n++) {
783  for (int m=1;m<=hist->GetNbinsY();m++) {
784  fprintf(out,"%f\n",hist->GetBinContent(n,m));
785  }
786  }
787  fclose(out);
788  }
789  }
790  if (1 == nPtCents) {
791  centName[0] = "b0-3";
792  } else if (2 == nPtCents) {
793  centName[0] = "50-100";
794  centName[1] = "0-50";
795  } else if (4 == nPtCents) {
796  centName[0] = "0-25";
797  centName[1] = "25-50";
798  centName[2] = "50-75";
799  centName[3] = "75-100";
800  } else if (5 == nPtCents) {
801  centName[0] = "1-3";
802  centName[1] = "4-6";
803  centName[2] = "7-9";
804  centName[3] = "10-12";
805  centName[4] = "13-24";
806  } else if (6 == nCents) {
807  centName[0] = "1-17";
808  centName[1] = "17-33";
809  centName[2] = "33-50";
810  centName[3] = "50-67";
811  centName[4] = "67-83";
812  centName[5] = "83-100";
813  } else if (7 == nPtCents) {
814  centName[0] = "80-100";
815  centName[1] = "64-80";
816  centName[2] = "48-64";
817  centName[3] = "32-48";
818  centName[4] = "16-32";
819  centName[5] = "8-16";
820  centName[6] = "0-8";
821  } else if (11 == nPtCents) {
822  centName[0] = "90-100";
823  centName[1] = "80-90";
824  centName[2] = "70-80";
825  centName[3] = "60-70";
826  centName[4] = "50-60";
827  centName[5] = "40-50";
828  centName[6] = "30-40";
829  centName[7] = "20-30";
830  centName[8] = "10-20";
831  centName[9] = "5-10";
832  centName[10] = "0-5";
833  }
834  for (int iD=0;iD<16;iD++) {
835  for (int iC=0;iC<nPtCents;iC++) {
836  for (int iP=0;iP<nPts;iP++) {
837  sprintf( outFileName, "%s/txt/fluct_%s_%s_%s.txt", mDir, outName[iD], centName[iC], ptRange[iP] );
838  sprintf(buffer,"%s%s_%i_%i", measName[iD], index[iD], iC, iP );
839  out = fopen(outFileName,"w");
840  hist = (TH2D *) gDirectory->Get(buffer);
841  for (int n=1;n<=hist->GetNbinsX();n++) {
842  for (int m=1;m<=hist->GetNbinsY();m++) {
843  fprintf(out,"%f\n",hist->GetBinContent(n,m));
844  }
845  }
846  fclose(out);
847 
848  sprintf( outFileName, "%s/txt/fluct_%s_%s_%s_err.txt", mDir, outName[iD], centName[iC], ptRange[iP] );
849  sprintf(buffer,"%sErrors_%i_%i", measName[iD], iC, iP );
850  out = fopen(outFileName,"w");
851  hist = (TH2D *) gDirectory->Get(buffer);
852  for (int n=1;n<=hist->GetNbinsX();n++) {
853  for (int m=1;m<=hist->GetNbinsY();m++) {
854  fprintf(out,"%f\n",hist->GetBinContent(n,m));
855  }
856  }
857  fclose(out);
858  }
859  }
860  }
861 }
862 void postProcess::importFromMatlab() {
863 }
864 int postProcess::getNCents () {
865  char hist[1024];
866  int nCents = 0;
867  for (int iB=0;iB<100;iB++) {
868  sprintf(hist,"NSig_%i", iB );
869  if (gDirectory->Get(hist)) {
870  nCents++;
871  } else{
872  break;
873  }
874  }
875  return nCents;
876 }
877 int postProcess::getNPtCents () {
878  char hist[1024];
879  int nPtCents = 0;
880  for (int iC=0;iC<100;iC++) {
881  sprintf(hist,"NSig_%i_0", iC );
882  if (gDirectory->Get(hist)) {
883  nPtCents++;
884  } else{
885  break;
886  }
887  }
888  return nPtCents;
889 }
890 int postProcess::getNPts () {
891  char hist[1024];
892  int nPts = 0;
893  for (int iP=0;iP<100;iP++) {
894  sprintf(hist,"NSig_0_%i", iP );
895  if (gDirectory->Get(hist)) {
896  nPts++;
897  } else{
898  break;
899  }
900  }
901  return nPts;
902 }
903 int postProcess::getNX( int nCent ) {
904  int nx = 0;
905  if (nCent < 3) {
906  nx = 1;
907  } else if (nCent < 7) {
908  nx = 2;
909  } else if (nCent < 13) {
910  nx = 3;
911  } else {
912  nx = 4;
913  }
914  return nx;
915 }
916 int postProcess::getNY( int nCent ) {
917  int ny = 0;
918  if (nCent < 2) {
919  ny = 1;
920  } else if (nCent < 5) {
921  ny = 2;
922  } else if (nCent < 10) {
923  ny = 3;
924  } else {
925  ny = 4;
926  }
927  return ny;
928 }
929 
930 
931 void postProcess::plotMCInfo( int print ) {
932  char buffer[1024];
933  sprintf(buffer,"%s/QA/QA.root",mDir);
934  if (mFile) {
935  delete mFile;
936  }
937  mFile = new TFile(buffer);
938 
939  if (!mCanType || strcmp(mCanType,"plotMCInfo")) {
940  mCan = new TCanvas("c","c",2*400,2*400);
941  mCan->SetWindowPosition(1,1);
942  mCan->Divide(2,2);
943  mCanType = "plotMCInfo";
944  }
945 
946 
947  TH1D *hist;
948  TAxis *x, *y;
949 
950  hist = (TH1D *) gDirectory->Get("impact");
951  x = hist->GetXaxis();
952  y = hist->GetYaxis();
953  x->SetTitleSize(0.075);
954  x->SetTitleColor(1);
955  x->SetTitleOffset(1);
956  x->SetNdivisions(505);
957  x->SetLabelSize(0.05);
958  x->SetTitle("fm");
959  y->SetTitleSize(0.075);
960  y->SetTitleOffset(1);
961  y->SetTitleColor(1);
962  y->SetNdivisions(505);
963  y->SetLabelSize(0.05);
964  y->SetTitle("Events");
965  hist->SetTitle("Impact Parameter");
966  mCan->cd(1);
967  hist->Draw();
968 
969  hist = (TH1D *) gDirectory->Get("binary");
970  x = hist->GetXaxis();
971  y = hist->GetYaxis();
972  x->SetTitleSize(0.075);
973  x->SetTitleColor(1);
974  x->SetTitleOffset(1);
975  x->SetNdivisions(505);
976  x->SetLabelSize(0.05);
977  x->SetTitle("Number");
978  y->SetTitleSize(0.075);
979  y->SetTitleOffset(1);
980  y->SetTitleColor(1);
981  y->SetNdivisions(505);
982  y->SetLabelSize(0.05);
983  y->SetTitle("Events");
984  hist->SetTitle("Binary collisions");
985  mCan->cd(2);
986  gPad->SetLogy();
987  hist->Draw();
988 
989  hist = (TH1D *) gDirectory->Get("participant");
990  x = hist->GetXaxis();
991  y = hist->GetYaxis();
992  x->SetTitleSize(0.075);
993  x->SetTitleColor(1);
994  x->SetTitleOffset(1);
995  x->SetNdivisions(505);
996  x->SetLabelSize(0.05);
997  x->SetTitle("Number");
998  y->SetTitleSize(0.075);
999  y->SetTitleOffset(1);
1000  y->SetTitleColor(1);
1001  y->SetNdivisions(505);
1002  y->SetLabelSize(0.05);
1003  y->SetTitle("Events");
1004  hist->SetTitle("Participants");
1005  mCan->cd(3);
1006  gPad->SetLogy();
1007  hist->Draw();
1008 
1009  if (print) {
1010  sprintf(buffer,"%s/txt/MCInfo.eps",mDir);
1011  mCan->Print(buffer);
1012  }
1013 }
1014 void postProcess::plotBinaryCollisions( int print ) {
1015  char buffer[1024];
1016  sprintf(buffer,"%s/QA/QA.root",mDir);
1017  if (mFile) {
1018  delete mFile;
1019  }
1020  mFile = new TFile(buffer);
1021 
1022  // Check for number of centrality bins by trying to read histograms until failure.
1023  mnCents = getNCents();
1024  mnx = getNX(mnCents+1);
1025  mny = getNY(mnCents+1);
1026 
1027  if (!mCanType ||
1028  (strcmp(mCanType,"plotBinaryCollisions") &&
1029  strcmp(mCanType,"plotParticipants"))) {
1030  mCan = new TCanvas("c","c",mnx*250,mny*250);
1031  mCan->SetWindowPosition(1,1);
1032  mCan->Divide(2,3);
1033  mCanType = "plotBinaryCollisions";
1034  }
1035 
1036 
1037  TH1D *bin;
1038  for (int iB=0;iB<mnCents;iB++) {
1039  sprintf(buffer,"binary_%i", iB );
1040  bin = (TH1D *) gDirectory->Get(buffer);
1041 
1042  TAxis *x = bin->GetXaxis();
1043  TAxis *y = bin->GetYaxis();
1044  x->SetTitleSize(0.075);
1045  x->SetTitleColor(1);
1046  x->SetTitleOffset(1);
1047  x->SetNdivisions(505);
1048  x->SetLabelSize(0.05);
1049  x->SetTitle("binary collisions");
1050  y->SetTitleSize(0.075);
1051  y->SetTitleOffset(1);
1052  y->SetTitleColor(1);
1053  y->SetNdivisions(505);
1054  y->SetLabelSize(0.05);
1055  y->SetTitle("Events");
1056  sprintf(buffer,"impact class %i",iB);
1057  bin->SetTitle(buffer);
1058 
1059  gPad->SetLogy();
1060  mCan->cd(iB+1);
1061  bin->Draw();
1062  }
1063  if (print) {
1064  sprintf(buffer,"%s/txt/binaryCollisions.eps",mDir);
1065  mCan->Print(buffer);
1066  }
1067 }
1068 void postProcess::plotParticipants( int print ) {
1069  char buffer[1024];
1070  sprintf(buffer,"%s/QA/QA.root",mDir);
1071  if (mFile) {
1072  delete mFile;
1073  }
1074  mFile = new TFile(buffer);
1075 
1076  // Check for number of centrality bins by trying to read histograms until failure.
1077  mnCents = getNCents();
1078  mnx = getNX(mnCents+1);
1079  mny = getNY(mnCents+1);
1080 
1081  if (!mCanType ||
1082  (strcmp(mCanType,"plotBinaryCollisions") &&
1083  strcmp(mCanType,"plotParticipants"))) {
1084  mCan = new TCanvas("c","c",mnx*250,mny*250);
1085  mCan->SetWindowPosition(1,1);
1086  mCan->Divide(2,3);
1087  mCanType = "plotBinaryCollisions";
1088  }
1089 
1090  TH1D *par;
1091  for (int iB=0;iB<mnCents;iB++) {
1092  sprintf(buffer,"participant_%i", iB );
1093  par = (TH1D *) gDirectory->Get(buffer);
1094 
1095  TAxis *x = par->GetXaxis();
1096  TAxis *y = par->GetYaxis();
1097  x->SetTitleSize(0.075);
1098  x->SetTitleColor(1);
1099  x->SetTitleOffset(1);
1100  x->SetNdivisions(505);
1101  x->SetLabelSize(0.05);
1102  x->SetTitle("participants");
1103  y->SetTitleSize(0.075);
1104  y->SetTitleOffset(1);
1105  y->SetTitleColor(1);
1106  y->SetNdivisions(505);
1107  y->SetLabelSize(0.05);
1108  y->SetTitle("Events");
1109  sprintf(buffer,"impact class %i",iB);
1110  par->SetTitle(buffer);
1111 
1112  mCan->cd(iB+1);
1113  gPad->SetLogy();
1114  par->Draw();
1115  }
1116  if (print) {
1117  sprintf(buffer,"%s/txt/participants.eps",mDir);
1118  mCan->Print(buffer);
1119  }
1120 }
1121 void postProcess::exportMCInfo() {
1122  char buffer[1024];
1123  sprintf(buffer,"%s/QA/QA.root",mDir);
1124  if (mFile) {
1125  delete mFile;
1126  }
1127  mFile = new TFile(buffer);
1128 
1129  // Check for number of centrality bins by trying to read histograms until failure.
1130  mnCents = getNCents();
1131 
1132  // Write binary collisions and participants.
1133  // Each centrality to its own file.
1134  FILE *out;
1135  char outFileName[1024];
1136  TH1D *hist;
1137  for (int iC=0;iC<mnCents;iC++) {
1138  sprintf( outFileName, "%s/txt/binaryCollisions_%i.txt", mDir, iC );
1139  out = fopen(outFileName,"w");
1140  sprintf(buffer,"binary_%i", iC );
1141  hist = (TH1D *) gDirectory->Get(buffer);
1142  for (int n=1;n<=hist->GetNbinsX();n++) {
1143  fprintf(out,"%i %6.f\n", n, hist->GetBinContent(n));
1144  }
1145  fclose(out);
1146 
1147  sprintf( outFileName, "%s/txt/participants_%i.txt", mDir, iC );
1148  out = fopen(outFileName,"w");
1149  sprintf(buffer,"participant_%i", iC );
1150  hist = (TH1D *) gDirectory->Get(buffer);
1151  for (int n=1;n<=hist->GetNbinsX();n++) {
1152  fprintf(out,"%i %6.f\n", n, hist->GetBinContent(n));
1153  }
1154  fclose(out);
1155  }
1156 
1157  sprintf( outFileName, "%s/txt/impactParameter.txt", mDir );
1158  out = fopen(outFileName,"w");
1159  hist = (TH1D *) gDirectory->Get("impact");
1160  for (int n=1;n<=hist->GetNbinsX();n++) {
1161  fprintf(out,"%8.4f %6.f\n",hist->GetBinCenter(n), hist->GetBinContent(n));
1162  }
1163  fclose(out);
1164 }
Generated by   doxygen 1.8.5