StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
calculateCumulant.C
1 #include "float.h"
2 #include "calculateCumulant.h"
3 
4 void calculateCumulant(const char* histFileName){
5 
6  // gSystem->Load("myProjectProfile_C.so"); //user defined integration
7 
8  double r0Sq=r0*r0;
9  double perCent=0.01;
10 
11 
12  TString xLabel = "Pseudorapidity";
13  if (isPidFlow) { xLabel = "Rapidity"; }
14 
15 
16  TFile* f = new TFile(histFileName,"UPDATE");
17 
18  struct histFullHars {
19 
20  TProfile2D** mHistCumul2D;
21  TProfile** mHistCumulEta;
22  TProfile** mHistCumulPt;
23 
24  TProfile2D* mHistCumulMix2D;
25  TProfile* mHistCumulMixEta;
26  TProfile* mHistCumulMixPt;
27 
28  TH2D** mHist_v2D;
29  TH1D** mHist_vEta;
30  TH1D** mHist_vPt;
31 
32 
33 
34  TH2D* mHistMix_v2D;
35  TH1D* mHistMix_vEta;
36  TH1D* mHistMix_vPt;
37 
38 
39 
40 
41  Double_t mCumulIntegG0[nCumulIntegOrders*nCumulInteg_qMax];
42  TH1D* mHistMultSum; //histo for holding mMultSum.
43  TH1D* mHistNEvent;
44 
45  TH1D* mHistMeanWgtSqrSum; // sum <w^2> over evts.
46 
47  Double_t mCumulIntegG0Mix[nCumulMixHar_pMax*nCumulMixHar_qMax];
48 
49  TH1D** mHistCumulIntegG0Sum; //summation of G value.
50  TH1D** mHistCumulIntegG0MixSum; //summation of G value.
51 
52  };
53 
54  struct histFulls {
55  TProfile** mHistCumul; //integrated cumulant from differential
56  TH1D** mHist_v; //integrated flow from differential
57 
58  TProfile* mHistCumulMix; //cumulant for mix 1st and 2nd hars.
59  TH1D* mHistMix_v; //integrated flow from differential
60 
61  struct histFullHars histFullHar[nHars];
62  };
63 
64  struct histFulls histFull[nSels];
65 
66 
67  TString* histTitle;
68  TString* histTitle2;
69 
70 
71  //================ fill in pointers ===================
72  //for each selection
73  for (int k = 0; k < nSels; k++) {
74  char countSels[2];
75  sprintf(countSels,"%d",k+1);
76 
77  histFull[k].mHistCumul = new TProfile*[nCumulDiffOrders];
78  histFull[k].mHist_v = new TH1D*[nCumulDiffOrders];
79 
80  for (int ord = 0; ord < nCumulDiffOrders; ord++) {
81  char theCumulOrderChar[2]; // if >10, need to use char*
82  sprintf(theCumulOrderChar,"%d",(ord+1)*2);
83  histTitle = new TString("Flow_Cumul_Order");
84  histTitle->Append(*theCumulOrderChar);
85  histTitle->Append("_Sel");
86  histTitle->Append(*countSels);
87  histFull[k].mHistCumul[ord] = (TProfile *)f->Get(histTitle->Data());
88  delete histTitle;
89 
90  histTitle = new TString("Flow_Cumul_v_Order");
91  histTitle->Append(*theCumulOrderChar);
92  histTitle->Append("_Sel");
93  histTitle->Append(*countSels);
94  histFull[k].mHist_v[ord] =
95  new TH1D(*(histFull[k].mHistCumul[ord]->ProjectionX(histTitle->Data(),"e")));
96  histFull[k].mHist_v[ord]->SetTitle(histTitle->Data());
97  histFull[k].mHist_v[ord]->SetXTitle("harmonic");
98  histFull[k].mHist_v[ord]->SetYTitle("v (%)");
99  if (ord>0) histFull[k].mHist_v[ord]->Scale(1./profScale);
100  delete histTitle;
101 
102  }
103 
104 
106 
107  histTitle = new TString("Flow_CumulMix");
108  histTitle->Append("_Sel");
109  histTitle->Append(*countSels);
110  histFull[k].mHistCumulMix = (TProfile *)f->Get(histTitle->Data());
111  delete histTitle;
112 
113 
114  histTitle = new TString("Flow_CumulMix_v_Sel");
115  histTitle->Append(*countSels);
116  histFull[k].mHistMix_v =
117  new TH1D(*(histFull[k].mHistCumulMix->ProjectionX(histTitle->Data(),"e")));
118  histFull[k].mHistMix_v->SetTitle(histTitle->Data());
119  histFull[k].mHistMix_v->SetXTitle("place for v1 from mixed harmonic");
120  histFull[k].mHistMix_v->SetYTitle("v (%)");
121  histFull[k].mHistMix_v->Scale(1./profScale);
122  delete histTitle;
124 
125 
126 
127 
128  // for each harmonic
129  for (int j = 0; j < nHars; j++) {
130  char countHars[2];
131  sprintf(countHars,"%d",j+1);
132 
133  histFull[k].histFullHar[j].mHistCumul2D =
134  new TProfile2D*[nCumulDiffOrders];
135  histFull[k].histFullHar[j].mHistCumulEta =
136  new TProfile*[nCumulDiffOrders];
137  histFull[k].histFullHar[j].mHistCumulPt =
138  new TProfile*[nCumulDiffOrders];
139 
140  histFull[k].histFullHar[j].mHist_v2D =
141  new TH2D*[nCumulDiffOrders];
142  histFull[k].histFullHar[j].mHist_vEta =
143  new TH1D*[nCumulDiffOrders];
144  histFull[k].histFullHar[j].mHist_vPt =
145  new TH1D*[nCumulDiffOrders];
146 
147 
149  if (j==0){ //only j==0 makes sense here.
150 
151  histTitle = new TString("Flow_CumulMix2D");
152  histTitle->Append("_Sel");
153  histTitle->Append(*countSels);
154  histTitle->Append("_Har");
155  histTitle->Append(*countHars);
156  histFull[k].histFullHar[j].mHistCumulMix2D =
157  (TProfile2D *)f->Get(histTitle->Data());
158  delete histTitle;
159 
160 
161 
162  histTitle = new TString("Flow_CumulMixEta");
163  histTitle->Append("_Sel");
164  histTitle->Append(*countSels);
165  histTitle->Append("_Har");
166  histTitle->Append(*countHars);
167  histFull[k].histFullHar[j].mHistCumulMixEta =
168  (TProfile *)f->Get(histTitle->Data());
169  delete histTitle;
170 
171  histTitle = new TString("Flow_CumulMixPt");
172  histTitle->Append("_Sel");
173  histTitle->Append(*countSels);
174  histTitle->Append("_Har");
175  histTitle->Append(*countHars);
176  histFull[k].histFullHar[j].mHistCumulMixPt =
177  (TProfile *)f->Get(histTitle->Data());
178  delete histTitle;
179 
180 
181 
182  histFull[k].histFullHar[j].mHistCumulIntegG0MixSum =
183  new TH1D*[nCumulMixHar_pMax*nCumulMixHar_qMax];
184 
185 
186  for (int pq = 0; pq < nCumulMixHar_pMax*nCumulMixHar_qMax ; pq++) {
187  int pIndex = (pq/nCumulMixHar_qMax);
188  int qIndex = pq%nCumulMixHar_qMax;
189 
190  char pIndexChar[2];
191  char qIndexChar[2]; // if >10, need to use char*
192  sprintf(pIndexChar,"%d",pIndex);
193  sprintf(qIndexChar,"%d",qIndex);
194 
195  histTitle = new TString("Flow_CumulIntegG0MixSum");
196  histTitle->Append("_GenFunIdxp");
197  histTitle->Append(*pIndexChar);
198  histTitle->Append("_GenFunIdxq");
199  histTitle->Append(*qIndexChar);
200  histTitle->Append("_Sel");
201  histTitle->Append(*countSels);
202  histTitle->Append("_Har");
203  histTitle->Append(*countHars);
204  histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq]=
205  (TH1D *)f->Get(histTitle->Data());
206  delete histTitle;
207 
208  histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] = 0.;
209  }
210  }
211 
212 
214 
215 
216 
217  // For each cumulant order
218  for (int ord = 0; ord < nCumulDiffOrders; ord++) {
219  char theCumulOrderChar[2]; // if >10, need to use char*
220  sprintf(theCumulOrderChar,"%d",(ord+1)*2);
221 
222  histTitle = new TString("Flow_Cumul2D_Order");
223  histTitle->Append(*theCumulOrderChar);
224  histTitle->Append("_Sel");
225  histTitle->Append(*countSels);
226  histTitle->Append("_Har");
227  histTitle->Append(*countHars);
228  histFull[k].histFullHar[j].mHistCumul2D[ord] =
229  (TProfile2D *)f->Get(histTitle->Data());
230  delete histTitle;
231 
232  histTitle = new TString("Flow_CumulEta_Order");
233  histTitle->Append(*theCumulOrderChar);
234  histTitle->Append("_Sel");
235  histTitle->Append(*countSels);
236  histTitle->Append("_Har");
237  histTitle->Append(*countHars);
238  histFull[k].histFullHar[j].mHistCumulEta[ord]=
239  (TProfile *)f->Get(histTitle->Data());
240  delete histTitle;
241 
242 
243  histTitle = new TString("Flow_CumulPt_Order");
244  histTitle->Append(*theCumulOrderChar);
245  histTitle->Append("_Sel");
246  histTitle->Append(*countSels);
247  histTitle->Append("_Har");
248  histTitle->Append(*countHars);
249  histFull[k].histFullHar[j].mHistCumulPt[ord] =
250  (TProfile *)f->Get(histTitle->Data());
251  delete histTitle;
252 
253  }
254 
255 
256  histFull[k].histFullHar[j].mHistCumulIntegG0Sum =
257  new TH1D*[nCumulIntegOrders*nCumulInteg_qMax];
258 
259  for (int pq = 0; pq < nCumulIntegOrders*nCumulInteg_qMax; pq++) {
260  int cumulIndex = (pq/nCumulInteg_qMax) + 1; // like 1,2,3.
261  //not begining with 0. That is "p" in Eq. (B1, PG5).
262  int qIndex = pq%(nCumulInteg_qMax); // like 0,1,..qMax-1
263  //begining with 0. Just like Eq. (B3, PG5).
264 
265  char theCumulOrderChar[2];
266  char qIndexOrderChar[2]; // if >10, need to use char*
267  sprintf(theCumulOrderChar,"%d",cumulIndex*2);
268  sprintf(qIndexOrderChar,"%d",qIndex);
269 
270  histTitle = new TString("Flow_CumulIntegG0Sum_Order");
271  histTitle->Append(*theCumulOrderChar);
272  histTitle->Append("_GenFunIdx");
273  histTitle->Append(*qIndexOrderChar);
274  histTitle->Append("_Sel");
275  histTitle->Append(*countSels);
276  histTitle->Append("_Har");
277  histTitle->Append(*countHars);
278  histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq] =
279  (TH1D *)f->Get(histTitle->Data());
280  delete histTitle;
281 
282  histFull[k].histFullHar[j].mCumulIntegG0[pq] = 0.;
283 
284  }
285 
286 
287 
288 
289  histTitle = new TString("Flow_CumulMultSum_Sel");
290  histTitle->Append(*countSels);
291  histTitle->Append("_Har");
292  histTitle->Append(*countHars);
293  histFull[k].histFullHar[j].mHistMultSum =
294  (TH1D *)f->Get(histTitle->Data());
295  delete histTitle;
296 
297 
298  histTitle = new TString("Flow_CumulNEvent_Sel");
299  histTitle->Append(*countSels);
300  histTitle->Append("_Har");
301  histTitle->Append(*countHars);
302  histFull[k].histFullHar[j].mHistNEvent =
303  (TH1D *)f->Get(histTitle->Data());
304  delete histTitle;
305 
306 
307  histTitle = new TString("Flow_CumulMeanWgtSqrSum_Sel");
308  histTitle->Append(*countSels);
309  histTitle->Append("_Har");
310  histTitle->Append(*countHars);
311  histFull[k].histFullHar[j].mHistMeanWgtSqrSum =
312  (TH1D *)f->Get(histTitle->Data());
313  delete histTitle;
314 
315 
316  }
317  }
318 
319  //================ end of fill in pointers ===================
320 
321 
322 
323 
324  for (int k = 0; k < nSels; k++) {
325 
326  char countSels[2];
327  sprintf(countSels,"%d",k+1);
328 
329  double meanIntegV[nHars]; // V**1
330  double meanIntegV2[nHars]; // V**2
331  double meanIntegV3[nHars]; // V**3
332  double meanIntegV4[nHars]; // V**4
333  double cumulInteg1[nHars]; // outside of harmonic loop
334  double cumulInteg2[nHars];
335  double cumulInteg3[nHars];
336  double q2[nHars]; // for old method. <Q>**2 in (74) of old paper.
337  double q4[nHars];
338  double q6[nHars];
339 
340 
341  double cumulIntegMix[nHars];
342  double meanIntegVMix[nHars];
343 
344  for (int j = 0; j < nHars; j++) {
345  meanIntegV[j] = 0.;
346  meanIntegV2[j] = 0.;
347  meanIntegV3[j] = 0.;
348  meanIntegV4[j] = 0.;
349  cumulInteg1[j] = 0.;
350  cumulInteg2[j] = 0.;
351  cumulInteg3[j] = 0.;
352  cumulIntegMix[j] = 0.;
353  meanIntegVMix[j] = 0.;
354  }
355 
356  for (int j = 0; j < nHars; j++) {
357 
358 
359  char countHars[2];
360  sprintf(countHars,"%d",j+1);
361 
362  cout<<" ========== Sel"<<k+1<<" Har"<<j+1<<" ==========="<<endl;
363 
364  cout<<" event # "<<histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1)<<endl;
365 
366  double mAvMult =
367  histFull[k].histFullHar[j].mHistMultSum->GetBinContent(1)/
368  histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
369 
370  cout<<"mAvMult Sel"<<k<<" har"<<j<<" "<<mAvMult<<endl;
371 
372 
373  double CpInteg[nCumulIntegOrders]; // Cp in (B4, PG6)
374 
375  for (int pq = 0; pq < nCumulIntegOrders; pq ++) CpInteg[pq] = 0.;
376  for (int pq = 0; pq < nCumulIntegOrders*nCumulInteg_qMax; pq++) {
377 
378  int theCumulOrder = (pq/nCumulInteg_qMax) + 1; // like 1,2,3.
379 
380  histFull[k].histFullHar[j].mCumulIntegG0[pq] =
381  histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq]->GetBinContent(1)/
382  histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
383 
384  // cout<<" --- histFull[k].histFullHar[j].mCumulIntegG0[pq] "<<histFull[k].histFullHar[j].mCumulIntegG0[pq]<<endl;
385 
386  if (!isNewMethod){
387  CpInteg[theCumulOrder-1] +=
388  (log(histFull[k].histFullHar[j].mCumulIntegG0[pq]) /
389  ((float)nCumulInteg_qMax));
390  } else {
391  CpInteg[theCumulOrder-1] +=
392  (mAvMult*(pow(histFull[k].histFullHar[j].mCumulIntegG0[pq], 1./mAvMult) -1.) /
393  float(nCumulInteg_qMax)); // (B3, PG6)
394  // cout<<" CpInteg[theCumulOrder-1] "<<CpInteg[theCumulOrder-1]<<endl;
395  }
396  }
397 
398  cumulInteg1[j] = // (B5, PG7)
399  (3.*CpInteg[1-1] - (3./2.)*CpInteg[2-1] + (1./3.)*CpInteg[3-1]) / r0Sq;
400 
401  cumulInteg2[j] = ((-10.*CpInteg[1-1]) + (8.*CpInteg[2-1]) -
402  (2.*CpInteg[3-1])) / (r0Sq*r0Sq);
403 
404  cumulInteg3[j] = ( (18.*CpInteg[1-1]) - (18.*CpInteg[2-1]) + (6.*CpInteg[3-1]))
405  / (r0Sq*r0Sq*r0Sq);
406 
407 
408  for (int ord = 0; ord < nCumulDiffOrders; ord++) {
409  char theCumulOrderChar[2]; // if >10, need to use char*
410  sprintf(theCumulOrderChar,"%d",(ord+1)*2);
411 
412 
413 
414  histTitle = new TString("Flow_Cumul_vEta_Order");
415  histTitle->Append(*theCumulOrderChar);
416  histTitle->Append("_Sel");
417  histTitle->Append(*countSels);
418  histTitle->Append("_Har");
419  histTitle->Append(*countHars);
420  histFull[k].histFullHar[j].mHist_vEta[ord] =
421  new TH1D(*(histFull[k].histFullHar[j].mHistCumulEta[ord]->
422  ProjectionX(histTitle->Data(),"e")));
423  histFull[k].histFullHar[j].mHist_vEta[ord]->SetTitle(histTitle->Data());
424  histFull[k].histFullHar[j].mHist_vEta[ord]->SetXTitle((char*)xLabel.Data());
425  histFull[k].histFullHar[j].mHist_vEta[ord]->SetYTitle("v (%)");
426  if (ord>0) histFull[k].histFullHar[j].mHist_vEta[ord]->Scale(1./profScale);
427  delete histTitle;
428 
429  histTitle = new TString("Flow_Cumul_vPt_Order");
430  histTitle->Append(*theCumulOrderChar);
431  histTitle->Append("_Sel");
432  histTitle->Append(*countSels);
433  histTitle->Append("_Har");
434  histTitle->Append(*countHars);
435 
436  // if you want integrate for only TPC/FTPC, you need to do
437  // 1. load myProjectProfile.C at the begining
438  // 2. append "_TPC" or "_FTPC" here to the title
439  // 3. use myProjectionY (X) to get mHist_vPt or mHist_vEta.
440  // histTitle->Append("_TPC");
441 
442  histTitle2 = new TString("Flow_Cumul2D_Order");
443  histTitle2->Append(*theCumulOrderChar);
444  histTitle2->Append("_Sel");
445  histTitle2->Append(*countSels);
446  histTitle2->Append("_Har");
447  histTitle2->Append(*countHars);
448 
449  histFull[k].histFullHar[j].mHist_vPt[ord] =
450  //uncomment this line if you want to integrate over different range
451  //myProjectionY((TProfile2D *)f->Get(histTitle2->Data()) , 0., 1., 1, ((j==0 || j== 2) ? 1 : 0) );
452  new TH1D(*(histFull[k].histFullHar[j].mHistCumulPt[ord]->ProjectionX(histTitle->Data(),"e")));
453  histFull[k].histFullHar[j].mHist_vPt[ord]->SetTitle(histTitle->Data());
454  histFull[k].histFullHar[j].mHist_vPt[ord]->SetName(histTitle->Data());
455  histFull[k].histFullHar[j].mHist_vPt[ord]->SetXTitle("Pt (GeV)");
456  histFull[k].histFullHar[j].mHist_vPt[ord]->SetYTitle("v (%)");
457  if (ord>0) histFull[k].histFullHar[j].mHist_vPt[ord]->Scale(1./profScale);
458  delete histTitle;
459  delete histTitle2;
460 
461 
462  //do not reverse the order of getting v2D and vPt !!
463 
464  histTitle = new TString("Flow_Cumul_v2D_Order");
465  histTitle->Append(*theCumulOrderChar);
466  histTitle->Append("_Sel");
467  histTitle->Append(*countSels);
468  histTitle->Append("_Har");
469  histTitle->Append(*countHars);
470  histFull[k].histFullHar[j].mHist_v2D[ord] =
471  new TH2D(*(histFull[k].histFullHar[j].mHistCumul2D[ord]->
472  ProjectionXY(histTitle->Data(),"e")));
473  histFull[k].histFullHar[j].mHist_v2D[ord]->SetTitle(histTitle->Data());
474  histFull[k].histFullHar[j].mHist_v2D[ord]->SetXTitle((char*)xLabel.Data());
475  histFull[k].histFullHar[j].mHist_v2D[ord]->SetYTitle("Pt (GeV)");
476  histFull[k].histFullHar[j].mHist_v2D[ord]->SetZTitle("v (%)");
477  if (ord>0) histFull[k].histFullHar[j].mHist_v2D[ord]->Scale(1./profScale);
478  delete histTitle;
479 
480 
481 
482  }
483 
484 
485  //In the case if cumulant fails, let <V> etc. to be 1. to avoid crash.
486  meanIntegV[j] = (cumulInteg1[j]>0.) ? sqrt(cumulInteg1[j]) : 1.; // <v> 2-part, m=1
487  meanIntegV2[j] = (cumulInteg1[j]>0.) ? cumulInteg1[j] : 1.; // <v**2> 2-part, m=2
488  meanIntegV3[j] = (cumulInteg2[j]<0.) ? pow(-1.*cumulInteg2[j], 3./4.) : 1.; // <v**3> 4-part, m=1
489  meanIntegV4[j] = (cumulInteg2[j]<0.) ? -1.*cumulInteg2[j] : 1.; // <v**4> 4-part, m=2
490 
491 
492 
493 
494  if (m_M==1) { // Eq. (PG14)
495 
496  histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
497  histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-1./(meanIntegV3[j]*perCent));// (34b)
498  histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
499  histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
500  histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
501  histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
502  } else if (m_M==2) {
503  histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
504  histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
505  histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
506  histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-0.5/(meanIntegV4[j]*perCent) ); // (35b)
507  histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
508  histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
509  }
510 
511  if (cumulInteg1[j]<0.) {
512  cout<<" Sel"<<k+1<<", <V**2> less than zero ! v"<<j+1<<" from 2 particle correlation failed."<<endl;
513  histFull[k].histFullHar[j].mHist_v2D[0]->Reset();
514  histFull[k].histFullHar[j].mHist_vEta[0]->Reset();
515  histFull[k].histFullHar[j].mHist_vPt[0]->Reset();
516  }
517 
518 
519  if (-1.*cumulInteg2[j]<0.) {
520  cout<<" Sel"<<k+1<<", <V**4> less than zero ! v"<<j+1<<" from 4 particle correlation failed."<<endl;
521  histFull[k].histFullHar[j].mHist_v2D[1]->Reset();
522  histFull[k].histFullHar[j].mHist_vEta[1]->Reset();
523  histFull[k].histFullHar[j].mHist_vPt[1]->Reset();
524  }
525 
526  }// end of for (int j = 0; j < nHars; j++)
527 
528 
529 
530 
532 
533 
534  for (int j = 0; j < nHars; j++) {
535 
536  if (j != 0) continue; //only j==0 makes sense for mixed har.
537 
538  char countHars[2];
539  sprintf(countHars,"%d",j+1);
540 
541  cout<<" ========== Sel"<<k+1<<" Har"<<j+1<<" ===== v1{3} ======"<<endl;
542 
543  double mAvMult =
544  histFull[k].histFullHar[j].mHistMultSum->GetBinContent(1)/
545  histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
546 
547  // cout<<" mAveMult "<<mAvMult<<endl;
548 
549  double mAveMeanWgtSqr = // <w**2>
550  histFull[k].histFullHar[j].mHistMeanWgtSqrSum->GetBinContent(1)/
551  histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
552 
553  cout<<" histFull[k].histFullHar[j].mHistMeanWgtSqrSum->GetBinContent(1) "<<histFull[k].histFullHar[j].mHistMeanWgtSqrSum->GetBinContent(1)<<endl;
554  cout<<" histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1) "<<histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1)<<endl;
555 
556 
557  double CpqMix[nCumulMixHar_pMax][nCumulMixHar_qMax];//(30)
558  for (int pq = 0; pq < nCumulMixHar_pMax*nCumulMixHar_qMax; pq++) {
559  int pIndex = pq/nCumulMixHar_qMax;
560  int qIndex = pq%nCumulMixHar_qMax;
561 
562 
563 
564  histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] =
565  histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq]->GetBinContent(1)/
566  histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1);
567 
568 
569  CpqMix[pIndex][qIndex]=mAvMult*(pow(histFull[k].histFullHar[j].mCumulIntegG0Mix[pq], 1./mAvMult) -1.);
570 
571 
572  }
573 
574  double CpxMix[nCumulMixHar_pMax];
575  double CpyMix[nCumulMixHar_pMax];
576 
577 
578  for (int p =0; p<nCumulMixHar_pMax; p++){
579  CpxMix[p] = (1./(4.*r0Mix))*(CpqMix[p][0] - CpqMix[p][2]);
580  CpyMix[p] = (1./(4.*r0Mix))*(CpqMix[p][3] - CpqMix[p][1]);
581 
582  }
583 
584 
585  cumulIntegMix[j] = (1./(4.*r0Mix*r0Mix))*(
586  CpxMix[0]-CpyMix[1]-CpxMix[2]+CpyMix[3]
587  +CpxMix[4]-CpyMix[5]-CpxMix[6]+CpyMix[7]);
588 
589  //use v2{4}
590  double tempMeanV = pow(-1.*cumulInteg2[1]*mAveMeanWgtSqr*mAveMeanWgtSqr,1./4.); // <wgt*v2>
591  //use v2{2}
592  // double tempMeanV = pow(cumulInteg1[1]*mAveMeanWgtSqr,1./2.); // <wgt*v2>
593 
594  cout<<"mAveMeanWgtSqr "<<mAveMeanWgtSqr<<endl;
595  cout<<" cmumulInteg2[1] "<<cumulInteg2[1]<<endl;
596  cout<<" tempMeanV "<<tempMeanV<<endl;
597  double tempMeanVMixSq = cumulIntegMix[j]/tempMeanV; //(24)
598 
599 
602  // double chi1 = sqrt(tempMeanVMixSq)*sqrt(mAvMult); //(36)
603  // double chi2 = tempMeanV*sqrt(mAvMult);
604  // double deltaC3Sqr = (2. + 4.*chi1*chi1 + 2.*chi2*chi2 + 4.*chi1*chi1*chi2*chi2+chi1*chi1*chi1*chi1)/(2.*mAvMult*mAvMult*mAvMult*(histFull[k].histFullHar[j].mHistNEvent->GetBinContent(1)));
605  //
606  // cout<<" C3 is "<<cumulIntegMix[j]<<endl;
607  // cout<<" delta C3 square is "<<deltaC3Sqr<<endl;
609 
610 
611  histTitle = new TString("Flow_CumulMix_v2D_Sel");
612  histTitle->Append(*countSels);
613  histTitle->Append("_Har");
614  histTitle->Append(*countHars);
615  histFull[k].histFullHar[j].mHistMix_v2D =
616  new TH2D(*(histFull[k].histFullHar[j].mHistCumulMix2D->
617  ProjectionXY(histTitle->Data(),"e")));
618  histFull[k].histFullHar[j].mHistMix_v2D->SetTitle(histTitle->Data());
619  histFull[k].histFullHar[j].mHistMix_v2D->SetXTitle((char*)xLabel.Data());
620  histFull[k].histFullHar[j].mHistMix_v2D->SetYTitle("Pt (GeV)");
621  histFull[k].histFullHar[j].mHistMix_v2D->SetZTitle("v (%)");
622  histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./profScale);
623 
624  delete histTitle;
625 
626 
627 
628  histTitle = new TString("Flow_CumulMix_vEta_Sel");
629  histTitle->Append(*countSels);
630  histTitle->Append("_Har");
631  histTitle->Append(*countHars);
632  histFull[k].histFullHar[j].mHistMix_vEta =
633  new TH1D(*(histFull[k].histFullHar[j].mHistCumulMixEta->
634  ProjectionX(histTitle->Data(),"e")));
635  histFull[k].histFullHar[j].mHistMix_vEta->SetTitle(histTitle->Data());
636  histFull[k].histFullHar[j].mHistMix_vEta->SetXTitle((char*)xLabel.Data());
637  histFull[k].histFullHar[j].mHistMix_vEta->SetYTitle("v (%)");
638  histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./profScale);
639 
640 
641  delete histTitle;
642 
643 
644 
645 
646 
647  histTitle = new TString("Flow_CumulMix_vPt_Sel");
648  histTitle->Append(*countSels);
649  histTitle->Append("_Har");
650  histTitle->Append(*countHars);
651 
652  // histTitle->Append("_TPC");
653 
654  histTitle2 = new TString("Flow_CumulMix2D_Sel");
655  histTitle2->Append(*countSels);
656  histTitle2->Append("_Har");
657  histTitle2->Append(*countHars);
658 
659 
660 
661 
662  histFull[k].histFullHar[j].mHistMix_vPt =
663  //myProjectionY((TProfile2D *)f->Get(histTitle2->Data()) , 0., 1. , 1, ((j==0)?1:0) );
664  new TH1D(*(histFull[k].histFullHar[j].mHistCumulMixPt->ProjectionX(histTitle->Data(),"e")));
665 
666  histFull[k].histFullHar[j].mHistMix_vPt->SetTitle(histTitle->Data());
667  histFull[k].histFullHar[j].mHistMix_vPt->SetXTitle("Pt (GeV)");
668  histFull[k].histFullHar[j].mHistMix_vPt->SetYTitle("v (%)");
669  histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./profScale);
670 
671  delete histTitle;
672  delete histTitle2;
673 
674 
675 
676  if (tempMeanVMixSq>0.){
677  meanIntegVMix[j] = sqrt(tempMeanVMixSq);//(24) <w1v1>
678 
679  histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
680  histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
681  histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
682  histFull[k].mHistMix_v->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
683  } else {
684  histFull[k].histFullHar[j].mHistMix_v2D->Reset();
685  histFull[k].histFullHar[j].mHistMix_vEta->Reset();
686  histFull[k].histFullHar[j].mHistMix_vPt->Reset();
687  histFull[k].mHistMix_v->Reset();
688  cout<<"### <wgt*v1>**2 = "<<tempMeanVMixSq<<" < 0. 3-part mixed har ana. failed "<<endl;
689  }
690 
691 
692 
693  } //end of Har loop for Mixed har.
694 
695 
696 
697 
698  if (m_M==1) {
699 
700  TH1D* histOfMeanIntegV = new TH1D(*(histFull[k].mHist_v[0]));
701  histOfMeanIntegV->Reset();
702 
703  TH1D* histOfMeanIntegV3 = new TH1D(*(histFull[k].mHist_v[1]));
704  histOfMeanIntegV3->Reset();
705 
706  for (int j = 1; j < nHars+1; j++) {
707  histOfMeanIntegV->SetBinContent(j, 1./(meanIntegV[j-1]*perCent));
708  histOfMeanIntegV->SetBinError(j,0.);
709  histOfMeanIntegV3->SetBinContent(j, -1./(meanIntegV3[j-1]*perCent));
710  histOfMeanIntegV3->SetBinError(j,0.);
711  }
712  histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV);
713  histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV3);
714 
715  //force to be zero if the value is "nan"
716  for (int ord = 0; ord < nCumulDiffOrders; ord++){
717  for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
718  if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
719  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ) {
720  histFull[k].mHist_v[ord]->SetBinContent(j,0.);
721  histFull[k].mHist_v[ord]->SetBinError(j,0.);
722  }
723  }
724  }
725 
726  for (int j = 1; j < nHars+1; j++) {
727  if (histFull[k].mHist_v[0]->GetBinContent(j)< FLT_MAX &&
728  histFull[k].mHist_v[0]->GetBinContent(j)> -1.*FLT_MAX)
729  cout << "##### 2-part v" << j << " = ("
730  << histFull[k].mHist_v[0]->GetBinContent(j)
731  <<" +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<" )"<<endl;
732  if (histFull[k].mHist_v[1]->GetBinContent(j)< FLT_MAX &&
733  histFull[k].mHist_v[1]->GetBinContent(j)> -1.*FLT_MAX)
734  cout << "##### 4-part v" << j << " = ("
735  << histFull[k].mHist_v[1]->GetBinContent(j)
736  <<" +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<" )"<<endl;
737  }
738 
739  delete histOfMeanIntegV;
740  delete histOfMeanIntegV3;
741 
742  } else if (m_M==2) {
743 
744  TH1D* histOfMeanIntegV2 = new TH1D(*(histFull[k].mHist_v[0]));
745  histOfMeanIntegV2->Reset();
746 
747  TH1D* histOfMeanIntegV4 = new TH1D(*(histFull[k].mHist_v[1]));
748  histOfMeanIntegV4->Reset();
749 
750  for (int j = 1; j < nHars+1; j++) {
751  histOfMeanIntegV2->SetBinContent(j, 1./(meanIntegV2[j-1]*perCent));
752  histOfMeanIntegV2->SetBinError(j,0.);
753  histOfMeanIntegV4->SetBinContent(j, -0.5/(meanIntegV4[j-1]*perCent));
754  histOfMeanIntegV4->SetBinError(j,0.);
755  }
756  histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV2);
757  histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV4);
758 
759  //force to be zero if the value is "nan"
760  for (int ord = 0; ord < nCumulDiffOrders; ord++){
761  for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
762  if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
763  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ){
764  histFull[k].mHist_v[ord]->SetBinContent(j,0.);
765  histFull[k].mHist_v[ord]->SetBinError(j,0.);
766  }
767  }
768  }
769 
770  for (int j = 1; j < nHars+1; j++) {
771 
772  if (histFull[k].mHist_v[0]->GetBinContent(j)< FLT_MAX &&
773  histFull[k].mHist_v[0]->GetBinContent(j)> -1.*FLT_MAX)
774  cout << "##### 2-part v" << j << " = ("
775  << histFull[k].mHist_v[0]->GetBinContent(j)
776  <<") +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<endl;
777  if (histFull[k].mHist_v[1]->GetBinContent(j)< FLT_MAX &&
778  histFull[k].mHist_v[1]->GetBinContent(j)> -1.*FLT_MAX)
779  cout << "##### 4-part v" << j << " = ("
780  << histFull[k].mHist_v[1]->GetBinContent(j)
781  <<") +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<endl;
782  }
783 
784  delete histOfMeanIntegV2;
785  delete histOfMeanIntegV4;
786  }
787 
788  }
789 
790 
791 
792  // ============== write out ====================
793 
794  f->cd();
795 
796  for (int k = 0; k < nSels; k++) {
797 
798  for (int ord = 0; ord < nCumulDiffOrders; ord++) {
799  histFull[k].mHist_v[ord]->Write(histFull[k].mHist_v[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
800  }
801 
802  histFull[k].mHistMix_v->Write(histFull[k].mHistMix_v->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
803 
804 
805  for (int j = 0; j < nHars; j++) {
806  cout<<" writting ........................."<<endl;
807  for (int ord = 0; ord < nCumulDiffOrders; ord++){
808  histFull[k].histFullHar[j].mHist_v2D[ord]->Write(histFull[k].histFullHar[j].mHist_v2D[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
809  histFull[k].histFullHar[j].mHist_vEta[ord]->Write(histFull[k].histFullHar[j].mHist_vEta[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
810  histFull[k].histFullHar[j].mHist_vPt[ord]->Write(histFull[k].histFullHar[j].mHist_vPt[ord]->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
811  }
812 
813 
814  if (j==0){
815  cout<<" j "<<j<<" k "<<k<<endl;
816  cout<<" histFull[k].histFullHar[j].mHistMix_v2D "<<histFull[k].histFullHar[j].mHistMix_v2D->GetTitle()<<endl;
817 
818  (histFull[k].histFullHar[j].mHistMix_v2D)->Write((histFull[k].histFullHar[j].mHistMix_v2D)->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
819  (histFull[k].histFullHar[j].mHistMix_vEta)->Write((histFull[k].histFullHar[j].mHistMix_vEta)->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
820  (histFull[k].histFullHar[j].mHistMix_vPt)->Write((histFull[k].histFullHar[j].mHistMix_vPt)->GetTitle(),TObject::kOverwrite | TObject::kSingleKey);
821 
822  }
823 
824  }
825  }
826  //=============================================
827 
828  f->Close();
829 
830 }
Generated by   doxygen 1.8.5