StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFlowCumulantMaker.cxx
1 //
3 // $Id: StFlowCumulantMaker.cxx,v 1.22 2006/02/22 19:12:29 posk Exp $
4 //
5 // Authors: Aihong Tang, Kent State U. Oct 2001
6 // Frame adopted from Art and Raimond's StFlowAnalysisMaker.
7 //
9 //
10 // Description: Maker to analyze Flow using the latest cumulant method.
11 // refer to Phy. Rev. C63 (2001) 054906 (old new method)
12 // and Phy. Rev. C64 (2001) 054901 (new new method)
13 // and nucl-ex/0110016 (Practical Guide)
14 // and PRy. Rev. C66 (2002) 014905 (directed flow from elliptic flow v1{3})
15 // Eq. numbers are from the PG or new new if not specified.
16 //
17 // Anything that has "Mix" in it is for v1{3} calculation
18 //
20 
21 #include <Stiostream.h>
22 #include <stdlib.h>
23 #include <math.h>
24 #include <float.h>
25 #include "StMaker.h"
26 #include "StFlowCumulantMaker.h"
27 #include "StFlowMaker/StFlowMaker.h"
28 #include "StFlowMaker/StFlowEvent.h"
29 #include "StFlowMaker/StFlowConstants.h"
30 #include "StFlowMaker/StFlowSelection.h"
31 #include "StEnumerations.h"
32 #include "PhysicalConstants.h"
33 #include "SystemOfUnits.h"
34 #include "TFile.h"
35 #include "TString.h"
36 #include "TObjString.h"
37 #include "TVectorD.h"
38 #include "TH1.h"
39 #include "TH2.h"
40 #include "TProfile.h"
41 #include "TProfile2D.h"
42 #include "TOrdCollection.h"
43 #include "StMessMgr.h"
44 #include "TMath.h"
45 
46 #define PR(x) cout << "##### FlowCumulantAnalysis: " << (#x) << " = " << (x) << endl;
47 
48 ClassImp(StFlowCumulantMaker)
49 
50 //-----------------------------------------------------------------------
51 
52 StFlowCumulantMaker::StFlowCumulantMaker(const Char_t* name): StMaker(name),
53  MakerName(name) {
54  pFlowSelect = new StFlowSelection();
55  SetHistoRanges();
56 }
57 
58 StFlowCumulantMaker::StFlowCumulantMaker(const Char_t* name,
59  const StFlowSelection& flowSelect) :
60  StMaker(name), MakerName(name) {
61  pFlowSelect = new StFlowSelection(flowSelect); // copy constructor
62  SetHistoRanges();
63 }
64 
65 //-----------------------------------------------------------------------
66 
67 StFlowCumulantMaker::~StFlowCumulantMaker() {
68 }
69 
70 //-----------------------------------------------------------------------
71 
72 Int_t StFlowCumulantMaker::Make() {
73  // Fill histograms
74 
75  // Get a pointer to StFlowEvent
76  StFlowMaker* pFlowMaker = NULL;
77  pFlowMaker = (StFlowMaker*)GetMaker("Flow");
78  if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
79  if (pFlowEvent && pFlowSelect->Select(pFlowEvent)) { // event selected
80  FillFromFlowEvent(); // get event quantities
81  FillEventHistograms(); // fill from FlowEvent
82  if (pFlowEvent) FillParticleHistograms(); // fill particle histograms
83  if (Debug()) StMaker::PrintInfo();
84  } else {
85  gMessMgr->Info("##### FlowCumulantMaker: FlowEvent pointer null");
86  }
87 
88  return kStOK;
89 }
90 
91 //-----------------------------------------------------------------------
92 
93 Int_t StFlowCumulantMaker::Init() {
94  // Book histograms
95 
96  float ptMaxPart = Flow::ptMaxPart;
97  if (pFlowSelect->PtMaxPart()) {
98  ptMaxPart = pFlowSelect->PtMaxPart();
99  }
100 
101  nPtBinsPart = Flow::nPtBinsPart;
102  if (pFlowSelect->PtBinsPart()) {
103  nPtBinsPart = pFlowSelect->PtBinsPart();
104  }
105 
106  xLabel = "Pseudorapidity";
107  if (strlen(pFlowSelect->PidPart()) != 0) { xLabel = "Rapidity"; }
108 
109  profScale=1000.;// in new root version, error does not show correctly
110  //if the number is too small. this is a trick to walkaround.
111 
112  r0 = 1.5; // this number should be small, but it could bring numerical
113  // error if it is too small.
114  r0Sq = r0 * r0;
115  r0Mix= 3.; //r0 for v1{3} calculation
116 
117  m_M = 1; // if m_M = 2, what measured is v2{v1v1v2}, v4{v2v2v4}, v6{v3v3v6}... etc.
118 
119  bool noDenomFileWarned = kFALSE;
120  TString* histTitle;
121 
122  for (int k = 0; k < Flow::nSels; k++) {
123 
124  // for each selection
125  histFull[k].mHistCumul = new TProfile*[Flow::nCumulDiffOrders];
126 
127  // mixed cumulant. (1st har mix with 2nd har. for v1 analysis).
128  //we study 3-part v1 only, so we do not need sth. like nCumulDiffOrders as above.
129  histTitle = new TString("Flow_CumulMix");
130  histTitle->Append("_Sel");
131  *histTitle +=k+1;
132  histFull[k].mHistCumulMix =
133  new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5,
134  (float)(Flow::nHars) + 0.5, -1.*FLT_MAX, FLT_MAX, "");
135  histFull[k].mHistCumulMix->SetXTitle("place for saving mixed cumulant");
136  delete histTitle;
137 
138 
139  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
140  char theCumulOrder[2]; // if >10, need to use char*
141  sprintf(theCumulOrder,"%d",(ord+1)*2);
142  histTitle = new TString("Flow_Cumul_Order");
143  *histTitle->Append(*theCumulOrder);
144  histTitle->Append("_Sel");
145  *histTitle += k+1;
146  histFull[k].mHistCumul[ord] =
147  new TProfile(histTitle->Data(), histTitle->Data(), Flow::nHars, 0.5,
148  (float)(Flow::nHars) + 0.5, -1.*FLT_MAX, FLT_MAX, "");
149  histFull[k].mHistCumul[ord]->SetXTitle("harmonic");
150  delete histTitle;
151  }
152 
153  // for each harmonic
154  for (int j = 0; j < Flow::nHars; j++) {
155 
156  histTitle = new TString("Flow_CumulMultSum_Sel");
157  *histTitle += k+1;
158  histTitle->Append("_Har");
159  *histTitle += j+1;
160  histFull[k].histFullHar[j].mHistMultSum =
161  new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
162  delete histTitle;
163 
164  histTitle = new TString("Flow_CumulMeanWgtSqrSum_Sel");
165  *histTitle +=k+1;
166  histTitle->Append("_Har");
167  *histTitle +=j+1;
168  histFull[k].histFullHar[j].mHistMeanWgtSqrSum =
169  new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
170  delete histTitle;
171 
172 
173  histTitle = new TString("Flow_CumulNEvent_Sel");
174  *histTitle += k+1;
175  histTitle->Append("_Har");
176  *histTitle += j+1;
177  histFull[k].histFullHar[j].mHistNEvent =
178  new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
179  delete histTitle;
180 
181  // **** for differential flow ****
182 
183  // cumulant d_p/n{k}
184  // p is the harmonic of the differential flow
185  // n is the harmonic of the integrated flow used for the differential flow
186  // measurment. p = mn = either n or 2n.
187 
188  // if m_M=1, k=2 dn/n{2} : cumulant from 2-part corr.
189  // if m_M=1, k=4 dn/n{4} : cumulant from 4-part corr.
190  // if m_M=2, k=2 d2n/n{3} : cumulant from 3-part corr., mixed harmonic
191  // if m_M=2, k=4 d2n/n{5} : cumulant from 5-part corr., mixed harmonic
192  // where {2},{3} corresponds to theCumulIndex=1 below.
193  // {4},{5} corresponds to theCumulIndex=2 below.
194 
195  histFull[k].histFullHar[j].mHistCumul2D =
196  new TProfile2D*[Flow::nCumulDiffOrders];
197  histFull[k].histFullHar[j].mHistCumulEta =
198  new TProfile*[Flow::nCumulDiffOrders];
199  histFull[k].histFullHar[j].mHistCumulPt =
200  new TProfile*[Flow::nCumulDiffOrders];
201 
202  // For each cumulant order
203  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
204  char theCumulOrder[2]; // if >10, need to use char*
205  sprintf(theCumulOrder,"%d",(ord+1)*2);
206 
207  histTitle = new TString("Flow_Cumul2D_Order");
208  histTitle->Append(*theCumulOrder);
209  histTitle->Append("_Sel");
210  *histTitle += k+1;
211  histTitle->Append("_Har");
212  *histTitle += j+1;
213  histFull[k].histFullHar[j].mHistCumul2D[ord] =
214  new TProfile2D(histTitle->Data(),histTitle->Data(), mNEtaBins,
215  mEtaMin, mEtaMax, nPtBinsPart, Flow::ptMin,
216  ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
217  histFull[k].histFullHar[j].mHistCumul2D[ord]->SetXTitle((char*)xLabel.Data());
218  histFull[k].histFullHar[j].mHistCumul2D[ord]->SetYTitle("Pt (GeV/c)");
219  delete histTitle;
220 
221  histTitle = new TString("Flow_CumulEta_Order");
222  histTitle->Append(*theCumulOrder);
223  histTitle->Append("_Sel");
224  *histTitle += k+1;
225  histTitle->Append("_Har");
226  *histTitle += j+1;
227  histFull[k].histFullHar[j].mHistCumulEta[ord] =
228  new TProfile(histTitle->Data(),histTitle->Data(), mNEtaBins,
229  mEtaMin, mEtaMax, -1.*FLT_MAX, FLT_MAX, "");
230  histFull[k].histFullHar[j].mHistCumulEta[ord]->SetXTitle((char*)xLabel.Data());
231  delete histTitle;
232 
233  histTitle = new TString("Flow_CumulPt_Order");
234  histTitle->Append(*theCumulOrder);
235  histTitle->Append("_Sel");
236  *histTitle += k+1;
237  histTitle->Append("_Har");
238  *histTitle += j+1;
239  histFull[k].histFullHar[j].mHistCumulPt[ord] =
240  new TProfile(histTitle->Data(), histTitle->Data(), nPtBinsPart,
241  Flow::ptMin, ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
242  histFull[k].histFullHar[j].mHistCumulPt[ord]->SetXTitle("Pt (GeV/c)");
243  delete histTitle;
244 
245  }
246 
247 
248  // for mixing 1st and 2nd harmonic. only coutHars=1 is meaningful.
249  histTitle = new TString("Flow_CumulMix2D");
250  histTitle->Append("_Sel");
251  *histTitle +=k+1;
252  histTitle->Append("_Har");
253  *histTitle +=j+1;
254  histFull[k].histFullHar[j].mHistCumulMix2D =
255  new TProfile2D(histTitle->Data(),histTitle->Data(), Flow::nEtaBins,
256  Flow::etaMin, Flow::etaMax, nPtBinsPart, Flow::ptMin,
257  ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
258  histFull[k].histFullHar[j].mHistCumulMix2D->SetXTitle((char*)xLabel.Data());
259  histFull[k].histFullHar[j].mHistCumulMix2D->SetYTitle("Pt (GeV)");
260  delete histTitle;
261 
262 
263  histTitle = new TString("Flow_CumulMixEta");
264  histTitle->Append("_Sel");
265  *histTitle +=k+1;
266  histTitle->Append("_Har");
267  *histTitle +=j+1;
268  histFull[k].histFullHar[j].mHistCumulMixEta =
269  new TProfile(histTitle->Data(),histTitle->Data(), Flow::nEtaBins,
270  Flow::etaMin, Flow::etaMax, -1.*FLT_MAX, FLT_MAX, "");
271  histFull[k].histFullHar[j].mHistCumulMixEta->SetXTitle((char*)xLabel.Data());
272  delete histTitle;
273 
274 
275  histTitle = new TString("Flow_CumulMixPt");
276  histTitle->Append("_Sel");
277  *histTitle +=k+1;
278  histTitle->Append("_Har");
279  *histTitle +=j+1;
280  histFull[k].histFullHar[j].mHistCumulMixPt =
281  new TProfile(histTitle->Data(), histTitle->Data(), nPtBinsPart,
282  Flow::ptMin, ptMaxPart, -1.*FLT_MAX, FLT_MAX, "");
283  histFull[k].histFullHar[j].mHistCumulMixPt->SetXTitle("Pt (GeV)");
284  delete histTitle;
285 
286 
287  histFull[k].histFullHar[j].mCumulG0Denom =
288  new TProfile*[Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax];
289 
290  // For each cumulant order * qMax orders
291  for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
292 
293  TString* histTitleIntegDenom; //for read in
294 
295  int cumulIndex = (pq/Flow::nCumulDiff_qMax) + 1; // like 1,2,3.
296  // for cumulant order 2,4,6 not begining with 0. This is "p" in Eq. (B1, PG5).
297  int qIndex = pq%Flow::nCumulDiff_qMax; // like 0,1,..._qMax-1
298  // begin with 0. This is "q" in Eq. (B3, PG5).
299 
300  char theCumulOrderChar[2];
301  char qIndexOrderChar[2]; // if >10, need to use char*
302  sprintf(theCumulOrderChar,"%d",cumulIndex*2);
303  sprintf(qIndexOrderChar,"%d",qIndex);
304 
305  histTitle = new TString("Flow_CumulDenom_Order");
306  histTitle->Append(*theCumulOrderChar);
307  histTitle->Append("_GenFcnqIdx");
308  histTitle->Append(*qIndexOrderChar);
309  histTitle->Append("_Sel");
310  *histTitle += k+1;
311  histTitle->Append("_Har");
312  *histTitle += j+1;
313  histFull[k].histFullHar[j].mCumulG0Denom[pq] =
314  new TProfile(histTitle->Data(),histTitle->Data(), 1,0., 1., -1.*FLT_MAX, FLT_MAX, "");
315  histFull[k].histFullHar[j].mCumulG0Denom[pq]->SetYTitle("<G>");
316  histTitleIntegDenom = new TString(histTitle->Data());
317  delete histTitle;
318 
319  // Open the file and get the histograms
320  // do not modify this section [ROOT bug]
321  TFile f("denominator.root","R");
322  if (f.IsOpen()) {
323 
324  f.cd();
325  TProfile* tempDenomProfile =
326  dynamic_cast<TProfile*>(f.Get(histTitleIntegDenom->Data()));
327  if (!tempDenomProfile) {
328  cout << "##### FlowCumulantAnalysis: can not find " <<
329  histTitleIntegDenom->Data() << endl;
330  return kFALSE;
331  }
332  delete histTitleIntegDenom;
333 
334  histFull[k].histFullHar[j].mCumulG0DenomRead[pq]
335  = tempDenomProfile->GetBinContent(1);
336 
337  f.Close();
338 
339  } else {
340 
341  if(!noDenomFileWarned) {
342  gMessMgr->Info("##### FlowCumulantAnalysis:denominator.root is not present, assumming this run is just for producing denominator.root. That means cumulant flow result in flow.cumulant.root is nonsense for this run. ");
343  noDenomFileWarned = kTRUE;
344  }
345  histFull[k].histFullHar[j].mCumulG0DenomRead[pq] = 1.;
346 
347  }
348 
349  double theTempPhi = twopi*((double)qIndex) / // Eq. (PG5)
350  ((double)Flow::nCumulDiff_qMax);
351  double theRz = r0*::sqrt((double)cumulIndex);
352  histFull[k].histFullHar[j].mDiffXz[pq] = theRz*cos(theTempPhi);
353  histFull[k].histFullHar[j].mDiffYz[pq] = theRz*sin(theTempPhi);
354  }
355 
356 
357  histFull[k].histFullHar[j].mCumulG0MixDenom =
358  new TProfile*[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
359 
360  histFull[k].histFullHar[j].mHistCumulIntegG0MixSum =
361  new TH1D*[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
362 
363 
364  for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) {
365 
366 
367  TString* histTitleIntegDenomMix; //for read in
368 
369  int pIndex = (pq/Flow::nCumulMixHar_qMax);
370  int qIndex = pq%Flow::nCumulMixHar_qMax;
371 
372  char pIndexChar[2];
373  char qIndexChar[2]; // if >10, need to use char*
374  sprintf(pIndexChar,"%d",pIndex);
375  sprintf(qIndexChar,"%d",qIndex);
376 
377 
378  histTitle = new TString("Flow_CumulMixDenom");
379  histTitle->Append("_GenFunIdxp");
380  histTitle->Append(*pIndexChar);
381  histTitle->Append("_GenFunIdxq");
382  histTitle->Append(*qIndexChar);
383  histTitle->Append("_Sel");
384  *histTitle +=k+1;
385  histTitle->Append("_Har");
386  *histTitle +=j+1;
387  histFull[k].histFullHar[j].mCumulG0MixDenom[pq] =
388  new TProfile(histTitle->Data(),histTitle->Data(), 1,0., 1., -1.*FLT_MAX, FLT_MAX, "");
389  histFull[k].histFullHar[j].mCumulG0MixDenom[pq]->SetYTitle("<G>");
390  histTitleIntegDenomMix = new TString(histTitle->Data());
391  delete histTitle;
392 
393 
394  histTitle = new TString("Flow_CumulIntegG0MixSum");
395  histTitle->Append("_GenFunIdxp");
396  histTitle->Append(*pIndexChar);
397  histTitle->Append("_GenFunIdxq");
398  histTitle->Append(*qIndexChar);
399  histTitle->Append("_Sel");
400  *histTitle +=k+1;
401  histTitle->Append("_Har");
402  *histTitle +=j+1;
403  histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq] =
404  new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
405  delete histTitle;
406 
407 
408  histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] = 0.;
409 
410 
411  TFile f("denominator.root","R");
412  if (f.IsOpen()) {
413 
414  f.cd();
415  TProfile* tempDenomMixProfile =
416  dynamic_cast<TProfile*>(f.Get(histTitleIntegDenomMix->Data()));
417  if (!tempDenomMixProfile) {
418  cout << "##### FlowCumulantAnalysis: can not find " <<
419  histTitleIntegDenomMix->Data() << endl;
420  return kFALSE;
421  }
422  delete histTitleIntegDenomMix;
423 
424  histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq]
425  = tempDenomMixProfile->GetBinContent(1);
426 
427  f.Close();
428 
429  } else {
430 
431  if(!noDenomFileWarned) {
432  gMessMgr->Info("##### FlowCumulantAnalysis:denominator.root is not present, assumming this run is just for producing denominator.root. That means cumulant flow result with mixed harmonics in flow.cumulant.root is nonsense for this run. ");
433  noDenomFileWarned = kTRUE;
434  }
435  histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq] = 1.;
436 
437  }
438 
439  }
440 
441 
442  for (int p = 0; p < Flow::nCumulMixHar_pMax; p++){ //v1{3} paper Eq (29)
443  histFull[k].histFullHar[j].mMixX1z[p]=r0Mix*cos( pi*((double)p)/4. );
444  histFull[k].histFullHar[j].mMixY1z[p]=r0Mix*sin( pi*((double)p)/4. );
445  }
446 
447  for (int q = 0; q < Flow::nCumulMixHar_qMax; q++){
448  histFull[k].histFullHar[j].mMixX2z[q]=r0Mix*cos( pi*((double)q)/2. );
449  histFull[k].histFullHar[j].mMixY2z[q]=r0Mix*sin( pi*((double)q)/2. );
450  }
451 
452 
453  // *** for integrated flow ***
454 
455  histFull[k].histFullHar[j].mHistCumulIntegG0Sum =
456  new TH1D*[Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax];
457 
458  for (int pq = 0; pq < Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax; pq++) {
459  int cumulIndex = (pq/Flow::nCumulInteg_qMax) + 1; // like 1,2,3.
460  //not begining with 0. That is "p" in Eq. (B1, PG5).
461  int qIndex = pq%(Flow::nCumulInteg_qMax); // like 0,1,..qMax-1
462  //begining with 0. Just like Eq. (B3, PG5).
463 
464  char theCumulOrderChar[2];
465  char qIndexOrderChar[2]; // if >10, need to use char*
466  sprintf(theCumulOrderChar,"%d",cumulIndex*2);
467  sprintf(qIndexOrderChar,"%d",qIndex);
468 
469  double theTempPhi = twopi*((double)qIndex) /
470  ((double)Flow::nCumulInteg_qMax);
471  double theRz = r0*::sqrt((double)cumulIndex);
472  histFull[k].histFullHar[j].mIntegXz[pq] = theRz*cos(theTempPhi);
473  histFull[k].histFullHar[j].mIntegYz[pq] = theRz*sin(theTempPhi);
474 
475  histFull[k].histFullHar[j].mCumulIntegG0[pq] = 0.;
476 
477  histTitle = new TString("Flow_CumulIntegG0Sum_Order");
478  histTitle->Append(*theCumulOrderChar);
479  histTitle->Append("_GenFunIdx");
480  histTitle->Append(*qIndexOrderChar);
481  histTitle->Append("_Sel");
482  *histTitle += k+1;
483  histTitle->Append("_Har");
484  *histTitle += j+1;
485  histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq] =
486  new TH1D(histTitle->Data(),histTitle->Data(),1,0.,1.);
487  delete histTitle;
488 
489  }
490 
491  histFull[k].histFullHar[j].mMultSum = 0.;
492  histFull[k].histFullHar[j].mNEvent = 0;
493  histFull[k].histFullHar[j].mMeanWgtSqrSum = 0.;
494 
495  }
496  }
497 
498  gMessMgr->SetLimit("##### FlowCumulantAnalysis", 2);
499  gMessMgr->Info("##### FlowCumulantAnalysis: $Id: StFlowCumulantMaker.cxx,v 1.22 2006/02/22 19:12:29 posk Exp $");
500 
501  return StMaker::Init();
502 }
503 
504 //-----------------------------------------------------------------------
505 
506 void StFlowCumulantMaker::FillFromFlowEvent() {
507  // Get event quantities from StFlowEvent
508 
509  for (int k = 0; k < Flow::nSels; k++) {
510  pFlowSelect->SetSelection(k);
511  for (int j = 0; j < Flow::nHars; j++) {
512  pFlowSelect->SetHarmonic(j);
513  pFlowSelect->SetSubevent(-1);
514 
515  // full event quantities
516  mMult[k][j] = pFlowEvent->Mult(pFlowSelect);
517  mSqrtOfSumWgtSqr[k][j] = ::sqrt(pFlowEvent->SumWeightSquare(pFlowSelect));
518 
519  }
520  }
521 }
522 
523 //-----------------------------------------------------------------------
524 
525 void StFlowCumulantMaker::FillEventHistograms() {
526  // Fill histograms with event quantities
527 
528  for (int k = 0; k < Flow::nSels; k++) {
529  pFlowSelect->SetSelection(k);
530  for (int j = 0; j < Flow::nHars; j++) {
531  pFlowSelect->SetHarmonic(j);
532 
533  histFull[k].histFullHar[j].mMultSum += (float)mMult[k][j];
534  histFull[k].histFullHar[j].mNEvent++;
535  histFull[k].histFullHar[j].mMeanWgtSqrSum +=
536  (mSqrtOfSumWgtSqr[k][j]*mSqrtOfSumWgtSqr[k][j])/(float)mMult[k][j];
537 
538 
539  for (int pq = 0; pq < Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax; pq++) {
540  histFull[k].histFullHar[j].mCumulIntegG0[pq] +=
541  pFlowEvent->G_New( pFlowSelect,
542  histFull[k].histFullHar[j].mIntegXz[pq],
543  histFull[k].histFullHar[j].mIntegYz[pq] );
544  }
545 
546  for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) {
547  int pIndex = pq/Flow::nCumulMixHar_qMax;
548  int qIndex = pq%Flow::nCumulMixHar_qMax;
549  if (j ==0) //only for v1
550  histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] +=
551  pFlowEvent->G_Mix( pFlowSelect,
552  histFull[k].histFullHar[j].mMixX1z[pIndex],
553  histFull[k].histFullHar[j].mMixY1z[pIndex],
554  histFull[k].histFullHar[j].mMixX2z[qIndex],
555  histFull[k].histFullHar[j].mMixY2z[qIndex] );
556  }
557 
558  }
559  }
560 
561 }
562 
563 //-----------------------------------------------------------------------
564 
565 void StFlowCumulantMaker::FillParticleHistograms() {
566  // Fill histograms from the particles
567 
568  // Initialize Iterator
569  StFlowTrackCollection* pFlowTracks = pFlowEvent->TrackCollection();
570  StFlowTrackIterator itr;
571 
572  // In the view of good coding, the following block should be placed inside
573  // the track loop. It is put here to reduce run time.
574  // If somebody changes the pFlowSelect to select tracks for v-part,
575  // he/she has to change pFlowSelect in this block for consistency.
576  // Selections have to be consistent between here and the loop.
577 
578  double* evtG[Flow::nSels][Flow::nHars];
579  double* theCrossterm[Flow::nSels][Flow::nHars];
580 
581  double* evtGMix[Flow::nSels][Flow::nHars];
582  double* theCrosstermMix[Flow::nSels][Flow::nHars];
583 
584  for (int k = 0; k < Flow::nSels; k++) {
585  for (int j = 0; j < Flow::nHars; j++) {
586  evtG[k][j] =
587  new double[Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax];
588  theCrossterm[k][j] =
589  new double[Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax];
590  evtGMix[k][j] =
591  new double[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
592  theCrosstermMix[k][j] =
593  new double[Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax];
594  }
595  }
596 
597  for (int k = 0; k < Flow::nSels; k++) {
598  for (int j = 0; j < Flow::nHars; j++) {
599  pFlowSelect->SetSelection(k);
600  pFlowSelect->SetHarmonic(j);
601 
602 
603  for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
604 
605  evtG[k][j][pq] =
606  (pFlowEvent->G_New( pFlowSelect,
607  histFull[k].histFullHar[j].mDiffXz[pq],
608  histFull[k].histFullHar[j].mDiffYz[pq] )) ;
609  theCrossterm[k][j][pq] =evtG[k][j][pq];
610  histFull[k].histFullHar[j].mCumulG0Denom[pq]->Fill(0.5,evtG[k][j][pq]);
611  }
612 
613  for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) {
614  int pIndex = pq/Flow::nCumulMixHar_qMax;
615  int qIndex = pq%Flow::nCumulMixHar_qMax;
616  if (j ==0) //only for v1
617  evtGMix[k][j][pq] =
618  pFlowEvent->G_Mix( pFlowSelect,
619  histFull[k].histFullHar[j].mMixX1z[pIndex],
620  histFull[k].histFullHar[j].mMixY1z[pIndex],
621  histFull[k].histFullHar[j].mMixX2z[qIndex],
622  histFull[k].histFullHar[j].mMixY2z[qIndex] );
623  theCrosstermMix[k][j][pq] =evtGMix[k][j][pq];
624  histFull[k].histFullHar[j].mCumulG0MixDenom[pq]->Fill(0.5,evtGMix[k][j][pq]);
625  }
626 
627  }
628  }
629 
630  for (itr = pFlowTracks->begin(); itr != pFlowTracks->end(); itr++) {
631  StFlowTrack* pFlowTrack = *itr;
632 
633  float phi = pFlowTrack->Phi();
634  if (phi < 0.) phi += twopi;
635  float eta = pFlowTrack->Eta();
636  float pt = pFlowTrack->Pt();
637 
638  for (int k = 0; k < Flow::nSels; k++) {
639  pFlowSelect->SetSelection(k);
640  double cumuTemp[Flow::nCumulDiffOrders];
641  double cumuTempFlip[Flow::nCumulDiffOrders];
642  double cumuTempMix;
643  double cumuTempMixFlip;
644  double order = 0.;
645  double phiWgt = 1.;
646  double phiWgtRaw = 1.; //no need raw phi weight for cumulants.
647 
648  for (int j = 0; j < Flow::nHars; j++) {
649  bool oddHar = (j+1) % 2;
650  pFlowSelect->SetHarmonic(j);
651  order = (double)(j+1);
652 
653  if (pFlowSelect->Select(pFlowTrack)) {
654  // Get phiWgt
655  phiWgt = pFlowEvent->Weight(k, j, pFlowTrack);
656  }
657 
658  // Caculate v for all particles selected for correlation analysis
659  if (pFlowSelect->SelectPart(pFlowTrack)) {
660 
661  float yOrEta =
662  (strlen(pFlowSelect->PidPart()) != 0) ? pFlowTrack->Y() : eta;
663 
664  double Dp[Flow::nCumulDiffOrders]; // the Dp in (B7, PG11)
665  for (int pq = 0; pq < Flow::nCumulDiffOrders; pq++) {
666  Dp[pq] = 0.; }
667 
668  for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
669  int theCumulOrder = (pq/Flow::nCumulDiff_qMax) + 1; // like 1,2.
670  // not begin with 0. which is "p" in (B6, PG5)
671  int qIndex = pq%(Flow::nCumulDiff_qMax); // like 0,1,..qMax-1
672  // begin with 0. which is "q" in (B6, PG5)
673 
674  //Calculate PG11
675  double theCoeff = ::pow(r0*::sqrt(double(theCumulOrder)), double(m_M)) /
676  float(Flow::nCumulDiff_qMax); // first term in (B7)
677  double theCosTerm = cos(twopi*float(qIndex)*float(m_M) /
678  float(Flow::nCumulDiff_qMax)); // cos() in (B7)
679  double theSinTerm = sin(twopi*float(qIndex)*float(m_M) /
680  float(Flow::nCumulDiff_qMax)); // sin() in (B7)
681 
682  if ( (pFlowSelect->SelectPart(pFlowTrack)) &&
683  (pFlowSelect->Select(pFlowTrack)) ) { // remove autocorrelation
684 
685  theCrossterm[k][j][pq] = evtG[k][j][pq] /
686  (1. + (phiWgt/mMult[k][j]) *
687  (2.*histFull[k].histFullHar[j].mDiffXz[pq]*cos(phi*order) +
688  2.*histFull[k].histFullHar[j].mDiffYz[pq]*sin(phi * order) ) );
689  // the argument in the last paragraph in page 9.
690 
691  }
692 
693 
694  double theXpq = (theCrossterm[k][j][pq]*cos(float(m_M) * order * phi)) / // (B6)
695  histFull[k].histFullHar[j].mCumulG0DenomRead[pq];
696  double theYpq = (theCrossterm[k][j][pq]*sin(float(m_M) * order * phi)) / // (B6)
697  histFull[k].histFullHar[j].mCumulG0DenomRead[pq];
698 
699  Dp[theCumulOrder-1] +=
700  theCoeff*(theCosTerm*theXpq + theSinTerm*theYpq); // (B7, PG11)
701 
702  }
703 
704 
706 
707  if (j==0){ //Har 1 only
708 
709  double DpxMix[Flow::nCumulMixHar_pMax];
710  double DpyMix[Flow::nCumulMixHar_pMax];
711 
712  double DpqxMix[Flow::nCumulMixHar_pMax][Flow::nCumulMixHar_qMax];
713  double DpqyMix[Flow::nCumulMixHar_pMax][Flow::nCumulMixHar_qMax];
714 
715 
716  for (int p = 0; p < Flow::nCumulMixHar_pMax; p++) {
717  DpxMix[p] = 0.; DpyMix[p] = 0.;}
718 
719  for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax; pq++) {
720 
721  int pIndex = pq/Flow::nCumulMixHar_qMax;
722  int qIndex = pq%Flow::nCumulMixHar_qMax;
723 
724  if ( (pFlowSelect->SelectPart(pFlowTrack)) &&
725  (pFlowSelect->Select(pFlowTrack)) ) { // rm autocorrelation
726 
727  double etaWgt = (oddHar) ? ( (eta>0) ? (pFlowEvent->EtaAbsWgtValue(eta)) : (-1.*(pFlowEvent->EtaAbsWgtValue(eta))) ) : 1.;
728 
729  double ptWgt = pFlowEvent->PtAbsWgtValue(pt);
730 
731  double detectorV1Wgt = 1.;
732 
733  if (pFlowTrack->TopologyMap().hasHitInDetector(kTpcId) ||
734  (pFlowTrack->TopologyMap().data(0) == 0 &&
735  pFlowTrack->TopologyMap().data(1) == 0)) {
736  // hits in Tpc or TopologyMap not available
737  detectorV1Wgt =pFlowEvent->V1TPCDetctWgtG_Mix(k);
738  } else if (pFlowTrack->TopologyMap().trackFtpcEast() ) {
739  detectorV1Wgt =pFlowEvent->V1FtpcEastDetctWgtG_Mix(k);
740  } else if (pFlowTrack->TopologyMap().trackFtpcWest() ) {
741  detectorV1Wgt =pFlowEvent->V1FtpcWestDetctWgtG_Mix(k);
742  }
743 
744 
745  double detectorV2Wgt = 1.;
746 
747  if (pFlowTrack->TopologyMap().hasHitInDetector(kTpcId) ||
748  (pFlowTrack->TopologyMap().data(0) == 0 &&
749  pFlowTrack->TopologyMap().data(1) == 0)) {
750  // hits in Tpc or TopologyMap not available
751  detectorV2Wgt =pFlowEvent->V2TPCDetctWgtG_Mix(k);
752  } else if (pFlowTrack->TopologyMap().trackFtpcEast() ) {
753  detectorV2Wgt =pFlowEvent->V2FtpcEastDetctWgtG_Mix(k);
754  } else if (pFlowTrack->TopologyMap().trackFtpcWest() ) {
755  detectorV2Wgt =pFlowEvent->V2FtpcWestDetctWgtG_Mix(k);
756  }
757 
758 
759  theCrosstermMix[k][j][pq] = evtGMix[k][j][pq] /
760  (1. + (phiWgtRaw*etaWgt*detectorV1Wgt/mMult[k][j]) *
761  (2.* histFull[k].histFullHar[j].mMixX1z[pIndex] * cos(phi * order) +
762  2.* histFull[k].histFullHar[j].mMixY1z[pIndex] * sin(phi * order) )
763  + (phiWgtRaw*ptWgt*detectorV2Wgt/mMult[k][j]) *
764  (2.* histFull[k].histFullHar[j].mMixX2z[qIndex] * cos(phi * order*2.) +
765  2.* histFull[k].histFullHar[j].mMixY2z[qIndex] * sin(phi * order*2.) ) );
766 
767  }
768 
769  // for writting out <G(p,q)>, the denominator
770 
771  DpqxMix[pIndex][qIndex] = theCrosstermMix[k][j][pq] * cos(order * phi)/
772  histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq];
773  DpqyMix[pIndex][qIndex] = theCrosstermMix[k][j][pq] * sin(order * phi)/
774  histFull[k].histFullHar[j].mCumulG0MixDenomRead[pq];
775 
776  }
777 
778 
779  for (int p=0; p< Flow::nCumulMixHar_pMax; p++){
780  DpxMix[p]=(1./(4.*r0Mix))*( DpqxMix[p][0]-DpqxMix[p][2]+DpqyMix[p][1]-DpqyMix[p][3]);
781  DpyMix[p]=(1./(4.*r0Mix))*( DpqyMix[p][0]-DpqyMix[p][2]+DpqxMix[p][3]-DpqxMix[p][1]);
782  }
783 
784 
785  cumuTempMix = (1./(4.*r0Mix))*(DpxMix[0]-DpyMix[2]-DpxMix[4]+DpyMix[6]);
786 
787  cumuTempMixFlip = cumuTempMix;
788 
789  if (eta < 0 && oddHar) {
790  cumuTempMixFlip *= -1.;
791  }
792 
793 
794  histFull[k].histFullHar[j].mHistCumulMix2D->
795  Fill(yOrEta, pt, cumuTempMix*profScale);
796  histFull[k].histFullHar[j].mHistCumulMixEta->
797  Fill(yOrEta, cumuTempMix*profScale);
798  histFull[k].histFullHar[j].mHistCumulMixPt->
799  Fill(pt, cumuTempMixFlip*profScale);
800 
801  histFull[k].mHistCumulMix->Fill(order,cumuTempMixFlip*profScale); //only order ==1 get filled anyway.
802 
803  }
804 
806 
807 
808  if (m_M==1) {
809  cumuTemp[0] = ((2.*Dp[1-1])-(0.5*Dp[2-1]))/r0Sq; // (B9, PG12)
810  cumuTemp[1] = ((-2.*Dp[1-1])+Dp[2-1])/(r0Sq*r0Sq);
811  } else if (m_M==2) {
812  cumuTemp[0] = ((4.*Dp[1-1])-(0.5*Dp[2-1]))/(r0Sq*r0Sq); // (B10)
813  cumuTemp[1] = ((-6.*Dp[1-1])+(1.5*Dp[2-1]))/(r0Sq*r0Sq*r0Sq);
814  }
815 
816  cumuTempFlip[0] = cumuTemp[0];
817  cumuTempFlip[1] = cumuTemp[1];
818  if (eta < 0 && oddHar) {
819  cumuTempFlip[0] *= -1.;
820  cumuTempFlip[1] *= -1.;
821  }
822 
823  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
824  histFull[k].histFullHar[j].mHistCumul2D[ord]->
825  Fill(yOrEta, pt, ((ord>0) ? cumuTemp[ord]*profScale : cumuTemp[ord]));
826  histFull[k].histFullHar[j].mHistCumulEta[ord]->
827  Fill(yOrEta, ((ord>0) ? cumuTemp[ord]*profScale : cumuTemp[ord]));
828  histFull[k].histFullHar[j].mHistCumulPt[ord]->
829  Fill(pt, ((ord>0) ? cumuTempFlip[ord]*profScale : cumuTempFlip[ord]));
830  }
831 
832  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++)
833  histFull[k].mHistCumul[ord]->Fill(order, ((ord>0) ? cumuTempFlip[ord]*profScale : cumuTempFlip[ord]) );
834  }
835  }
836  }
837  }
838 
839 }
840 
841 //-----------------------------------------------------------------------
842 
843 Int_t StFlowCumulantMaker::Finish() {
844 
845  TString* histTitle;
846 
847  TOrdCollection* XpqYpqDenomNames = new TOrdCollection(Flow::nSels*Flow::nHars);
848  TOrdCollection* savedHistNames = new TOrdCollection(Flow::nSels*Flow::nHars*Flow::nCumulDiffOrders);
849 
850  cout << endl << "##### Cumulant Maker:" << endl;
851  for (int k = 0; k < Flow::nSels; k++) {
852 
853  cout << "##### selection "<< k+1 <<" #### "<<endl;
854 
855  // integrated flow from cumulant
856  histFull[k].mHist_v = new TH1D*[Flow::nCumulDiffOrders];
857 
858  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
859  char theCumulOrder[2]; // if >10, need to use char*
860  sprintf(theCumulOrder,"%d",(ord+1)*2);
861 
862  histTitle = new TString("Flow_Cumul_v_Order");
863  histTitle->Append(*theCumulOrder);
864  histTitle->Append("_Sel");
865  *histTitle += k+1;
866  histFull[k].mHist_v[ord] =
867  new TH1D(*(histFull[k].mHistCumul[ord]->ProjectionX(histTitle->Data(),"e")));
868  histFull[k].mHist_v[ord]->SetTitle(histTitle->Data());
869  histFull[k].mHist_v[ord]->SetXTitle("harmonic");
870  histFull[k].mHist_v[ord]->SetYTitle("v (%)");
871  delete histTitle;
872  if (ord>0) histFull[k].mHist_v[ord]->Scale(1./profScale);
873  savedHistNames->AddLast(histFull[k].mHist_v[ord]);
874  }
875 
877  histTitle = new TString("Flow_CumulMix_v_Sel");
878  *histTitle +=k+1;
879  histFull[k].mHistMix_v =
880  new TH1D(*(histFull[k].mHistCumulMix->ProjectionX(histTitle->Data(),"e")));
881  histFull[k].mHistMix_v->SetTitle(histTitle->Data());
882  histFull[k].mHistMix_v->SetXTitle("place for v1 from mixed harmonic");
883  histFull[k].mHistMix_v->SetYTitle("v (%)");
884  delete histTitle;
885 
886  histFull[k].mHistMix_v->Scale(1./profScale);
887  savedHistNames->AddLast(histFull[k].mHistMix_v);
888 
890 
891  double meanIntegV[Flow::nHars]; // V**1
892  double meanIntegV2[Flow::nHars]; // V**2
893  double meanIntegV3[Flow::nHars]; // V**3
894  double meanIntegV4[Flow::nHars]; // V**4
895  double cumulInteg1[Flow::nHars]; // outside of harmonic loop
896  double cumulInteg2[Flow::nHars];
897  double cumulInteg3[Flow::nHars];
898 
899  double cumulIntegMix[Flow::nHars];
900  double meanIntegVMix[Flow::nHars];
901 
902  for (int j = 0; j < Flow::nHars; j++) {
903  meanIntegV[j] = 0.;
904  meanIntegV2[j] = 0.;
905  meanIntegV3[j] = 0.;
906  meanIntegV4[j] = 0.;
907  cumulInteg1[j] = 0.;
908  cumulInteg2[j] = 0.;
909  cumulInteg3[j] = 0.;
910  cumulIntegMix[j] = 0.;
911  meanIntegVMix[j] = 0.;
912  }
913 
914  for (int j = 0; j < Flow::nHars; j++) {
915 
916  double mAvMult = // average multiplicity
917  float(histFull[k].histFullHar[j].mMultSum)/
918  (float(histFull[k].histFullHar[j].mNEvent));
919 
920  histFull[k].histFullHar[j].mHistMultSum->
921  SetBinContent(1,double(histFull[k].histFullHar[j].mMultSum));
922  histFull[k].histFullHar[j].mHistNEvent->
923  SetBinContent(1,double(histFull[k].histFullHar[j].mNEvent));
924  histFull[k].histFullHar[j].mHistMeanWgtSqrSum->
925  SetBinContent(1,double(histFull[k].histFullHar[j].mMeanWgtSqrSum));
926 
927  double CpInteg[Flow::nCumulIntegOrders]; // Cp in (B4, PG6)
928 
929  for (int pq = 0; pq < Flow::nCumulIntegOrders; pq ++) CpInteg[pq] = 0.;
930  for (int pq = 0; pq < Flow::nCumulIntegOrders*Flow::nCumulInteg_qMax; pq++) {
931  int theCumulOrder = (pq/Flow::nCumulInteg_qMax) + 1; // like 1,2,3.
932  // not begining with 0. That is "p" in (B3, PG5)
933  //int qIndex = pq%(Flow::nCumulInteg_qMax); // like 0,1,..qMax-1
934  // begining with 0. "q" in (B3, PG5)
935 
936  histFull[k].histFullHar[j].mHistCumulIntegG0Sum[pq]-> // don't move it.
937  SetBinContent(1,histFull[k].histFullHar[j].mCumulIntegG0[pq]);
938 
939  //mCumulIntegG0 is the sum of IntegG0 until the code right below,
940  //where mCumulIntegG0 becomes <Gn(z)> :
941  histFull[k].histFullHar[j].mCumulIntegG0[pq] /=
942  float(histFull[k].histFullHar[j].mNEvent); // <Gn(z)> (PG 4)
943 
944  CpInteg[theCumulOrder-1] +=
945  (mAvMult*(::pow(histFull[k].histFullHar[j].mCumulIntegG0[pq], 1./mAvMult) -1.) /
946  float(Flow::nCumulInteg_qMax)); // (B3, PG6)
947 
948  }
949 
950  // add Xpq Ypq denominator to write out file list
951  for (int pq = 0; pq < Flow::nCumulDiffOrders*Flow::nCumulDiff_qMax; pq++) {
952  XpqYpqDenomNames->AddLast(histFull[k].histFullHar[j].mCumulG0Denom[pq]);
953  }
954 
955  cumulInteg1[j] = // (B5, PG7)
956  (3.*CpInteg[1-1] - (3./2.)*CpInteg[2-1] + (1./3.)*CpInteg[3-1]) / r0Sq;
957 
958  cumulInteg2[j] = ((-10.*CpInteg[1-1]) + (8.*CpInteg[2-1]) -
959  (2.*CpInteg[3-1])) / (r0Sq*r0Sq);
960 
961  cumulInteg3[j] = ( (18.*CpInteg[1-1]) - (18.*CpInteg[2-1]) + (6.*CpInteg[3-1]))
962  / (r0Sq*r0Sq*r0Sq);
963 
964  // now histograms for flow results:
965  histFull[k].histFullHar[j].mHist_v2D = new TH2D*[Flow::nCumulDiffOrders];
966  histFull[k].histFullHar[j].mHist_vEta = new TH1D*[Flow::nCumulDiffOrders];
967  histFull[k].histFullHar[j].mHist_vPt = new TH1D*[Flow::nCumulDiffOrders];
968 
969  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
970  char theCumulOrder[2]; // if >10, need to use char*
971  sprintf(theCumulOrder,"%d",(ord+1)*2);
972 
973  histTitle = new TString("Flow_Cumul_v2D_Order");
974  histTitle->Append(*theCumulOrder);
975  histTitle->Append("_Sel");
976  *histTitle += k+1;
977  histTitle->Append("_Har");
978  *histTitle += j+1;
979  histFull[k].histFullHar[j].mHist_v2D[ord] =
980  new TH2D(*(histFull[k].histFullHar[j].mHistCumul2D[ord]->
981  ProjectionXY(histTitle->Data(),"e")));
982  histFull[k].histFullHar[j].mHist_v2D[ord]->SetTitle(histTitle->Data());
983  histFull[k].histFullHar[j].mHist_v2D[ord]->SetXTitle((char*)xLabel.Data());
984  histFull[k].histFullHar[j].mHist_v2D[ord]->SetYTitle("Pt (GeV/c)");
985  histFull[k].histFullHar[j].mHist_v2D[ord]->SetZTitle("v (%)");
986  delete histTitle;
987 
988  if (ord>0) histFull[k].histFullHar[j].mHist_v2D[ord]->Scale(1./profScale);
989 
990  histTitle = new TString("Flow_Cumul_vEta_Order");
991  histTitle->Append(*theCumulOrder);
992  histTitle->Append("_Sel");
993  *histTitle += k+1;
994  histTitle->Append("_Har");
995  *histTitle += j+1;
996  histFull[k].histFullHar[j].mHist_vEta[ord] =
997  new TH1D(*(histFull[k].histFullHar[j].mHistCumulEta[ord]->
998  ProjectionX(histTitle->Data(),"e")));
999  histFull[k].histFullHar[j].mHist_vEta[ord]->SetTitle(histTitle->Data());
1000  histFull[k].histFullHar[j].mHist_vEta[ord]->SetXTitle((char*)xLabel.Data());
1001  histFull[k].histFullHar[j].mHist_vEta[ord]->SetYTitle("v (%)");
1002  delete histTitle;
1003 
1004  if (ord>0) histFull[k].histFullHar[j].mHist_vEta[ord]->Scale(1./profScale);
1005 
1006  histTitle = new TString("Flow_Cumul_vPt_Order");
1007  histTitle->Append(*theCumulOrder);
1008  histTitle->Append("_Sel");
1009  *histTitle += k+1;
1010  histTitle->Append("_Har");
1011  *histTitle += j+1;
1012  histFull[k].histFullHar[j].mHist_vPt[ord] =
1013  new TH1D(*(histFull[k].histFullHar[j].mHistCumulPt[ord]->
1014  ProjectionX(histTitle->Data(),"e")));
1015  histFull[k].histFullHar[j].mHist_vPt[ord]->SetTitle(histTitle->Data());
1016  histFull[k].histFullHar[j].mHist_vPt[ord]->SetXTitle("Pt (GeV/c)");
1017  histFull[k].histFullHar[j].mHist_vPt[ord]->SetYTitle("v (%)");
1018  delete histTitle;
1019 
1020  if (ord>0) histFull[k].histFullHar[j].mHist_vPt[ord]->Scale(1./profScale);
1021 
1022  }
1023 
1024  // new method, Eq. (PG8)
1025  meanIntegV[j] = ::sqrt(cumulInteg1[j]); // <v> 2-part, m=1
1026  meanIntegV2[j] = cumulInteg1[j]; // <v**2> 2-part, m=2
1027  meanIntegV3[j] = ::pow(-1.*cumulInteg2[j], 3./4.); // <v**3> 4-part, m=1
1028  meanIntegV4[j] = -1.*cumulInteg2[j]; // <v**4> 4-part, m=2
1029 
1030  if (meanIntegV2[j]<0.) cout<<" Sel"<<k+1<<", <V**2> less than zero ! v"
1031  <<j+1<<" from 2 particle correlation failed."<<endl;
1032  if (meanIntegV4[j]<0.) cout<<" Sel"<<k+1<<", <V**4> less than zero ! v"
1033  <<j+1<<" from 4 particle correlation failed."<<endl;
1034 
1035 
1036  if (m_M==1) { // Eq. (PG14)
1037  histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
1038  histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-1./(meanIntegV3[j]*perCent));// (34b)
1039  histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
1040  histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
1041  histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV[j]*perCent)); // (34a)
1042  histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-1./(meanIntegV3[j]*perCent)); // (34b)
1043  } else if (m_M==2) {
1044  histFull[k].histFullHar[j].mHist_v2D[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
1045  histFull[k].histFullHar[j].mHist_v2D[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
1046  histFull[k].histFullHar[j].mHist_vEta[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
1047  histFull[k].histFullHar[j].mHist_vEta[1]->Scale(-0.5/(meanIntegV4[j]*perCent) ); // (35b)
1048  histFull[k].histFullHar[j].mHist_vPt[0]->Scale(1./(meanIntegV2[j]*perCent)); // (35a)
1049  histFull[k].histFullHar[j].mHist_vPt[1]->Scale(-0.5/(meanIntegV4[j]*perCent)); // (35b)
1050  }
1051 
1052  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++) {
1053  savedHistNames->AddLast(histFull[k].histFullHar[j].mHist_v2D[ord]);
1054  savedHistNames->AddLast(histFull[k].histFullHar[j].mHist_vEta[ord]);
1055  savedHistNames->AddLast(histFull[k].histFullHar[j].mHist_vPt[ord]);
1056  }
1057  }
1058 
1059 
1061 
1062  for (int j = 0; j < Flow::nHars; j++) {
1063  if (j != 0) continue; //only j==0 makes sense for mixed har.
1064 
1065  double mAvMult = // average multiplicity
1066  float(histFull[k].histFullHar[j].mMultSum)/
1067  (float(histFull[k].histFullHar[j].mNEvent));
1068 
1069  double mAveMeanWgtSqr = // <w**2>
1070  float(histFull[k].histFullHar[j].mMeanWgtSqrSum)/
1071  (float(histFull[k].histFullHar[j].mNEvent));
1072 
1073  double CpqMix[Flow::nCumulMixHar_pMax][Flow::nCumulMixHar_qMax];//(30)
1074  for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax; pq++) {
1075  int pIndex = pq/Flow::nCumulMixHar_qMax;
1076  int qIndex = pq%Flow::nCumulMixHar_qMax;
1077 
1078  histFull[k].histFullHar[j].mHistCumulIntegG0MixSum[pq]->
1079  SetBinContent(1,histFull[k].histFullHar[j].mCumulIntegG0Mix[pq]);
1080  histFull[k].histFullHar[j].mCumulIntegG0Mix[pq] /= //<G(z1,z2)>
1081  float(histFull[k].histFullHar[j].mNEvent);
1082 
1083  CpqMix[pIndex][qIndex]=mAvMult*(pow(histFull[k].histFullHar[j].mCumulIntegG0Mix[pq], 1./mAvMult) -1.); //Mix (21)
1084 
1085  }
1086 
1087  double CpxMix[Flow::nCumulMixHar_pMax];
1088  double CpyMix[Flow::nCumulMixHar_pMax];
1089 
1090  for (int p =0; p<Flow::nCumulMixHar_pMax; p++){//Mix (31)
1091  CpxMix[p] = (1./(4.*r0Mix))*(CpqMix[p][0] - CpqMix[p][2]);
1092  CpyMix[p] = (1./(4.*r0Mix))*(CpqMix[p][3] - CpqMix[p][1]);
1093  }
1094 
1095  cumulIntegMix[j] = (1./(4.*r0Mix*r0Mix))*( // Mix (32)
1096  CpxMix[0]-CpyMix[1]-CpxMix[2]+CpyMix[3]
1097  +CpxMix[4]-CpyMix[5]-CpxMix[6]+CpyMix[7]);
1098 
1099 
1100  double tempMeanV = pow(-1.*cumulInteg2[1]*mAveMeanWgtSqr*mAveMeanWgtSqr,1./4.); // <wgt*v2>
1101  double tempMeanVMixSq = cumulIntegMix[j]/tempMeanV; //(24)
1102 
1103  if (tempMeanVMixSq>0.)
1104  meanIntegVMix[j] = sqrt(tempMeanVMixSq);//(24)
1105  else cout<<"### <wgt*v1>**2 = "<<tempMeanVMixSq<<" < 0. failed "<<endl;
1106 
1107 
1108  histTitle = new TString("Flow_CumulMix_v2D_Sel");
1109  *histTitle +=k+1;
1110  histTitle->Append("_Har");
1111  *histTitle +=j+1;
1112  histFull[k].histFullHar[j].mHistMix_v2D =
1113  new TH2D(*(histFull[k].histFullHar[j].mHistCumulMix2D->
1114  ProjectionXY(histTitle->Data(),"e")));
1115  histFull[k].histFullHar[j].mHistMix_v2D->SetTitle(histTitle->Data());
1116  histFull[k].histFullHar[j].mHistMix_v2D->SetXTitle((char*)xLabel.Data());
1117  histFull[k].histFullHar[j].mHistMix_v2D->SetYTitle("Pt (GeV)");
1118  histFull[k].histFullHar[j].mHistMix_v2D->SetZTitle("v (%)");
1119  delete histTitle;
1120 
1121  histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./profScale);
1122 
1123  histTitle = new TString("Flow_CumulMix_vEta_Sel");
1124  *histTitle +=k+1;
1125  histTitle->Append("_Har");
1126  *histTitle +=j+1;
1127  histFull[k].histFullHar[j].mHistMix_vEta =
1128  new TH1D(*(histFull[k].histFullHar[j].mHistCumulMixEta->
1129  ProjectionX(histTitle->Data(),"e")));
1130  histFull[k].histFullHar[j].mHistMix_vEta->SetTitle(histTitle->Data());
1131  histFull[k].histFullHar[j].mHistMix_vEta->SetXTitle((char*)xLabel.Data());
1132  histFull[k].histFullHar[j].mHistMix_vEta->SetYTitle("v (%)");
1133  delete histTitle;
1134 
1135  histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./profScale);
1136 
1137 
1138  histTitle = new TString("Flow_CumulMix_vPt_Sel");
1139  *histTitle +=k+1;
1140  histTitle->Append("_Har");
1141  *histTitle +=j+1;
1142  histFull[k].histFullHar[j].mHistMix_vPt =
1143  new TH1D(*(histFull[k].histFullHar[j].mHistCumulMixPt->
1144  ProjectionX(histTitle->Data(),"e")));
1145  histFull[k].histFullHar[j].mHistMix_vPt->SetTitle(histTitle->Data());
1146  histFull[k].histFullHar[j].mHistMix_vPt->SetXTitle("Pt (GeV)");
1147  histFull[k].histFullHar[j].mHistMix_vPt->SetYTitle("v (%)");
1148  delete histTitle;
1149 
1150  histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./profScale);
1151 
1152  histFull[k].histFullHar[j].mHistMix_v2D->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
1153  histFull[k].histFullHar[j].mHistMix_vEta->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
1154  histFull[k].histFullHar[j].mHistMix_vPt->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
1155  histFull[k].mHistMix_v->Scale(1./(tempMeanV*meanIntegVMix[j]*perCent));
1156 
1157  savedHistNames->AddLast(histFull[k].histFullHar[j].mHistMix_v2D);
1158  savedHistNames->AddLast(histFull[k].histFullHar[j].mHistMix_vEta);
1159  savedHistNames->AddLast(histFull[k].histFullHar[j].mHistMix_vPt);
1160  savedHistNames->AddLast(histFull[k].mHistMix_v);
1161 
1162  }
1164 
1165 
1166  if (m_M==1) {
1167 
1168  TH1D* histOfMeanIntegV = new TH1D(*(histFull[k].mHist_v[0]));
1169  histOfMeanIntegV->Reset();
1170 
1171  TH1D* histOfMeanIntegV3 = new TH1D(*(histFull[k].mHist_v[1]));
1172  histOfMeanIntegV3->Reset();
1173 
1174  for (int j = 1; j < Flow::nHars+1; j++) {
1175  histOfMeanIntegV->SetBinContent(j, 1./(meanIntegV[j-1]*perCent));
1176  histOfMeanIntegV->SetBinError(j,0.);
1177  histOfMeanIntegV3->SetBinContent(j, -1./(meanIntegV3[j-1]*perCent));
1178  histOfMeanIntegV3->SetBinError(j,0.);
1179  }
1180  histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV);
1181  histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV3);
1182 
1183  //force to be zero if the value is "nan"
1184  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++){
1185  for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
1186  if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
1187  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ) {
1188  histFull[k].mHist_v[ord]->SetBinContent(j,0.);
1189  histFull[k].mHist_v[ord]->SetBinError(j,0.);
1190  }
1191  }
1192  }
1193 
1194  for (int j = 1; j < Flow::nHars+1; j++) {
1195  cout << "##### 2-part v" << j << " = ("
1196  << histFull[k].mHist_v[0]->GetBinContent(j)
1197  <<" +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<" )"<<endl;
1198  cout << "##### 4-part v" << j << " = ("
1199  << histFull[k].mHist_v[1]->GetBinContent(j)
1200  <<" +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<" )"<<endl;
1201  }
1202 
1203  delete histOfMeanIntegV;
1204  delete histOfMeanIntegV3;
1205 
1206  } else if (m_M==2) {
1207 
1208  TH1D* histOfMeanIntegV2 = new TH1D(*(histFull[k].mHist_v[0]));
1209  histOfMeanIntegV2->Reset();
1210 
1211  TH1D* histOfMeanIntegV4 = new TH1D(*(histFull[k].mHist_v[1]));
1212  histOfMeanIntegV4->Reset();
1213 
1214  for (int j = 1; j < Flow::nHars+1; j++) {
1215  histOfMeanIntegV2->SetBinContent(j, 1./(meanIntegV2[j-1]*perCent));
1216  histOfMeanIntegV2->SetBinError(j,0.);
1217  histOfMeanIntegV4->SetBinContent(j, -0.5/(meanIntegV4[j-1]*perCent));
1218  histOfMeanIntegV4->SetBinError(j,0.);
1219  }
1220  histFull[k].mHist_v[0]->Multiply(histOfMeanIntegV2);
1221  histFull[k].mHist_v[1]->Multiply(histOfMeanIntegV4);
1222 
1223  //force to be zero if the value is "nan"
1224  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++){
1225  for (int j=0; j<histFull[k].mHist_v[ord]->GetNbinsX(); j++) {
1226  if ( !(histFull[k].mHist_v[ord]->GetBinContent(j) < FLT_MAX &&
1227  histFull[k].mHist_v[ord]->GetBinContent(j) > -1.*FLT_MAX) ){
1228  histFull[k].mHist_v[ord]->SetBinContent(j,0.);
1229  histFull[k].mHist_v[ord]->SetBinError(j,0.);
1230  }
1231  }
1232  }
1233 
1234  for (int j = 1; j < Flow::nHars+1; j++) {
1235  cout << "##### 2-part v" << j << " = ("
1236  << histFull[k].mHist_v[0]->GetBinContent(j)
1237  <<") +/- "<< histFull[k].mHist_v[0]->GetBinError(j)<<endl;
1238  cout << "##### 4-part v" << j << " = ("
1239  << histFull[k].mHist_v[1]->GetBinContent(j)
1240  <<") +/- "<< histFull[k].mHist_v[1]->GetBinError(j)<<endl;
1241  }
1242 
1243  delete histOfMeanIntegV2;
1244  delete histOfMeanIntegV4;
1245  }
1246 
1247  for (int ord = 0; ord < Flow::nCumulDiffOrders; ord++)
1248  savedHistNames->AddLast(histFull[k].mHist_v[ord]);
1249 
1250  }
1251 
1252  // GetHistList()->ls();
1253 
1254  // Write most histograms
1255  TFile histFile("flow.cumulant.root", "RECREATE");
1256  for (int k = 0; k < Flow::nSels; k++) {
1257  for (int j = 0; j < Flow::nHars; j++) {
1258  for (int pq = 0; pq < Flow::nCumulMixHar_pMax*Flow::nCumulMixHar_qMax ; pq++) {
1259  XpqYpqDenomNames->AddLast(histFull[k].histFullHar[j].mCumulG0MixDenom[pq]);
1260  }
1261  }
1262  }
1263 
1264  TVectorD* cumulConstants = new TVectorD(30);
1265  (*cumulConstants)(0)=double(Flow::nHars);
1266  (*cumulConstants)(1)=double(Flow::nSels);
1267  (*cumulConstants)(2)=double(Flow::nSubs);
1268  (*cumulConstants)(3)=double(Flow::nPhiBins);
1269  (*cumulConstants)(4)=double(Flow::nPhiBinsFtpc);
1270  (*cumulConstants)(5)=double(mNEtaBins);
1271  (*cumulConstants)(6)=double(nPtBinsPart);
1272  (*cumulConstants)(7)=double(Flow::nCumulIntegOrders);
1273  (*cumulConstants)(8)=double(Flow::nCumulInteg_qMax);
1274  (*cumulConstants)(9)=double(Flow::nCumulDiffOrders);
1275  (*cumulConstants)(10)=double(Flow::nCumulDiff_qMax);
1276  (*cumulConstants)(11)=double(Flow::nCumulMixHar_pMax);
1277  (*cumulConstants)(12)=double(Flow::nCumulMixHar_qMax);
1278  (*cumulConstants)(13)=r0;
1279  (*cumulConstants)(14)=m_M;
1280  (*cumulConstants)(15)=(strlen(pFlowSelect->PidPart()) != 0) ? 1 : 0;//1,pidflow
1281  (*cumulConstants)(16)=r0Mix;
1282  (*cumulConstants)(17)=double(profScale);
1283 
1284  cumulConstants->Write("CumulConstants",TObject::kOverwrite | TObject::kSingleKey);
1285  savedHistNames->AddLast(cumulConstants);
1286 
1287  TObjString* cumulMethodTag = new TObjString( "cumulNew" );
1288  cumulMethodTag->Write("CumulMethodTag",TObject::kOverwrite | TObject::kSingleKey);
1289  savedHistNames->AddLast(cumulMethodTag);
1290 
1291  savedHistNames->Write();
1292 
1293  histFile.Close();
1294 
1295  // write profile for <G>, the denominator of XpqYpq.
1296  TFile XpqYpqDenomNewFile("denominatorNew.root","RECREATE");
1297  XpqYpqDenomNames->Write();
1298  XpqYpqDenomNewFile.Close();
1299  delete XpqYpqDenomNames;
1300 
1301  delete pFlowSelect;
1302 
1303  return StMaker::Finish();
1304 }
1305 
1306 //-----------------------------------------------------------------------
1307 
1308 void StFlowCumulantMaker::SetHistoRanges(Bool_t ftpc_included) {
1309 
1310  if (ftpc_included) {
1311  mEtaMin = Flow::etaMin;
1312  mEtaMax = Flow::etaMax;
1313  mNEtaBins = Flow::nEtaBins;
1314  }
1315  else {
1316  mEtaMin = Flow::etaMinTpcOnly;
1317  mEtaMax = Flow::etaMaxTpcOnly;
1318  mNEtaBins = Flow::nEtaBinsTpcOnly;
1319  }
1320 
1321  return;
1322 }
1323 
1325 //
1326 // $Log: StFlowCumulantMaker.cxx,v $
1327 // Revision 1.22 2006/02/22 19:12:29 posk
1328 // Reduced number of output histograms.
1329 //
1330 // Revision 1.21 2005/10/27 17:53:18 aihong
1331 // changes made so that the cumulant method won't pick up the raw phi weight
1332 //
1333 // Revision 1.20 2004/12/17 15:50:09 aihong
1334 // check in v1{3} code
1335 //
1336 // Revision 1.19 2004/12/09 23:47:08 posk
1337 // Minor changes in code formatting.
1338 // Added hist for TPC primary dca to AnalysisMaker.
1339 //
1340 // Revision 1.18 2004/11/16 21:22:22 aihong
1341 // removed old cumulant method
1342 //
1343 // Revision 1.17 2004/08/24 20:22:39 oldi
1344 // Minor modifications to avoid compiler warnings.
1345 //
1346 // Revision 1.16 2003/09/02 17:58:10 perev
1347 // gcc 3.2 updates + WarnOff
1348 //
1349 // Revision 1.15 2003/07/07 21:58:19 posk
1350 // Made units of momentum GeV/c instead of GeV.
1351 //
1352 // Revision 1.14 2003/03/03 16:24:36 aihong
1353 // blow up 4-part cumulant by 1000 in order to let error bars calculated by ROOT
1354 //
1355 // Revision 1.13 2003/01/10 16:40:36 oldi
1356 // Several changes to comply with FTPC tracks:
1357 // - Switch to include/exclude FTPC tracks introduced.
1358 // The same switch changes the range of the eta histograms.
1359 // - Eta symmetry plots for FTPC tracks added and separated from TPC plots.
1360 // - PhiWgts and related histograms for FTPC tracks split in FarEast, East,
1361 // West, FarWest (depending on vertex.z()).
1362 // - Psi_Diff plots for 2 different selections and the first 2 harmonics added.
1363 // - Cut to exclude mu-events with no primary vertex introduced.
1364 // (This is possible for UPC events and FTPC tracks.)
1365 // - Global DCA cut for FTPC tracks added.
1366 // - Global DCA cuts for event plane selection separated for TPC and FTPC tracks.
1367 // - Charge cut for FTPC tracks added.
1368 //
1369 // Revision 1.12 2002/05/19 18:58:00 aihong
1370 // speed up cumulant
1371 //
1372 // Revision 1.11 2002/02/19 14:42:18 jeromel
1373 // Added float.h for Linux 7.2
1374 //
1375 // Revision 1.10 2002/02/02 01:10:12 posk
1376 // Added documentation
1377 //
1378 // Revision 1.9 2002/01/24 20:53:51 aihong
1379 // add histograms for combining results from parallel jobs
1380 //
1381 // Revision 1.8 2002/01/14 23:42:36 posk
1382 // Renamed ScalerProd histograms. Moved print commands to FlowMaker::Finish().
1383 //
1384 // Revision 1.7 2001/12/21 17:01:59 aihong
1385 // minor changes
1386 //
1387 // Revision 1.6 2001/12/18 19:27:27 posk
1388 // "proton" and "antiproton" replaced by "pr+" and "pr-".
1389 //
1390 // Revision 1.5 2001/12/11 22:04:01 posk
1391 // Four sets of phiWgt histograms.
1392 // StFlowMaker StFlowEvent::PhiWeight() changes.
1393 // Cumulant histogram names changed.
1394 //
1395 // Revision 1.4 2001/12/06 01:21:14 jeromel
1396 // Mandatory correction : Extraneous comma removed.
1397 //
1398 // Revision 1.3 2001/11/09 21:14:50 posk
1399 // Switched from CERNLIB to TMath. Using global dca instead of dca.
1400 //
1401 // Revision 1.2 2001/11/08 03:12:24 aihong
1402 // clean up redundant histograms
1403 //
1404 // Revision 1.1 2001/11/02 04:47:42 aihong
1405 // install cumulant maker
1406 //
kStOK
Definition: Stypes.h:40
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
Generated by   doxygen 1.8.5