StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
directCumulants_v2.C
1 //
3 // $Id: directCumulants_v2.C,v 1.3 2010/06/10 16:33:57 posk Exp $
4 //
5 // Authors: Dhevan Gangadharan, UCLA, Jan 2010
6 //
8 //
9 // Description: Macro to sum the Terms from the output of StFlowDirectCumulantMaker,
10 // plot v_2{2] and v_2{4}, and output the graphs to flow.diCrumulant.graphs.root.
11 // You can plot them again with plotDirectCumulants_v2pt.C .
12 // cent=1 is the most central, cent=0 is the yield weighted average of all the others.
13 //
15 
16 #include <math.h>
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <Riostream.h>
20 
21 #include "TFile.h"
22 #include "TString.h"
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TGraphErrors.h"
26 #include "TArray.h"
27 #include "TLegend.h"
28 #include "TStyle.h"
29 
30 void directCumulants_v2(){
31 
32  TCanvas* cOld = (TCanvas*)gROOT->GetListOfCanvases(); // delete old canvas
33  if (cOld) cOld->Delete();
34 
35  bool ptGraphs = kFALSE;
36 
37  //TFile *myfile = new TFile("Y4_test4.root","READ");
38  TFile *myfile = new TFile("flow.dirCumulant.root","READ");
39  //TFile *myfile = new TFile("~dhevang/test.root","READ");
40  const int CENTBINS= 9;
41  float x[CENTBINS] = {2.5, 7.5, 15., 25., 35., 45., 55., 65., 75.};
42 
43  gROOT->SetStyle("Plain"); // set style
44  //gStyle->SetOptDate(0);
45  gStyle->SetOptStat(0);
46 
47  const int TERMS= 10; // terms needed to calculate the cumulants
48  const int TYPES= 2; // differential or integrated
49  const int PHASES= 2; // sin or cos
50  const int SPECIES= 1; // charged particles = 0
51  const int PTBINS= 62; // max pt bins
52  const int PTINTBINS=19; // max pt bins for pt integration
53  int ptOverflow = PTBINS-1;
54  int ptCu = PTBINS-2; // integrated cumulants
55 
56  double yields[SPECIES][PTBINS]={{0}};
57  double cent_yields[SPECIES][CENTBINS][PTBINS]={{{0}}};
58 
59  // centrality 0 is min bias
60  double v2_2[SPECIES][CENTBINS+1][PTBINS]={{{0.}}}; // differential
61  double v2_2_e[SPECIES][CENTBINS+1][PTBINS]={{{0.}}};
62  double v2_4[SPECIES][CENTBINS+1][PTBINS]={{{0.}}};
63  double v2_4_e[SPECIES][CENTBINS+1][PTBINS]={{{0.}}};
64  double v2IntPt_2[SPECIES][CENTBINS+1]={{0.}}; // integral
65  double v2IntPt_2_e[SPECIES][CENTBINS+1]={{0.}};
66  double v2IntPt_4[SPECIES][CENTBINS+1]={{0.}};
67  double v2IntPt_4_e[SPECIES][CENTBINS+1]={{0.}};
68  double used_cent_yields2[SPECIES][CENTBINS+1]={{0.}};
69  double used_cent_yields4[SPECIES][CENTBINS+1]={{0.}};
70  double usedInt_cent_yields2[SPECIES][CENTBINS+1]={{0.}};
71  double usedInt_cent_yields4[SPECIES][CENTBINS+1]={{0.}};
72 
73  TString *name;
74  TString *name1;
75  TString *name2;
76  TString *name3;
77  TString *name4;
78 
79  double Terms[TERMS][TYPES][PHASES][SPECIES][CENTBINS][PTBINS]={{{{{{0}}}}}};
80  double Terms_e[TERMS][TYPES][PHASES][SPECIES][CENTBINS][PTBINS]={{{{{{0}}}}}};
81 
83  // This part takes all the info out of the input file and puts it into the "Terms" arrays
84 
85  for(int species=0; species<SPECIES; species++){
86 
87  name = new TString("Cent_Pt_yields_");
88  *name += species;
89 
90  for(int cent=0; cent<CENTBINS; cent++){
91 
92  for(int ptbin=0; ptbin<PTBINS; ptbin++){
93 
94  cent_yields[species][cent][ptbin] = ((TH2D*)myfile->Get(name->Data()))->GetBinContent(cent+1,ptbin+1);
95  yields[species][ptbin] += cent_yields[species][cent][ptbin];
96 
97  }
98  }
99  }
100 
101  for(int term=0; term<TERMS; term++){
102 
103  name1 = new TString("Term_");
104  *name1 += term+1;
105 
106  for(int type=0; type<TYPES; type++){
107 
108  name2 = new TString();
109  if(type==0) name2->Append("_D_"); // differential
110  else name2->Append("_I_"); // integral
111 
112  for(int phase=0; phase<PHASES; phase++){
113 
114  name3 = new TString();
115  if(phase==0) name3->Append("cos_");
116  else name3->Append("sin_");
117  for(int species=0; species<SPECIES; species++){
118 
119  name4 = new TString();
120  *name4 += species;
121 
122  for(int cent=0; cent<CENTBINS; cent++){
123 
124  for(int ptbin=0; ptbin<PTBINS; ptbin++){
125 
126  if(cent_yields[species][cent][ptbin] < 1) continue;
127 
128  name = new TString();
129  name->Append(name1->Data());
130  name->Append(name2->Data());
131  name->Append(name3->Data());
132  name->Append(name4->Data());
133 
134  Terms[term][type][phase][species][cent][ptbin] = ((TH2D*)myfile->Get(name->Data()))->GetBinContent(cent+1,ptbin+1)/cent_yields[species][cent][ptbin];
135 
136  name->Append("_Sq");
137  Terms_e[term][type][phase][species][cent][ptbin] = ((TH2D*)myfile->Get(name->Data()))->GetBinContent(cent+1,ptbin+1)/cent_yields[species][cent][ptbin];
138  Terms_e[term][type][phase][species][cent][ptbin] -= pow(Terms[term][type][phase][species][cent][ptbin],2);
139  Terms_e[term][type][phase][species][cent][ptbin] /= cent_yields[species][cent][ptbin];
140  Terms_e[term][type][phase][species][cent][ptbin] = sqrt(Terms_e[term][type][phase][species][cent][ptbin]);
141 
142  }
143  }
144  }
145  }
146  }
147  }
148 
149 
151  int COS=0, SIN=1, I=1;
152 
153  double Cumulant_term[11][TYPES][SPECIES][CENTBINS][PTBINS]={{{{{0}}}}};
154  double Cumulant_term_e[11][TYPES][SPECIES][CENTBINS][PTBINS]={{{{{0}}}}};
155  double Cumulant_4[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
156  double Cumulant_4_e[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
157  double Cumulant_2[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
158  double Cumulant_2_e[TYPES][SPECIES][CENTBINS][PTBINS]={{{{0}}}};
159 
160  // This part reconstructs all the Cumulant terms from "Terms"
161  for(int type=0; type<TYPES; type++){
162 
163  for(int species=0; species<SPECIES; species++){
164 
165  for(int cent=0; cent<CENTBINS; cent++){
166 
167  for(int ptbin=0; ptbin<PTBINS; ptbin++){
168 
169  Cumulant_term[0][type][species][cent][ptbin] = Terms[0][type][COS][species][cent][ptbin];
170  Cumulant_term_e[0][type][species][cent][ptbin] = Terms_e[0][type][COS][species][cent][ptbin];
171 
172  Cumulant_term[1][type][species][cent][ptbin] = Terms[1][type][COS][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin];
173  Cumulant_term[1][type][species][cent][ptbin] -= Terms[1][type][SIN][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin];
174  Cumulant_term_e[1][type][species][cent][ptbin] = pow(Terms_e[1][type][COS][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin],2);
175  Cumulant_term_e[1][type][species][cent][ptbin] += pow(Terms[1][type][COS][species][cent][ptbin]*Terms_e[2][I][COS][species][cent][ptbin],2);
176  Cumulant_term_e[1][type][species][cent][ptbin] += pow(Terms_e[1][type][SIN][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin],2);
177  Cumulant_term_e[1][type][species][cent][ptbin] += pow(Terms[1][type][SIN][species][cent][ptbin]*Terms_e[2][I][SIN][species][cent][ptbin],2);
178  Cumulant_term_e[1][type][species][cent][ptbin] = sqrt(Cumulant_term_e[1][type][species][cent][ptbin]);
179 
180  Cumulant_term[2][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[4][I][COS][species][cent][ptbin];
181  Cumulant_term[2][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[4][I][SIN][species][cent][ptbin];
182 
183  Cumulant_term[3][type][species][cent][ptbin] = Terms[5][I][COS][species][cent][ptbin]*Terms[6][type][COS][species][cent][ptbin];
184  Cumulant_term[3][type][species][cent][ptbin] -= Terms[5][I][SIN][species][cent][ptbin]*Terms[6][type][SIN][species][cent][ptbin];
185 
186  Cumulant_term[4][type][species][cent][ptbin] = Terms[7][I][COS][species][cent][ptbin]*Terms[8][type][COS][species][cent][ptbin];
187  Cumulant_term[4][type][species][cent][ptbin] -= Terms[7][I][SIN][species][cent][ptbin]*Terms[8][type][SIN][species][cent][ptbin];
188 
189  Cumulant_term[5][type][species][cent][ptbin] = Terms[9][type][COS][species][cent][ptbin]*Terms[9][I][COS][species][cent][ptbin];
190  Cumulant_term[5][type][species][cent][ptbin] -= Terms[9][type][SIN][species][cent][ptbin]*Terms[9][I][SIN][species][cent][ptbin];
191 
192  Cumulant_term[6][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin];
193  Cumulant_term[6][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[2][I][COS][species][cent][ptbin];
194  Cumulant_term[6][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin];
195  Cumulant_term[6][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[2][I][SIN][species][cent][ptbin];
196 
197  Cumulant_term_e[6][type][species][cent][ptbin] = pow(Terms[3][type][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms_e[2][I][COS][species][cent][ptbin],2);
198  Cumulant_term_e[6][type][species][cent][ptbin] = sqrt(Cumulant_term_e[6][I][species][cent][ptbin]);
199 
200 
201  Cumulant_term[7][type][species][cent][ptbin] = Terms[3][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[1][type][COS][species][cent][ptbin];
202  Cumulant_term[7][type][species][cent][ptbin] -= Terms[3][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[1][type][COS][species][cent][ptbin];
203  Cumulant_term[7][type][species][cent][ptbin] -= Terms[3][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[1][type][SIN][species][cent][ptbin];
204  Cumulant_term[7][type][species][cent][ptbin] -= Terms[3][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[1][type][SIN][species][cent][ptbin];
205 
206  Cumulant_term_e[7][type][species][cent][ptbin] = pow(Terms[3][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms_e[1][type][COS][species][cent][ptbin],2);
207  Cumulant_term_e[7][type][species][cent][ptbin] = sqrt(Cumulant_term_e[7][type][species][cent][ptbin]);
208 
209 
210  Cumulant_term[8][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[9][I][COS][species][cent][ptbin];
211  Cumulant_term[8][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[9][I][COS][species][cent][ptbin];
212  Cumulant_term[8][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[9][I][SIN][species][cent][ptbin];
213  Cumulant_term[8][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[9][I][SIN][species][cent][ptbin];
214 
215  Cumulant_term[9][type][species][cent][ptbin] = Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[9][type][COS][species][cent][ptbin];
216  Cumulant_term[9][type][species][cent][ptbin] -= Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[9][type][COS][species][cent][ptbin];
217  Cumulant_term[9][type][species][cent][ptbin] -= Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[9][type][SIN][species][cent][ptbin];
218  Cumulant_term[9][type][species][cent][ptbin] -= Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[9][type][SIN][species][cent][ptbin];
219 
220  Cumulant_term[10][type][species][cent][ptbin] = Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
221  Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
222  Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
223  Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
224  Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin];
225  Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
226  Cumulant_term[10][type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[5][I][COS][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
227  Cumulant_term[10][type][species][cent][ptbin] += Terms[3][type][SIN][species][cent][ptbin]*Terms[5][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin]*Terms[7][I][SIN][species][cent][ptbin];
228 
229 
230  Cumulant_4[type][species][cent][ptbin] = Cumulant_term[0][type][species][cent][ptbin] -2*Cumulant_term[1][type][species][cent][ptbin];
231  Cumulant_4[type][species][cent][ptbin] -= Cumulant_term[2][type][species][cent][ptbin] + Cumulant_term[3][type][species][cent][ptbin];
232  Cumulant_4[type][species][cent][ptbin] -= 2*Cumulant_term[4][type][species][cent][ptbin];
233  Cumulant_4[type][species][cent][ptbin] -= Cumulant_term[5][type][species][cent][ptbin];
234  Cumulant_4[type][species][cent][ptbin] += 4*Cumulant_term[6][type][species][cent][ptbin];
235  Cumulant_4[type][species][cent][ptbin] += 4*Cumulant_term[7][type][species][cent][ptbin];
236  Cumulant_4[type][species][cent][ptbin] += 2*Cumulant_term[8][type][species][cent][ptbin];
237  Cumulant_4[type][species][cent][ptbin] += 2*Cumulant_term[9][type][species][cent][ptbin];
238  Cumulant_4[type][species][cent][ptbin] -= 6*Cumulant_term[10][type][species][cent][ptbin];
239 
241 // if (ptbin==ptCu && type==1) { // integrated cumulants
242 // cout << endl << "centrality= " << cent+1 << endl;
243 // double noAcc = Cumulant_term[0][type][species][cent][ptbin] -2*Cumulant_term[1][type][species][cent][ptbin];
244 // cout << "noAcc: " << setprecision(4) << noAcc << endl;
245 // cout << " 2: " << -Cumulant_term[2][type][species][cent][ptbin] << endl;
246 // cout << " 3: " << -Cumulant_term[3][type][species][cent][ptbin] << endl;
247 // cout << " 4: " << -2*Cumulant_term[4][type][species][cent][ptbin] << endl;
248 // cout << " 5: " << -Cumulant_term[5][type][species][cent][ptbin] << endl;
249 // cout << " 6: " << 4*Cumulant_term[6][type][species][cent][ptbin] << endl;
250 // cout << " 7: " << 4*Cumulant_term[7][type][species][cent][ptbin] << endl;
251 // cout << " 8: " << 2*Cumulant_term[8][type][species][cent][ptbin] << endl;
252 // cout << " 9: " << 2*Cumulant_term[9][type][species][cent][ptbin] << endl;
253 // cout << "10: " << -6*Cumulant_term[10][type][species][cent][ptbin] << endl;
254 // cout << "Acc/noAcc: " << Cumulant_4[type][species][cent][ptbin] / noAcc << endl;
255 // if(Cumulant_4[1][0][cent][ptCu]<0 && noAcc<0) {
256 // cout << "v2NoAcc: " << pow(-noAcc,0.25) << endl;
257 // cout << "v2Acc: " << pow(-Cumulant_4[1][0][cent][ptCu],0.25) << endl;
258 // cout << "v2Acc/noAcc: " << pow(-Cumulant_4[1][0][cent][ptCu],0.25) / pow(-noAcc,0.25) << endl;
259 // }
260 // }
262 
263  Cumulant_4_e[type][species][cent][ptbin] = Cumulant_term_e[0][type][species][cent][ptbin]**2;
264  Cumulant_4_e[type][species][cent][ptbin] += 2*(Cumulant_term_e[1][type][species][cent][ptbin])**2;
265  Cumulant_4_e[type][species][cent][ptbin] += 4*(Cumulant_term_e[6][type][species][cent][ptbin])**2;
266  Cumulant_4_e[type][species][cent][ptbin] += 4*(Cumulant_term_e[7][type][species][cent][ptbin])**2;
267  Cumulant_4_e[type][species][cent][ptbin] = sqrt(Cumulant_4_e[type][species][cent][ptbin]);
268 
269 
270  Cumulant_2[type][species][cent][ptbin] = Terms[1][type][COS][species][cent][ptbin];
271  Cumulant_2[type][species][cent][ptbin] -= Terms[3][type][COS][species][cent][ptbin]*Terms[7][type][COS][species][cent][ptbin];
272 
273  Cumulant_2_e[type][species][cent][ptbin] = Terms_e[1][type][COS][species][cent][ptbin]**2;
274  Cumulant_2_e[type][species][cent][ptbin] += (Terms_e[3][type][COS][species][cent][ptbin]*Terms[7][type][COS][species][cent][ptbin])**2;
275  Cumulant_2_e[type][species][cent][ptbin] += (Terms[3][type][COS][species][cent][ptbin]*Terms_e[7][type][COS][species][cent][ptbin])**2;
276  Cumulant_2_e[type][species][cent][ptbin] = sqrt(Cumulant_2_e[type][species][cent][ptbin]);
277 
278 
279  }
280  }
281  }
282  }
283 
285  // each pt bin is calculated for all centralities, then the next pt bin
286  int ptcount[SPECIES]={0}; // index of pt bins
287  double used_pt_yields4[CENTBINS+1] = {0.};
288  double used_pt_yields2[CENTBINS+1] = {0.};
289  double v2IntCu_2[CENTBINS] = {0.};
290  double v2IntCu_2Err[CENTBINS] = {0.};
291  double v2IntCu_4[CENTBINS] = {0.};
292  double v2IntCu_4Err[CENTBINS] = {0.};
293 
294  // This part reconstructs the elliptic flow from the Cumulant terms
295  for(int species=0; species<SPECIES; species++){
296 
297  cout << endl << "cent: \t v2IntCu_2 +/- err,\t v2IntCu_4 +/- err" << endl;
298 
299  for(int ptbin=0; ptbin<PTBINS; ptbin++){
300 
301  // defaults for min bias set to an unphysical result
302  v2_2[species][0][ptcount[species]] = 10.;
303  v2_2_e[species][0][ptcount[species]] = 10.;
304  v2_4[species][0][ptcount[species]] = 10.;
305  v2_4_e[species][0][ptcount[species]] = 10.;
306  used_pt_yields4[0] = 0.; // reset for each pt and species
307  used_pt_yields2[0] = 0.;
308 
309  for(int cent=0; cent<CENTBINS; cent++){
310  if (ptbin==ptCu) { // integrated cumulants
311  if(Cumulant_2[1][0][cent][ptCu]>0) {
312  v2IntCu_2[cent] = pow(Cumulant_2[1][0][cent][ptCu],0.5);
313  v2IntCu_2Err[cent] = pow((.5/(Cumulant_2[1][0][cent][ptCu])),0.5)*Cumulant_2_e[1][0][cent][ptCu];
314  }
315  if(Cumulant_4[1][0][cent][ptCu]<0) {
316  v2IntCu_4[cent] = pow(-Cumulant_4[1][0][cent][ptCu],0.25);
317  v2IntCu_4Err[cent] = pow((.25/(-Cumulant_4[1][0][cent][ptCu])),0.75)*Cumulant_4_e[1][0][cent][ptCu];
318  }
319  cout << cent+1 << ": \t" << setprecision(3) << v2IntCu_2[cent]*100. <<
320  " \t+/- " << v2IntCu_2Err[cent]*100. << ", \t" << v2IntCu_4[cent]*100. <<
321  " \t+/- " << v2IntCu_4Err[cent]*100. << endl;
322  }
323 
324  // defaults for each cent set to an unphysical result
325  v2_2[species][cent+1][ptcount[species]] = 10.;
326  v2_2_e[species][cent+1][ptcount[species]] = 10.;
327  v2_4[species][cent+1][ptcount[species]] = 10.;
328  v2_4_e[species][cent+1][ptcount[species]] = 10.;
329  used_pt_yields4[cent+1] = 0.;
330  used_pt_yields2[cent+1] = 0.;
331 
332  // 4-particle cumulant v2
333  if(Cumulant_4[1][species][cent][ptbin] < 0 && Cumulant_4[0][species][cent][ptbin] < 0) { // both must be negative
334  v2_4[species][cent+1][ptcount[species]] += -cent_yields[species][cent][ptbin]*Cumulant_4[0][species][cent][ptbin]/(pow(-Cumulant_4[1][species][cent][ptbin],.75));
335  v2_4[species][0][ptcount[species]] += v2_4[species][cent+1][ptcount[species]];
336  if(ptbin<PTINTBINS) { // pt integration
337  v2IntPt_4[species][cent+1] += v2_4[species][cent+1][ptcount[species]];
338  usedInt_cent_yields4[species][cent+1] += cent_yields[species][cent][ptbin];
339  }
340  v2_4_e[species][cent+1][ptcount[species]] += pow(cent_yields[species][cent][ptbin]*Cumulant_4_e[0][species][cent][ptbin]/(pow(-Cumulant_4[1][species][cent][ptbin],.75)),2);
341  v2_4_e[species][cent+1][ptcount[species]] += pow(.75*cent_yields[species][cent][ptbin]*Cumulant_4[0][species][cent][ptbin]*Cumulant_4_e[1][species][cent][ptbin]/(pow(-Cumulant_4[1][species][cent][ptbin],1.75)),2);
342  v2_4_e[species][0][ptcount[species]] += v2_4_e[species][cent+1][ptcount[species]];
343  if(ptbin<PTINTBINS) v2IntPt_4_e[species][cent+1] += v2_4_e[species][cent+1][ptcount[species]];
344 
345  used_pt_yields4[cent+1] = cent_yields[species][cent][ptbin];
346  used_cent_yields4[species][cent+1] += cent_yields[species][cent][ptbin];
347  used_pt_yields4[0] += used_pt_yields4[cent+1]; // min bias
348  }//if
349 
350  // 2-particle cumulant v2
351  if(Cumulant_2[1][species][cent][ptbin] > 0) { // must be positive
352  v2_2[species][cent+1][ptcount[species]] += cent_yields[species][cent][ptbin]*Cumulant_2[0][species][cent][ptbin]/(pow(Cumulant_2[1][species][cent][ptbin],.5));
353  v2_2[species][0][ptcount[species]] += v2_2[species][cent+1][ptcount[species]];
354  if(ptbin<PTINTBINS) { // pt integration
355  v2IntPt_2[species][cent+1] += v2_2[species][cent+1][ptcount[species]];
356  usedInt_cent_yields2[species][cent+1] += cent_yields[species][cent][ptbin];
357  }
358  v2_2_e[species][cent+1][ptcount[species]] += pow(cent_yields[species][cent][ptbin]*Cumulant_2_e[0][species][cent][ptbin]/(pow(Cumulant_2[1][species][cent][ptbin],.5)),2);
359  v2_2_e[species][cent+1][ptcount[species]] += pow(.5*cent_yields[species][cent][ptbin]*Cumulant_2[0][species][cent][ptbin]*Cumulant_2_e[1][species][cent][ptbin]/(pow(Cumulant_2[1][species][cent][ptbin],1.5)),2);
360  v2_2_e[species][0][ptcount[species]] += v2_2_e[species][cent+1][ptcount[species]];
361  if(ptbin<PTINTBINS) v2IntPt_2_e[species][cent+1] += v2_2_e[species][cent+1][ptcount[species]];
362 
363  used_pt_yields2[cent+1] = cent_yields[species][cent][ptbin];
364  used_cent_yields2[species][cent+1] += cent_yields[species][cent][ptbin];
365  used_pt_yields2[0] += used_pt_yields2[cent+1]; // min bias
366  }//if
367 
368  }//cent
369 
370  for(int cent=0; cent<CENTBINS+1; cent++){ // includes cent=0
371 
372  if(used_pt_yields2[cent] > 0){
373  v2_2[species][cent][ptcount[species]] /= used_pt_yields2[cent];
374  v2_2_e[species][cent][ptcount[species]] = sqrt(v2_2_e[species][cent][ptcount[species]]);
375  v2_2_e[species][cent][ptcount[species]] /= used_pt_yields2[cent];
376  }
377  if(used_pt_yields4[cent] > 0){
378  v2_4[species][cent][ptcount[species]] /= used_pt_yields4[cent];
379  v2_4_e[species][cent][ptcount[species]] = sqrt(v2_4_e[species][cent][ptcount[species]]);
380  v2_4_e[species][cent][ptcount[species]] /= used_pt_yields4[cent];
381  }
382 
383  }//cent
384 
385  ptcount[species]++;
386  }//pt
387 
388  // For pt integrated v2
389  for(int cent=0; cent<CENTBINS+1; cent++){ // includes cent=0
390  if(usedInt_cent_yields2[species][cent] > 0){
391  v2IntPt_2[species][cent] /= usedInt_cent_yields2[species][cent];
392  v2IntPt_2_e[species][cent] = sqrt(v2IntPt_2_e[species][cent]);
393  v2IntPt_2_e[species][cent] /= usedInt_cent_yields2[species][cent];
394  }
395  if(usedInt_cent_yields4[species][cent] > 0){
396  v2IntPt_4[species][cent] /= usedInt_cent_yields4[species][cent];
397  v2IntPt_4_e[species][cent] = sqrt(v2IntPt_4_e[species][cent]);
398  v2IntPt_4_e[species][cent] /= usedInt_cent_yields4[species][cent];
399  }
400  }//cent
401 
402  }//species
403 
404  cout << endl;
405  for(int ptbin=0; ptbin<PTBINS; ptbin++){ // for checking
406  //cout << ptbin << ": " << v2_2[0][3][ptbin] << " +/- " << v2_2_e[0][3][ptbin] << endl;
407  //cout << ptbin << ": " << cent_yields[0][3][ptbin] << endl;
408  //cout << ptbin << ": " << Cumulant_2[1][0][5][ptbin] << endl;
409  }
410 
411  cout << endl << "cent: \t v2IntPt_2 +/- err,\t v2IntPt_4 +/- err" << endl;
412  for(int cent=0; cent<CENTBINS+1; cent++){
413 
414  //Throw away last point since its middle-pt value is not well defined.
415  //It represents all particles with pt > 6 GeV/c
416  //cout << ptOverflow << ": " << v2_2[0][cent][ptOverflow] << endl;
417  v2_2[0][cent][ptOverflow]= 10.;
418  v2_4[0][cent][ptOverflow]= 10.;
419  //Throw away next-to-last point since it contains the integrated cumulant.
420  v2_2[0][cent][ptCu]= 10.;
421  v2_4[0][cent][ptCu]= 10.;
422  //cout << ptOverflow << ": " << v2_2[0][cent][ptOverflow] << endl; // pt = 10?
423 
424  // output integrated over pt
425  cout << cent << ": \t" << setprecision(3) << v2IntPt_2[0][cent]*100. << " \t+/- " <<
426  v2IntPt_2_e[0][cent]*100. << ", \t" << v2IntPt_4[0][cent]*100. << " \t+/- " <<
427  v2IntPt_4_e[0][cent]*100. << endl;
428 
429  }//cent
430  cout << endl;
431 
433  // plot the results
434 
435  TFile graphFile("flow.dirCumulant.graphs.root", "RECREATE"); // output
436 
437  // For v2(pt)
438 
439  double middle_pt_points[PTBINS]={0};
440  for(int i=0; i<PTBINS; i++){
441  middle_pt_points[i]=double((i)/10. + 0.1); // GeV = bin/10.
442  }
443  double middle_pt_points_e[PTBINS]={0}; // no errors
444 
445  float v2max = 0.4;
446 
447  TGraphErrors *charged_v2_2[CENTBINS+1];
448  TGraphErrors *charged_v2_4[CENTBINS+1];
449  TCanvas* can[CENTBINS+1];
450  int canvasWidth = 780, canvasHeight = 600; // landscape
451  for(int cent=0; cent<CENTBINS+1; cent++){
452  TString* canName = new TString("centrality_");
453  *canName += (cent);
454  if(ptGraphs) {
455  cout << "plot " << canName->Data() << endl;
456  can[cent] = new TCanvas(canName->Data(), canName->Data(), canvasWidth, canvasHeight);
457  }
458  TLegend *legend = new TLegend(.1,.7,.46,.9,NULL,"brNDC");
459  legend->SetFillColor(kWhite);
460  TPaveLabel* title = new TPaveLabel(0.7,0.96,0.9,0.99,canName->Data());
461  if(ptGraphs) title->Draw();
462 
463  TString graphName("v22_");
464  graphName += (cent);
465  charged_v2_2[cent]= new TGraphErrors(PTBINS,middle_pt_points,v2_2[0][cent],middle_pt_points_e,v2_2_e[0][cent]);
466  charged_v2_2[cent]->SetMarkerStyle(20);
467  charged_v2_2[cent]->SetMarkerColor(2);
468  charged_v2_2[cent]->SetLineColor(2);
469  charged_v2_2[cent]->SetMinimum(0.);
470  charged_v2_2[cent]->SetMaximum(v2max);
471  charged_v2_2[cent]->GetXaxis()->SetTitle("P_{t} (GeV/c)");
472  charged_v2_2[cent]->GetYaxis()->SetTitle("v_{2}");
473  charged_v2_2[cent]->SetTitle("Elliptic Flow");
474  charged_v2_2[cent]->Write(graphName.Data());
475 
476  TString graphName("v24_");
477  graphName += (cent);
478  charged_v2_4[cent]= new TGraphErrors(PTBINS,middle_pt_points,v2_4[0][cent],middle_pt_points_e,v2_4_e[0][cent]);
479  charged_v2_4[cent]->SetMarkerStyle(20);
480  charged_v2_4[cent]->SetMarkerColor(4);
481  charged_v2_4[cent]->SetLineColor(4);
482  charged_v2_4[cent]->SetMinimum(0.);
483  charged_v2_4[cent]->SetMaximum(v2max);
484  charged_v2_4[cent]->GetXaxis()->SetTitle("P_{t} (GeV/c)");
485  charged_v2_4[cent]->GetYaxis()->SetTitle("v_{2}{4}");
486  charged_v2_4[cent]->SetTitle("Elliptic Flow");
487  charged_v2_4[cent]->Write(graphName.Data());
488 
489  if(ptGraphs) charged_v2_2[cent]->Draw("AP");
490  legend->AddEntry(charged_v2_2[cent],"charged hadron v_{2}{2}","p");
491  if(ptGraphs) charged_v2_4[cent]->Draw("P");
492  legend->AddEntry(charged_v2_4[cent],"charged hadron v_{2}{4}","p");
493  if(ptGraphs) legend->Draw("same");
494 
495  //can[cent]->Print(".pdf");
496 
497  }
498 
499  //Event_counter histograms
500 
501  TH1* Event_counter = dynamic_cast<TH1*>(myfile->Get("Event_counter"));
502  if (!Event_counter) {
503  cout << "### Can't find hist Event_counter" << endl;
504  } else {
505  // make the graph page
506  gStyle->SetOptStat("e");
507  TString* canEvtName = new TString("Event_counter");
508  cout << "plot " << canEvtName->Data() << endl;
509  TCanvas* canEvt = new TCanvas(canEvtName->Data(), canEvtName->Data(), 780, 600);
510  Event_counter->SetMinimum(0.);
511  Event_counter->GetXaxis()->SetTitle("centrality bin");
512  Event_counter->Draw();
513  Event_counter->Write(canEvtName->Data());
514  }
515  TH1* Event_counterWeighted = dynamic_cast<TH1*>(myfile->Get("Event_counterWeighted"));
516  if (!Event_counterWeighted) {
517  cout << "### Can't find hist Event_counterWeighted" << endl;
518  } else {
519  // make the graph page
520  //gStyle->SetOptStat("e");
521  TString* canEvtWgtName = new TString("Event_counterWeighted");
522  cout << "plot " << canEvtWgtName->Data() << endl;
523  TCanvas* canEvtWgt = new TCanvas(canEvtWgtName->Data(), canEvtWgtName->Data(), 780, 600);
524  Event_counterWeighted->SetMinimum(0.);
525  Event_counterWeighted->GetXaxis()->SetTitle("centrality bin");
526  Event_counterWeighted->Draw("hist");
527  Event_counterWeighted->Write(canEvtWgtName->Data());
528  }
529  //gStyle->SetOptStat(0);
530 
531  // v2 integrated
532 
533  float v2Intmax = 8.;
534 
535  // make the graph page
536  TString* canIntName = new TString("v2Int");
537  cout << "plot " << canIntName->Data() << endl;
538  TCanvas* canInt = new TCanvas(canIntName->Data(), canIntName->Data(), 600, 780);
539  TLegend* legendInt = new TLegend(0.35,0.1,0.6,0.35);
540  legendInt->SetFillColor(kWhite);
541 
542  // make a histogram
543  TString* histCenName = new TString("v2Int");
544  TH1F* histCen = new TH1F(histCenName->Data(), histCenName->Data(), 80, 0., 80.);
545  histCen->SetLineColor(kBlack);
546  histCen->Draw();
547 
548  TGraphErrors *grv2IntPt_2 = new TGraphErrors(CENTBINS);
549  TGraphErrors *grv2IntPt_4 = new TGraphErrors(CENTBINS);
550  TGraphErrors *grv2IntCu_2 = new TGraphErrors(CENTBINS);
551  TGraphErrors *grv2IntCu_4 = new TGraphErrors(CENTBINS);
552 
553  int m = 0;
554  for(int cent=0; cent<CENTBINS; cent++){
555  grv2IntPt_2->SetPoint(m, x[CENTBINS-cent-1], v2IntPt_2[0][cent+1]*100.);
556  grv2IntPt_2->SetPointError(m, 0., v2IntPt_2_e[0][cent+1]*100.);
557  grv2IntPt_4->SetPoint(m, x[CENTBINS-cent-1], v2IntPt_4[0][cent+1]*100.);
558  grv2IntPt_4->SetPointError(m, 0., v2IntPt_4_e[0][cent+1]*100.);
559  grv2IntCu_2->SetPoint(m, x[CENTBINS-cent-1], v2IntCu_2[cent]*100.);
560  grv2IntCu_2->SetPointError(m, 0., v2IntCu_2Err[cent]*100.);
561  grv2IntCu_4->SetPoint(m, x[CENTBINS-cent-1], v2IntCu_4[cent]*100.);
562  grv2IntCu_4->SetPointError(m, 0., v2IntCu_4Err[cent]*100.);
563 
564  m++;
565  }//cent
566 
567  grv2IntPt_2->SetMarkerStyle(kFullCircle);
568  grv2IntPt_2->SetMarkerColor(kRed);
569  grv2IntPt_2->SetLineColor(kRed);
570  grv2IntPt_2->SetMinimum(0.);
571  grv2IntPt_2->SetMaximum(v2Intmax);
572  grv2IntPt_2->SetTitle("Elliptic Flow");
573  grv2IntPt_2->Draw("AP");
574  legendInt->AddEntry(grv2IntPt_2,"v_{2}{2}(pt)","p");
575 
576  grv2IntPt_4->SetMarkerStyle(kFullSquare);
577  grv2IntPt_4->SetMarkerColor(kBlue);
578  grv2IntPt_4->SetLineColor(kBlue);
579  grv2IntPt_4->Draw("P");
580  legendInt->AddEntry(grv2IntPt_4,"v_{2}{4}(pt)","p");
581 
582  grv2IntCu_2->SetMarkerStyle(kOpenCircle);
583  grv2IntCu_2->SetMarkerColor(kRed);
584  grv2IntCu_2->SetLineColor(kRed);
585  grv2IntCu_2->Draw("P");
586  legendInt->AddEntry(grv2IntCu_2,"v_{2}{2}(cu)","p");
587 
588  grv2IntCu_4->SetMarkerStyle(kOpenSquare);
589  grv2IntCu_4->SetMarkerColor(kBlue);
590  grv2IntCu_4->SetLineColor(kBlue);
591  grv2IntCu_4->Draw("P");
592  legendInt->AddEntry(grv2IntCu_4,"v_{2}{4}(cu)","p");
593 
594  legendInt->Draw("same");
595  // output
596  grv2IntPt_2->Write("v22IntPt");
597  grv2IntPt_4->Write("v24IntPt");
598  grv2IntCu_2->Write("v22IntCu");
599  grv2IntCu_4->Write("v24IntCu");
600 
601  //can[cent]->Print(".pdf");
602 
603  // label the axes
604  TLatex l;
605  l.SetNDC();
606  l.SetTextSize(0.06);
607  l.DrawLatex(0.5,0.02,"% Most Central");
608  l.SetTextAngle(90);
609  l.DrawLatex(0.05,0.7,"v_{2} (%)" );
610 
611  graphFile.Close();
612 
613 }
Generated by   doxygen 1.8.5