StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StFlowDirectCumulantMaker.cxx
1 //
3 // $Id: StFlowDirectCumulantMaker.cxx,v 1.2 2010/06/10 16:33:53 posk Exp $
4 //
5 // Authors: Dhevan Gangadharan, UCLA, Dec 2009
6 //
8 //
9 // Description: Maker to analyze Flow using direct cumulants.
10 // It needs only one pass through the data.
11 // The outputs are 2D histograms for the needed terms, of events vs centrality and pt.
12 // The output file is flow.dirCumulant.root.
13 // The macro directCumulants_v2.C then converts these terms to cumulants
14 // and then v2(pt) for each centrality and min bias.
15 // cent=1 is the most central, cent=0 is the yield weighted average of all the others.
16 // The integrated v2 is also calculated by integrating over pt.
17 // Some documentation is at http://www.star.bnl.gov/protected/bulkcorr/posk/direct4part/index.html
18 //
20 
21 #include <Stiostream.h>
22 #include <stdlib.h>
23 #include <math.h>
24 #include "StMaker.h"
25 #include "StFlowDirectCumulantMaker.h"
26 #include "StFlowMaker/StFlowMaker.h"
27 #include "StFlowMaker/StFlowEvent.h"
28 #include "StFlowMaker/StFlowConstants.h"
29 #include "StFlowMaker/StFlowSelection.h"
30 #include "PhysicalConstants.h"
31 #include "SystemOfUnits.h"
32 #include "TFile.h"
33 #include "TString.h"
34 #include "TH2.h"
35 #include "StMessMgr.h"
36 #define PR(x) cout << "##### FlowAnalysis: " << (#x) << " = " << (x) << endl;
37 
38 ClassImp(StFlowDirectCumulantMaker)
39 
40 //-----------------------------------------------------------------------
41 StFlowDirectCumulantMaker::StFlowDirectCumulantMaker(const Char_t* name): StMaker(name),
42  MakerName(name) {
43  pFlowSelect = new StFlowSelection();
44 }
45 
46 StFlowDirectCumulantMaker::StFlowDirectCumulantMaker(const Char_t* name,
47  const StFlowSelection& flowSelect) :
48  StMaker(name), MakerName(name) {
49  pFlowSelect = new StFlowSelection(flowSelect); // copy constructor
50 }
51 
52 StFlowDirectCumulantMaker::~StFlowDirectCumulantMaker() {
53 }
54 
55 //-----------------------------------------------------------------------
56 Int_t StFlowDirectCumulantMaker::Make() {
57  // do events
58 
59  // Get a pointer to StFlowEvent
60  StFlowMaker* pFlowMaker = NULL;
61  pFlowMaker = (StFlowMaker*)GetMaker("Flow");
62  if (pFlowMaker) pFlowEvent = pFlowMaker->FlowEventPointer();
63  if (pFlowEvent) {
64  int runID = pFlowEvent->RunID();
65  //cout << "runID= " << runID << endl;
66 
67  if(runID < 8000000) { // not year 7, Mweights = 1
68  if (Mweights[0][0] != 1.) { // first event
69  gMessMgr->Info("##### FlowDirCumulant: Setting Multiplicity Weights to 1");
70  //cout << "Multiplicity Weights:" << endl;
71  for(int i=0; i<30; i++){
72  //cout << endl;
73  for(int j=0; j<1000; j++){
74  //cout << setprecision(3) << Mweights[i][j] << ", ";
75  Mweights[i][j]=1;
76  }
77  //cout << endl;
78  }
79  }
80  }
81 
82  FillParticleArrays(); // This fills the particle arrays from pFlowTrack
83  CalculateTerms();
84  Fill_Histograms(); // Fill Histograms and let things accumulate
85  } else {
86  gMessMgr->Info("##### FlowDirCumulant: FlowEvent pointer null");
87  return kStOK;
88  }
89 
90  if (Debug()) StMaker::PrintInfo();
91 
92  return kStOK;
93 }
94 
95 //-----------------------------------------------------------------------
96 Int_t StFlowDirectCumulantMaker::Init() {
97  // Book histograms
98  cout << endl <<"StFlowDirectCumulantMaker::Init()" << endl;
99 
100  GetMweights(); // Y7 mutiplicity weights; due to Y7 centrality bias
101 
102  TERM1=0;
103  TERM2=1;
104  TERM3=2;
105  TERM4=3;
106  TERM5=4;
107  TERM6=5;
108  TERM7=6;
109  TERM8=7;
110  TERM9=8;
111  TERM10=9;
112 
113  DIF=0;
114  INT=1;
115  COS=0;
116  SIN=1;
117 
118  TString *name;
119  const int ptbins = Flow::PTBINS - 1;
120 
121  for(int term=0; term<Flow::TERMS; term++){
122  for(int type=0; type<Flow::TYPES; type++){
123  for(int phase=0; phase<Flow::PHASES; phase++){
124  for(int species=0; species<Flow::SPECIES; species++){
125 
126  name = new TString("Term_");
127  *name += term+1;
128 
129  if(type==0) name->Append("_D");
130  else name->Append("_I");
131 
132  if(phase==0) name->Append("_cos_");
133  else name->Append("_sin_");
134 
135  *name += species;
136 
137  Term[term].Type[type].Phase[phase].Species[species].Sum_Singles = new TH2D(name->Data(),"Correlation Sum",Flow::nCents,.5,Flow::nCents+0.5,ptbins,0.05,ptbins/10.+0.05);
138  Term[term].Type[type].Phase[phase].Species[species].Sum_Singles->GetXaxis()->SetTitle("Centrality bin");
139  Term[term].Type[type].Phase[phase].Species[species].Sum_Singles->GetYaxis()->SetTitle("p_{t} bin/10");
140 
141  name->Append("_Sq");
142  Term[term].Type[type].Phase[phase].Species[species].Sum_Squares = new TH2D(name->Data(),"Squares Sum",Flow::nCents,.5,Flow::nCents+0.5,ptbins,0.05,ptbins/10.+0.05);
143  Term[term].Type[type].Phase[phase].Species[species].Sum_Squares->GetXaxis()->SetTitle("Centrality bin");
144  Term[term].Type[type].Phase[phase].Species[species].Sum_Squares->GetYaxis()->SetTitle("p_{t} bin/10");
145 
146  }
147  }
148  }
149  }
150 
151  for(int species=0; species<Flow::SPECIES; species++){
152  name = new TString("Cent_Pt_yields_");
153  *name += species;
154  Entries[species].cent_yields = new TH2D(name->Data(),"Pt_spectrum",Flow::nCents,.5,Flow::nCents+0.5,ptbins,0.05,ptbins/10.+0.05);
155  Entries[species].cent_yields->GetXaxis()->SetTitle("Centrality bin");
156  Entries[species].cent_yields->GetYaxis()->SetTitle("p_{t} bin");
157  }
158 
159  Event_counter = new TH1D("Event_counter","Event_counter",Flow::nCents,0.5,Flow::nCents+0.5);
160  Event_counterWeighted = new TH1D("Event_counterWeighted","Event_counterWeighted",Flow::nCents,0.5,Flow::nCents+0.5);
161  Event_counterWeighted->Sumw2();
162 
163  gMessMgr->SetLimit("##### FlowDirCumu", 2);
164  gMessMgr->Info("##### FlowDirCumu: $Id: StFlowDirectCumulantMaker.cxx,v 1.2 2010/06/10 16:33:53 posk Exp $");
165 
166  return StMaker::Init();
167 }
168 
169 //-----------------------------------------------------------------------
170 void StFlowDirectCumulantMaker::GetMweights() {
171  // Y7 mutiplicity weights; due to Y7 centrality bias
172 
173  TFile* WeightFile = new TFile("M_WeightsY7.root","READ");
174 
175  if(WeightFile->IsOpen()) {
176  for(int i=0; i<30; i++){ // vertex z
177  for(int j=0; j<1000; j++){ // ref mult
178 
179  Mweights[i][j] = ((TH2D*)WeightFile->Get("Weights"))->GetBinContent(i+1,j+1);
180 
181  }
182  }
183  WeightFile->Close();
184  } else {
185  gMessMgr->Info("##### FlowDirCumulant: No year 7 Multiplicity Weights");
186  for(int i=0; i<30; i++){
187  for(int j=0; j<1000; j++){
188  Mweights[i][j]=1;
189  }
190  }
191  }
192 }
193 
194 //----------------------------------------------------------------------
195 void StFlowDirectCumulantMaker::FillParticleArrays() {
196  // Fill _p primary quantities
197 
198  grefmult = pFlowEvent->MultEta();
199 
200  StThreeVectorF vertex = pFlowEvent->VertexPos();
201  zbin=0; // z_vertex bin used for multiplicity weights for the Y7 AuAu Centrality bias
202  for(int i=0; i<30; i++){
203  if( (vertex[2] > (2*i-30)) && (vertex[2] < (2*(i+1)-30)) ){
204  zbin=i+1;
205  break;
206  }
207  }
208 
209  // Clear previous event quantities
210  p_count=0;
211 
212  for(int i=0; i<Flow::PTBINS; i++) {
213  pt_bin_counter_p[i]=0;
214  }
215 
216  // Primary tracks
217  StFlowTrackCollection* pFlowTracks_primary = pFlowEvent->TrackCollection();
218  StFlowTrackIterator itr;
219 
221  // Primary track arrays
222  for (itr = pFlowTracks_primary->begin(); itr != pFlowTracks_primary->end(); itr++) {
223  StFlowTrack* pFlowTrack = *itr;
224  // cuts were made as StFlowTrack was filled
225 
226  phi_p[p_count] = pFlowTrack->Phi();
227  if(phi_p[p_count] < 0.) phi_p[p_count] += twopi;
228  id_p[p_count] = pFlowTrack->id();
229  pt_p[p_count] = pFlowTrack->Pt();
230 
231  weights_p[p_count] = 1.0;
232  if (pFlowEvent->PtWgt()) { // pt wgt
233  weights_p[p_count] *= (pt_p[p_count] < pFlowEvent->PtWgtSaturation()) ? pt_p[p_count] :
234  pFlowEvent->PtWgtSaturation(); // pt weighting going constant
235  }
236 
237  if(pt_p[p_count] <= .2 ) { ptbin_p[p_count] = 1; goto ptbin; }
238  else if(pt_p[p_count] <= .4 ) { ptbin_p[p_count] = 3; goto ptbin; }
239  else if(pt_p[p_count] <= .6 ) { ptbin_p[p_count] = 5; goto ptbin; }
240  else if(pt_p[p_count] <= .8 ) { ptbin_p[p_count] = 7; goto ptbin; }
241  else if(pt_p[p_count] <= 1.0) { ptbin_p[p_count] = 9; goto ptbin; }
242  else if(pt_p[p_count] <= 1.2) { ptbin_p[p_count] =11; goto ptbin; }
243  else if(pt_p[p_count] <= 1.4) { ptbin_p[p_count] =13; goto ptbin; }
244  else if(pt_p[p_count] <= 1.6) { ptbin_p[p_count] =15; goto ptbin; }
245  else if(pt_p[p_count] <= 1.8) { ptbin_p[p_count] =17; goto ptbin; }
246  else if(pt_p[p_count] <= 2.0) { ptbin_p[p_count] =19; goto ptbin; }
247  else if(pt_p[p_count] <= 2.2) { ptbin_p[p_count] =21; goto ptbin; }
248  else if(pt_p[p_count] <= 2.4) { ptbin_p[p_count] =23; goto ptbin; }
249  else if(pt_p[p_count] <= 2.6) { ptbin_p[p_count] =25; goto ptbin; }
250  else if(pt_p[p_count] <= 2.8) { ptbin_p[p_count] =27; goto ptbin; }
251  else if(pt_p[p_count] <= 3.2) { ptbin_p[p_count] =30; goto ptbin; }
252  else if(pt_p[p_count] <= 3.6) { ptbin_p[p_count] =34; goto ptbin; }
253  else if(pt_p[p_count] <= 4.0) { ptbin_p[p_count] =38; goto ptbin; }
254  else if(pt_p[p_count] <= 4.4) { ptbin_p[p_count] =42; goto ptbin; }
255  else if(pt_p[p_count] <= 4.8) { ptbin_p[p_count] =46; goto ptbin; }
256  else if(pt_p[p_count] <= 5.2) { ptbin_p[p_count] =50; goto ptbin; }
257  else if(pt_p[p_count] <= 5.6) { ptbin_p[p_count] =54; goto ptbin; }
258  else if(pt_p[p_count] <= 6.0) { ptbin_p[p_count] =58; goto ptbin; }
259  else { ptbin_p[p_count] =60; }
260 
261  ptbin:
262  pt_bin_counter_p[ptbin_p[p_count]]++;
263  p_count++;
264  if(p_count == Flow::MAXMULT) {
265  gMessMgr->Info("##### FlowDirCumulant: MAXMULT not big enough");
266  break;
267  }
268  }
269 
270 }
271 
272 //-----------------------------------------------------------------------
273 void StFlowDirectCumulantMaker::CalculateTerms() {
274  // Calculate terms
275 
276  if(p_count < 4) return; // Not a good event
277 
278  ClearTerms(); // Clears all the sum terms
279  CalculateSumTerms(kTRUE); // Fills the sum terms for the case of integrated flow
280  CalculateCorrelationsIntegrated(); // Calculates all the integrated correlations: 4-particle, 3-particle, 2-particle, and 1-particle of all needed types
281 
282  ClearTerms(); // Clears all the sum terms
283  CalculateSumTerms(kFALSE);// Fills the sum terms for the case of differential flow
284  CalculateCorrelationsDifferential(); // Calculates all the differential correlations: 4-particle, 3-particle, 2-particle,and 1-particle of all needed types
285 
286 }
287 
288 //----------------------------------------------------------------------
289 void StFlowDirectCumulantMaker::CalculateSumTerms(bool integrated) {
290  // Calculate component sum terms
291 
292  double phi_one = -1;
293  int ptbin_one = -1;
294  int id_one = -1;
295  double W_one = 1;
296  double W_two = 1;
297 
298  // species=0 == charged_hadron,
299  // All of the sum terms used below are components of the 4,3,2 particle correlations done with 2 loops.
300  // 1-particle correlations are of course still done with 1 loop.
301 
302  for(int species=0; species<Flow::SPECIES; species++){
304 
305  int max = p_count;
306  for(int one=0; one<max; one++){
308 
309  phi_one = phi_p[one];
310  ptbin_one = ptbin_p[one];
311  id_one = id_p[one];
312  if(integrated == kTRUE) W_one = weights_p[one];
313  else W_one = 1.;
314 
315  sum_1_cos_diff[species][ptbin_one] += cos(2*phi_one)*W_one;
316  sum_1_sin_diff[species][ptbin_one] += sin(2*phi_one)*W_one;
317  sum_2_cos_diff[species][ptbin_one] += cos(-2*phi_one)*W_one;
318  sum_2_sin_diff[species][ptbin_one] += sin(-2*phi_one)*W_one;
319 
320  sum_1_cos_full[species] += cos(2*phi_one)*W_one;
321  sum_1_sin_full[species] += sin(2*phi_one)*W_one;
322  sum_2_cos_full[species] += cos(-2*phi_one)*W_one;
323  sum_2_sin_full[species] += sin(-2*phi_one)*W_one;
324 
325  for(int two=0; two<p_count; two++){
327 
328  if(integrated == kTRUE || species==0) {if(id_p[two] == id_one) continue;}
329  W_two = weights_p[two];
330 
331  sum_3_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*W_one*W_two;
332  sum_3_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*W_one*W_two;
333 
334  sum_4_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*(phi_one-phi_p[two]))*W_one*W_two*W_one*W_two;
335  sum_4_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*sin(2*(phi_one-phi_p[two]))*W_one*W_two*W_one*W_two;
336 
337  sum_5_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*W_one*W_two*W_two;
338  sum_5_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*W_one*W_two*W_two;
339 
340  sum_6_cos_sin_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*W_one*W_two*W_two;
341  sum_6_sin_cos_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*W_one*W_two*W_two;
342 
343  sum_7_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*cos(2*phi_one)*W_one*W_two*W_one*W_two;
344 
345  sum_8_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*cos(2*phi_p[two])*cos(2*phi_p[two])*W_one*W_two*W_two*W_two;
346 
347  sum_9_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*sin(-2*phi_one)*W_one*W_two*W_one*W_two;
348 
349  sum_10_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one-phi_p[two]))*sin(2*phi_p[two])*sin(-2*phi_p[two])*W_one*W_two*W_two*W_two;
350 
351  sum_11_cos_diff[species][ptbin_one][ptbin_p[two]] += cos(2*(phi_one+phi_p[two]))*W_one*W_two;
352  sum_11_sin_diff[species][ptbin_one][ptbin_p[two]] += sin(2*(phi_one+phi_p[two]))*W_one*W_two;
353 
354  }
355  }
356 
357  int min_bin=0, max_bin=0;
358  for(int ptbin1=0; ptbin1<Flow::PTBINS; ptbin1++){
360  if(integrated == kTRUE) {min_bin=0; max_bin=Flow::PTBINS;}
361  else {min_bin = ptbin1; max_bin = ptbin1+1;}
362 
363  for(int ptbin2=min_bin; ptbin2<max_bin; ptbin2++){
365  if(integrated == kTRUE) {if(species==0) if(ptbin1==ptbin2) continue;}
366 
367  for(int ptbin3=0; ptbin3<Flow::PTBINS; ptbin3++){
369  if(species==0) if(ptbin1==ptbin3) continue;
370 
371  sum_3_cos[species][ptbin1] += sum_3_cos_diff[species][ptbin2][ptbin3];
372  sum_3_sin[species][ptbin1] += sum_3_sin_diff[species][ptbin2][ptbin3];
373 
374  sum_4_cos[species][ptbin1] += sum_4_cos_diff[species][ptbin2][ptbin3];
375  sum_4_sin[species][ptbin1] += sum_4_sin_diff[species][ptbin2][ptbin3];
376 
377  sum_5_cos[species][ptbin1] += sum_5_cos_diff[species][ptbin2][ptbin3];
378  sum_5_sin[species][ptbin1] += sum_5_sin_diff[species][ptbin2][ptbin3];
379 
380  sum_6_cos_sin[species][ptbin1] += sum_6_cos_sin_diff[species][ptbin2][ptbin3];
381  sum_6_sin_cos[species][ptbin1] += sum_6_sin_cos_diff[species][ptbin2][ptbin3];
382 
383  sum_7[species][ptbin1] += sum_7_diff[species][ptbin2][ptbin3];
384 
385  sum_8[species][ptbin1] += sum_8_diff[species][ptbin2][ptbin3];
386 
387  sum_9[species][ptbin1] += sum_9_diff[species][ptbin2][ptbin3];
388 
389  sum_10[species][ptbin1] += sum_10_diff[species][ptbin2][ptbin3];
390 
391  sum_11_cos[species][ptbin1] += sum_11_cos_diff[species][ptbin2][ptbin3];
392  sum_11_sin[species][ptbin1] += sum_11_sin_diff[species][ptbin2][ptbin3];
393 
394  }
395  }
396 
397  if(integrated == kTRUE){
398  sum_3_cos_reference[species][ptbin1] = sum_3_cos[species][ptbin1];
399  sum_3_sin_reference[species][ptbin1] = sum_3_sin[species][ptbin1];
400 
401  sum_11_cos_reference[species][ptbin1] = sum_11_cos[species][ptbin1];
402  sum_11_sin_reference[species][ptbin1] = sum_11_sin[species][ptbin1];
403 
404  if(species==0){
405  sum_1_cos[species][ptbin1] = sum_1_cos_full[species] - sum_1_cos_diff[species][ptbin1];
406  sum_1_sin[species][ptbin1] = sum_1_sin_full[species] - sum_1_sin_diff[species][ptbin1];
407  sum_2_cos[species][ptbin1] = sum_2_cos_full[species] - sum_2_cos_diff[species][ptbin1];
408  sum_2_sin[species][ptbin1] = sum_2_sin_full[species] - sum_2_sin_diff[species][ptbin1];
409 
410  sum_1_cos_reference[species][ptbin1] = sum_1_cos_full[species] - sum_1_cos_diff[species][ptbin1];
411  sum_1_sin_reference[species][ptbin1] = sum_1_sin_full[species] - sum_1_sin_diff[species][ptbin1];
412  sum_2_cos_reference[species][ptbin1] = sum_2_cos_full[species] - sum_2_cos_diff[species][ptbin1];
413  sum_2_sin_reference[species][ptbin1] = sum_2_sin_full[species] - sum_2_sin_diff[species][ptbin1];
414  }
415  else {
416  sum_1_cos[species][ptbin1] = sum_1_cos_full[species];
417  sum_1_sin[species][ptbin1] = sum_1_sin_full[species];
418  sum_2_cos[species][ptbin1] = sum_2_cos_full[species];
419  sum_2_sin[species][ptbin1] = sum_2_sin_full[species];
420 
421  sum_1_cos_reference[species][ptbin1] = sum_1_cos_full[species];
422  sum_1_sin_reference[species][ptbin1] = sum_1_sin_full[species];
423  sum_2_cos_reference[species][ptbin1] = sum_2_cos_full[species];
424  sum_2_sin_reference[species][ptbin1] = sum_2_sin_full[species];
425  }
426  }
427  else {
428  sum_1_cos[species][ptbin1] = sum_1_cos_diff[species][ptbin1];
429  sum_1_sin[species][ptbin1] = sum_1_sin_diff[species][ptbin1];
430  sum_2_cos[species][ptbin1] = sum_2_cos_diff[species][ptbin1];
431  sum_2_sin[species][ptbin1] = sum_2_sin_diff[species][ptbin1];
432  }
433 
434  }
435 
436  }
437 
438 }
439 
440 //----------------------------------------------------------------------
441 void StFlowDirectCumulantMaker::Combinations(Bool_t integrated, int species, int counter) {
442  // set the number of combinatinos to ref_combos_4,3,2,1,0
443 
444  ref_combos_4 = 0;
445  ref_combos_3 = 0;
446  ref_combos_2 = 0;
447  ref_combos_1 = 0;
448  ref_combos_0 = 0;
449 
450  if(integrated==kTRUE) {
451  if(species==0) {
452  ref_combos_4 = (p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1)*(p_count-pt_bin_counter_p[counter]-2);
453  ref_combos_4 *=(p_count-pt_bin_counter_p[counter]-3);
454  ref_combos_3 = (p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1)*(p_count-pt_bin_counter_p[counter]-2);
455  ref_combos_2 = (p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1);
456  ref_combos_1 = (p_count-pt_bin_counter_p[counter]);
457  ref_combos_0 = (p_count-pt_bin_counter_p[counter]);
458  }
459  }
460  else {
461  if(species==0) {
462  ref_combos_4 = pt_bin_counter_p[counter]*(p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1);
463  ref_combos_4 *= (p_count-pt_bin_counter_p[counter]-2);
464  ref_combos_3 = pt_bin_counter_p[counter]*(p_count-pt_bin_counter_p[counter])*(p_count-pt_bin_counter_p[counter]-1);
465  ref_combos_2 = pt_bin_counter_p[counter]*(p_count-pt_bin_counter_p[counter]);
466  ref_combos_1 = pt_bin_counter_p[counter];
467  ref_combos_0 = (p_count-pt_bin_counter_p[counter]);
468  }
469  }
470 
471 }
472 
473 //----------------------------------------------------------------------
474 void StFlowDirectCumulantMaker::CalculateCorrelationsDifferential() {
475  // Calculate terms for a fixed pt
476 
477  // species=0 == charged_hadron
478  // This part puts the component sum terms together to form the relevant correlations
479 
480  for(int species=0; species<Flow::SPECIES; species++){
481 
482  for(int ptbin=0; ptbin<Flow::PTBINS; ptbin++){
483 
484  Combinations(kFALSE,species,ptbin);
485 
486  if(ref_combos_4 == 0) continue;
487  if(ref_combos_3 == 0) continue;
488  if(ref_combos_2 == 0) continue;
489  if(ref_combos_1 == 0) continue;
490  if(ref_combos_0 == 0) continue;
491 
492  Cumulant_Terms[TERM1][DIF][COS][species][ptbin] = sum_3_cos[species][ptbin]*sum_3_cos_reference[species][ptbin];
493  Cumulant_Terms[TERM1][DIF][COS][species][ptbin] -= 2*(sum_5_cos[species][ptbin]*sum_2_cos_reference[species][ptbin] - sum_8[species][ptbin]);
494  Cumulant_Terms[TERM1][DIF][COS][species][ptbin] += 2*(sum_6_cos_sin[species][ptbin]*sum_2_sin_reference[species][ptbin] - sum_10[species][ptbin]);
495  Cumulant_Terms[TERM1][DIF][COS][species][ptbin] /= double(ref_combos_4);
496 
497  Cumulant_Terms[TERM2][DIF][COS][species][ptbin] = sum_3_cos[species][ptbin]/double(ref_combos_2);
498  Cumulant_Terms[TERM2][DIF][SIN][species][ptbin] = sum_3_sin[species][ptbin]/double(ref_combos_2);
499 
500  Cumulant_Terms[TERM4][DIF][COS][species][ptbin] = sum_1_cos[species][ptbin]/double(ref_combos_1);
501  Cumulant_Terms[TERM4][DIF][SIN][species][ptbin] = sum_1_sin[species][ptbin]/double(ref_combos_1);
502 
503  Cumulant_Terms[TERM7][DIF][COS][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_2_cos_reference[species][ptbin] - sum_5_cos[species][ptbin]));
504  Cumulant_Terms[TERM7][DIF][COS][species][ptbin] -= (sum_3_sin[species][ptbin]*sum_2_sin_reference[species][ptbin] + sum_5_sin[species][ptbin]);
505  Cumulant_Terms[TERM7][DIF][COS][species][ptbin] /= double(ref_combos_3);
506 
507  Cumulant_Terms[TERM7][DIF][SIN][species][ptbin] = (sum_3_sin[species][ptbin]*sum_2_cos_reference[species][ptbin]) + (sum_3_cos[species][ptbin]*sum_2_sin_reference[species][ptbin]);
508  Cumulant_Terms[TERM7][DIF][SIN][species][ptbin] -= (sum_6_sin_cos[species][ptbin] - sum_6_cos_sin[species][ptbin]);
509  Cumulant_Terms[TERM7][DIF][SIN][species][ptbin] /= double(ref_combos_3);
510 
511  Cumulant_Terms[TERM9][DIF][COS][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_1_cos_reference[species][ptbin] - sum_5_cos[species][ptbin]));
512  Cumulant_Terms[TERM9][DIF][COS][species][ptbin] -= (sum_3_sin[species][ptbin]*sum_1_sin_reference[species][ptbin] - sum_5_sin[species][ptbin]);
513  Cumulant_Terms[TERM9][DIF][COS][species][ptbin] /= double(ref_combos_3);
514 
515  Cumulant_Terms[TERM9][DIF][SIN][species][ptbin] = (sum_3_sin[species][ptbin]*sum_1_cos_reference[species][ptbin]) + (sum_3_cos[species][ptbin]*sum_1_sin_reference[species][ptbin]);
516  Cumulant_Terms[TERM9][DIF][SIN][species][ptbin] -= (sum_6_sin_cos[species][ptbin] + sum_6_cos_sin[species][ptbin]);
517  Cumulant_Terms[TERM9][DIF][SIN][species][ptbin] /= double(ref_combos_3);
518 
519  Cumulant_Terms[TERM10][DIF][COS][species][ptbin] = sum_11_cos[species][ptbin]/double(ref_combos_2);
520  Cumulant_Terms[TERM10][DIF][SIN][species][ptbin] = sum_11_sin[species][ptbin]/double(ref_combos_2);
521 
522  }
523  }
524 }
525 
526 //----------------------------------------------------------------------
527 void StFlowDirectCumulantMaker::CalculateCorrelationsIntegrated() {
528  // Calculate terms missing a pt to avoid autocorrelations
529 
530  // species=0 == charged_hadron
531  // Similar to CalculateCorrelationsDifferential but for the integrated flow case
532 
533  for(int species=0; species<Flow::SPECIES; species++){
534 
535  for(int ptbin=0; ptbin<Flow::PTBINS; ptbin++){
536 
537  Combinations(kTRUE,species,ptbin);
538 
539  if(ref_combos_4 == 0) continue;
540  if(ref_combos_3 == 0) continue;
541  if(ref_combos_2 == 0) continue;
542  if(ref_combos_1 == 0) continue;
543  if(ref_combos_0 == 0) continue;
544 
545  Cumulant_Terms[TERM1][INT][COS][species][ptbin] = sum_3_cos[species][ptbin]*sum_3_cos[species][ptbin] - 2*sum_4_cos[species][ptbin];
546  Cumulant_Terms[TERM1][INT][COS][species][ptbin] -= 4*((sum_5_cos[species][ptbin]*sum_2_cos[species][ptbin] - sum_8[species][ptbin] - sum_7[species][ptbin]) - (sum_6_cos_sin[species][ptbin]*sum_2_sin[species][ptbin] - sum_10[species][ptbin] - sum_9[species][ptbin]));
547  Cumulant_Terms[TERM1][INT][COS][species][ptbin] /= double(ref_combos_4);
548 
549  Cumulant_Terms[TERM2][INT][COS][species][ptbin] = sum_3_cos[species][ptbin]/double(ref_combos_2);
550  Cumulant_Terms[TERM2][INT][SIN][species][ptbin] = sum_3_sin[species][ptbin]/double(ref_combos_2);
551  Cumulant_Terms[TERM3][INT][COS][species][ptbin] = sum_3_cos[species][ptbin]/double(ref_combos_2);
552  Cumulant_Terms[TERM3][INT][SIN][species][ptbin] = sum_3_sin[species][ptbin]/double(ref_combos_2);
553 
554  Cumulant_Terms[TERM4][INT][COS][species][ptbin] = sum_1_cos[species][ptbin]/double(ref_combos_1);
555  Cumulant_Terms[TERM4][INT][SIN][species][ptbin] = sum_1_sin[species][ptbin]/double(ref_combos_1);
556 
557  Cumulant_Terms[TERM5][INT][COS][species][ptbin] = (sum_3_cos[species][ptbin]*sum_2_cos[species][ptbin] - sum_3_sin[species][ptbin]*sum_2_sin[species][ptbin]);
558  Cumulant_Terms[TERM5][INT][COS][species][ptbin] -= 2*(sum_5_cos[species][ptbin]);
559  Cumulant_Terms[TERM5][INT][COS][species][ptbin] /= double(ref_combos_3);
560  Cumulant_Terms[TERM5][INT][SIN][species][ptbin] = (sum_3_cos[species][ptbin]*sum_2_sin[species][ptbin] + sum_3_sin[species][ptbin]*sum_2_cos[species][ptbin]);
561 
562  Cumulant_Terms[TERM5][INT][SIN][species][ptbin] += 2*(sum_6_cos_sin[species][ptbin]);
563  Cumulant_Terms[TERM5][INT][SIN][species][ptbin] /= double(ref_combos_3);
564 
565  Cumulant_Terms[TERM6][INT][COS][species][ptbin] = sum_1_cos[species][ptbin]/double(ref_combos_1);
566  Cumulant_Terms[TERM6][INT][SIN][species][ptbin] = sum_1_sin[species][ptbin]/double(ref_combos_1);
567 
568  Cumulant_Terms[TERM7][INT][COS][species][ptbin] = Cumulant_Terms[TERM5][INT][COS][species][ptbin];
569  Cumulant_Terms[TERM7][INT][SIN][species][ptbin] = Cumulant_Terms[TERM5][INT][SIN][species][ptbin];
570 
571  Cumulant_Terms[TERM8][INT][COS][species][ptbin] = sum_2_cos[species][ptbin]/double(ref_combos_1);
572  Cumulant_Terms[TERM8][INT][SIN][species][ptbin] = sum_2_sin[species][ptbin]/double(ref_combos_1);
573 
574  Cumulant_Terms[TERM9][INT][COS][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_1_cos[species][ptbin] - 2*sum_5_cos[species][ptbin]));
575  Cumulant_Terms[TERM9][INT][COS][species][ptbin] -= (sum_3_sin[species][ptbin]*sum_1_sin[species][ptbin]);
576  Cumulant_Terms[TERM9][INT][COS][species][ptbin] /= double(ref_combos_3);
577  Cumulant_Terms[TERM9][INT][SIN][species][ptbin] = ((sum_3_cos[species][ptbin]*sum_1_sin[species][ptbin] - 2*sum_6_cos_sin[species][ptbin]));
578  Cumulant_Terms[TERM9][INT][SIN][species][ptbin] += (sum_3_sin[species][ptbin]*sum_1_cos[species][ptbin]);
579  Cumulant_Terms[TERM9][INT][SIN][species][ptbin] /= double(ref_combos_3);
580 
581  Cumulant_Terms[TERM10][INT][COS][species][ptbin] = sum_11_cos[species][ptbin]/double(ref_combos_2);
582  Cumulant_Terms[TERM10][INT][SIN][species][ptbin] = sum_11_sin[species][ptbin]/double(ref_combos_2);
583 
584  }
585  }
586 
587 }
588 
589 //----------------------------------------------------------------------
590 void StFlowDirectCumulantMaker::Fill_Histograms() {
591  // Fill term histograms
592 
593  Event_counter->Fill(pFlowEvent->Centrality());
594  if (Mweights[zbin-1][grefmult-1] < 12.2) {
595  Event_counterWeighted->Fill(pFlowEvent->Centrality(), Mweights[zbin-1][grefmult-1]);
596  }
597 
598  for(int ptbin=0; ptbin<Flow::PTBINS; ptbin++){
599  for(int species=0; species<Flow::SPECIES; species++){
600 
601  double the_pt = double(double(ptbin)/10.); // GeV = bin/10.
602  if(ptbin == Flow::PTBINS-1) Combinations(kTRUE,species,ptbin);
603  else Combinations(kFALSE,species,ptbin);
604 
605  if(ref_combos_4 == 0) continue;
606  if(ref_combos_3 == 0) continue;
607  if(ref_combos_2 == 0) continue;
608  if(ref_combos_1 == 0) continue;
609  if(ref_combos_0 == 0) continue;
610 
611  if(Mweights[zbin-1][grefmult-1] > 12.2) continue;
612  double PF = ref_combos_1*Mweights[zbin-1][grefmult-1]; // prefactor for Y7 centrality bias correction
613 
614  for(int phase=0; phase<Flow::PHASES; phase++){
615  for(int type=0; type<Flow::TYPES; type++){
616  for(int term=0; term<Flow::TERMS; term++){
617  if(fabs(Cumulant_Terms[term][type][phase][species][ptbin]) < 1.) {
618  Term[term].Type[type].Phase[phase].Species[species].Sum_Singles->Fill(pFlowEvent->Centrality(),the_pt, PF*Cumulant_Terms[term][type][phase][species][ptbin]);
619 
620  Term[term].Type[type].Phase[phase].Species[species].Sum_Squares->Fill(pFlowEvent->Centrality(),the_pt, PF*pow(Cumulant_Terms[term][type][phase][species][ptbin],2));
621  }
622 
623 
624  }
625  }
626  }
627  Entries[species].cent_yields->Fill(pFlowEvent->Centrality(), the_pt, ref_combos_1*Mweights[zbin-1][grefmult-1]);
628  }
629  }
630 
631 }
632 
633 //----------------------------------------------------------------------
634 void StFlowDirectCumulantMaker::ClearTerms() {
635  // All the component sum terms need to be cleared twice per event: 1 for the Differential calculation, 1 for the Integrated calculation.
636 
637  for(int species=0; species<Flow::SPECIES; species++){
638 
639  sum_1_cos_full[species]=0;
640  sum_1_sin_full[species]=0;
641  sum_2_cos_full[species]=0;
642  sum_2_sin_full[species]=0;
643 
644  for(int ptbin1=0; ptbin1<Flow::PTBINS; ptbin1++){
645 
646  for(int ptbin2=0; ptbin2<Flow::PTBINS; ptbin2++){
647 
648  sum_1_cos_diff[species][ptbin1]=0, sum_1_sin_diff[species][ptbin1]=0;
649  sum_2_cos_diff[species][ptbin1]=0, sum_2_sin_diff[species][ptbin1]=0;
650  sum_3_cos_diff[species][ptbin1][ptbin2]=0, sum_3_sin_diff[species][ptbin1][ptbin2]=0;
651  sum_4_cos_diff[species][ptbin1][ptbin2]=0, sum_4_sin_diff[species][ptbin1][ptbin2]=0;
652  sum_5_cos_diff[species][ptbin1][ptbin2]=0, sum_5_sin_diff[species][ptbin1][ptbin2]=0;
653  sum_6_cos_sin_diff[species][ptbin1][ptbin2]=0, sum_6_sin_cos_diff[species][ptbin1][ptbin2]=0;
654  sum_7_diff[species][ptbin1][ptbin2]=0, sum_8_diff[species][ptbin1][ptbin2]=0, sum_9_diff[species][ptbin1][ptbin2]=0, sum_10_diff[species][ptbin1][ptbin2]=0;
655  sum_11_cos_diff[species][ptbin1][ptbin2]=0, sum_11_sin_diff[species][ptbin1][ptbin2]=0;
656 
657  sum_1_cos[species][ptbin1]=0, sum_1_sin[species][ptbin1]=0;
658  sum_2_cos[species][ptbin1]=0, sum_2_sin[species][ptbin1]=0;
659  sum_3_cos[species][ptbin1]=0, sum_3_sin[species][ptbin1]=0;
660  sum_4_cos[species][ptbin1]=0, sum_4_sin[species][ptbin1]=0;
661  sum_5_cos[species][ptbin1]=0, sum_5_sin[species][ptbin1]=0;
662  sum_6_cos_sin[species][ptbin1]=0, sum_6_sin_cos[species][ptbin1]=0;
663  sum_7[species][ptbin1]=0, sum_8[species][ptbin1]=0, sum_9[species][ptbin1]=0, sum_10[species][ptbin1]=0;
664  sum_11_cos[species][ptbin1]=0, sum_11_sin[species][ptbin1]=0;
665 
666  }
667 
668  }
669  }
670 
671 }
672 
673 //-----------------------------------------------------------------------
674 Int_t StFlowDirectCumulantMaker::Finish() {
675  // Wtite out file of term histograms
676 
677  cout << endl << "##### Direct Cumulant Maker:" << endl;
678 
679  // Write all terms
680  const char *flowname = "flow.dirCumulant.root";
681  TFile histFile(flowname, "RECREATE");
682 
683  GetHistList()->Write();
684  //GetHistList()->ls();
685 
686  histFile.Close();
687 
688  delete pFlowSelect;
689 
690  return StMaker::Finish();
691 }
692 
693 //-----------------------------------------------------------------------
kStOK
Definition: Stypes.h:40
StFlowDirectCumulantMaker::StFlowDirectCumulantMaker
StFlowDirectCumulantMaker(const Char_t *name="FlowDirectCumulant")
Default constructor.
Definition: StFlowDirectCumulantMaker.cxx:41
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
Generated by   doxygen 1.8.5